流动湍流及其控制的数值模拟研究

1、流动湍流及其控制的数值模拟流动湍流及其控制的数值模拟唐国宁20051130内容提要内容提要流动湍流研究的简要回顾流动湍流研究的简要回顾湍流控制的一般方法湍流控制的一般方法直接数值模拟的常用方法直接数值模拟的常用方法流动湍流控制的数值模拟结果流动湍流控制的数值模拟结果(一)流动湍流研究的简要回顾流动湍流研究的简要回顾 层流层流 流体质点的运动轨迹光滑而有规则流体质点的运动轨迹光滑而有规则 湍流湍流 是一种有结构、不规则、多尺度的随是一种有结构、不规则、多尺度的随机流动机流动 最早由雷诺圆管实验最早由雷诺圆管实验, 人们认识到湍流的随人们认识到湍流的随机性机性 1894年年,O.Reynolds

2、首次将湍流流动分解首次将湍流流动分解为为平均量平均量与与脉动量脉动量之和之和 ，并由此导出湍流并由此导出湍流的统计方程的统计方程 ,uuuppp 21iiijijijjuuupuftxxxx 0iiux 雷诺平均方程雷诺平均方程 脉动方程脉动方程 21ijiiijijijjjuuuuupuftxxxxx 0iiux 21()iiiijjiijjjijjjuuuupuuu uu utxxxxxx 0iiux雷诺应力雷诺应力 封闭湍流统计方程的方法叫封闭湍流统计方程的方法叫湍流模式湍流模式 代数涡黏模式、代数涡黏模式、 涡黏模式（包括线性涡黏模式（包括线性和非线性），代数应力模式和雷诺应力模和非线

3、性），代数应力模式和雷诺应力模式（或称二价矩模式）式（或称二价矩模式） 20世纪初，一些杰出的流体力学家相继世纪初，一些杰出的流体力学家相继对对涡黏系数涡黏系数提出各种提出各种流体力学模型流体力学模型，如，如Taylor(1921年年)的的涡模型涡模型，Prandtl(1925年年)的的混合长度模型混合长度模型和和Vonkaman(1930年年)相似相似模型等模型等. 考虑连续介质不规则运动的特点运动考虑连续介质不规则运动的特点运动的的多尺度性多尺度性。导致产生描述多尺度的。导致产生描述多尺度的谱概谱概念和谱方法念和谱方法，并最终产生了，并最终产生了Kolmogorov的的局部各向同性的局部各

4、向同性的通用谱通用谱（即（即-5/3谱）谱） 近近50年来，大量的物理实验和数值模拟使年来，大量的物理实验和数值模拟使得人们对得人们对湍流本质湍流本质的认识越来越深入，人的认识越来越深入，人们注意到湍流并不是完全随机的，而是有们注意到湍流并不是完全随机的，而是有某种相干结构，某种相干结构，湍流湍流是一种有是一种有结构的不规结构的不规则多尺度流动则多尺度流动 实验和数值模拟表明：实验和数值模拟表明：无论是简单还是复无论是简单还是复杂湍流，都存在一定的大涡结构杂湍流，都存在一定的大涡结构。湍流中。湍流中大涡拟序结构对湍流生成和发展有主宰作大涡拟序结构对湍流生成和发展有主宰作用，因此抑制或消除大涡结

5、构可能抑制整用，因此抑制或消除大涡结构可能抑制整体的湍流强度，甚至使流动层流化，体的湍流强度，甚至使流动层流化，这是这是近代湍流和降噪的主要思想近代湍流和降噪的主要思想 湍流的有结构、不规则、多尺度流动导致湍流的有结构、不规则、多尺度流动导致了研究湍流的困难了研究湍流的困难 最好的湍流模型就是最好的湍流模型就是N-S方程方程 湍流的直接数值计算也存在局限性。湍流的直接数值计算也存在局限性。为了为了能模拟最小尺度的涡，计算网格的分辨率能模拟最小尺度的涡，计算网格的分辨率要足够高要足够高,而计算区域的尺寸应足以容纳最而计算区域的尺寸应足以容纳最大的涡，因此三维计算域的网格点数具有大的涡，因此三维计

6、算域的网格点数具有 ，例如当，例如当Re=5000时时,三维计算域三维计算域的网格点数至少有的网格点数至少有 2.1* 10*8, 由于计算机由于计算机内存的限制，内存的限制，目前直接数值模拟只能研究目前直接数值模拟只能研究低雷诺数湍流低雷诺数湍流Re5000.9/4xyzeN N NR 从从20世纪以来，科学家虽然对湍流的不规世纪以来，科学家虽然对湍流的不规则性进行了深入的理论探讨，则性进行了深入的理论探讨，以以Kolmogorov为代表的莫斯科学派和以为代表的莫斯科学派和以Taylor、batchelor为代表的剑桥学派对湍为代表的剑桥学派对湍流有杰出的贡献流有杰出的贡献，但是认真考察一百

7、多年，但是认真考察一百多年的湍流文献，人们发现推动湍流研究的理的湍流文献，人们发现推动湍流研究的理论屈指可数。论屈指可数。可以称做奠基性的理论只有可以称做奠基性的理论只有Kolmogorov的局部各向同性湍流理论的局部各向同性湍流理论。（二）湍流控制的一般方法湍流控制的一般方法在流体力学领域，在流体力学领域，对流体流动的控制就是对流体流动的控制就是要使壁湍流和无剪切流动发生有益的变化要使壁湍流和无剪切流动发生有益的变化，例如使例如使转捩延迟或提前；抑制或增强湍流；转捩延迟或提前；抑制或增强湍流；防止或促进分离防止或促进分离等。产生有用的结果是：等。产生有用的结果是：使阻力减少、升力增加、混合增

8、强和抑制使阻力减少、升力增加、混合增强和抑制流动引起的噪声流动引起的噪声。迄今人们提出了许多湍。迄今人们提出了许多湍流控制方法，有主动控制和被动控制两大流控制方法，有主动控制和被动控制两大类，主要控制方法有类，主要控制方法有： 壁面的抽吸或注入。壁面的抽吸或注入。 添加聚合物、表面话化剂、纤维、微型水添加聚合物、表面话化剂、纤维、微型水泡。泡。 用磁场控制流场。用磁场控制流场。 电磁力控制流场。电磁力控制流场。 化学反应改变流场温度或流体密度。化学反应改变流场温度或流体密度。 改变流场几何参数改变流场几何参数如在壁面开小尺度的如在壁面开小尺度的纵向糟或肋条，使壁面振荡等。纵向糟或肋条，使壁面振

9、荡等。 多孔壁。多孔壁。 在壁面中嵌入物体在壁面中嵌入物体。D型柱体绕流L型槽形肋条控制壁面的吹吸（三）直接数值模拟常用方法直接数值模拟常用方法湍流的数值模拟方法有两种：湍流的数值模拟方法有两种： 大涡模拟大涡模拟 直接数值模拟直接数值模拟. 大涡模拟的基本思想大涡模拟的基本思想是把包括脉动在内的湍流的是把包括脉动在内的湍流的瞬时运动通过某种滤波方法分解为瞬时运动通过某种滤波方法分解为大尺度运动大尺度运动和和小尺度运动小尺度运动两部分。两部分。 大尺度运动（大涡）强烈依赖于初始条件大尺度运动（大涡）强烈依赖于初始条件，其形，其形态与强度因平均流动的不同而不同。态与强度因平均流动的不同而不同。

10、大涡与平均运动之间有强烈的相互作用大涡与平均运动之间有强烈的相互作用，它直接，它直接由平均运动提供能量，反过来又对平均运动有强由平均运动提供能量，反过来又对平均运动有强烈的影响，大部分的质量、动量和能量的输运是烈的影响，大部分的质量、动量和能量的输运是由大涡运动引起。由大涡运动引起。 大尺度量要通过数值求解运动微分方程计算出大尺度量要通过数值求解运动微分方程计算出来来 根据经典的湍流统计理论，根据经典的湍流统计理论，认为小尺度运动（小认为小尺度运动（小涡）是完全随机的涡）是完全随机的，与平均运动和流场的边界形，与平均运动和流场的边界形状几乎没有关系，因此状几乎没有关系，因此对小尺度运动用建立模

11、型对小尺度运动用建立模型进行模拟进行模拟。 小尺度对大尺度运动的影响在运动方程中表现为小尺度对大尺度运动的影响在运动方程中表现为类似于雷诺应力项，称为类似于雷诺应力项，称为亚格子雷诺应力亚格子雷诺应力. 湍流的直接数值模拟湍流的直接数值模拟就是数值求解完整的、就是数值求解完整的、三维非定常的三维非定常的Navier-Stokes方程组，计算方程组，计算包括脉动运动在内的湍流所有瞬时流动量包括脉动运动在内的湍流所有瞬时流动量在三维流场中的时间演化在三维流场中的时间演化 最早的直接数值模拟工作开始于上个世纪最早的直接数值模拟工作开始于上个世纪70年代。年代。 到目前为止，在现有的计算能力限制下，到

12、目前为止，在现有的计算能力限制下，直接数值模拟还只能计算中等雷诺数以下直接数值模拟还只能计算中等雷诺数以下且具有简单几何边界条件的湍流流动，网且具有简单几何边界条件的湍流流动，网格可达到格可达到512512512。在这种情况下，。在这种情况下，数值模拟完全可以代替物理实验数值模拟完全可以代替物理实验。 对于湍流的直接数值模拟，对于湍流的直接数值模拟，数值方法主要数值方法主要有有有限元法、有限体积法、高阶有限差分有限元法、有限体积法、高阶有限差分法和谱方法等法和谱方法等。 对于直槽道湍流的直接模拟，由于其边界对于直槽道湍流的直接模拟，由于其边界简单，通常采用谱方法简单，通常采用谱方法 谱方法具有

13、以下优点谱方法具有以下优点：1）精度高）精度高2）具）具有准确的空间微分有准确的空间微分 谱方法的缺点是谱方法的缺点是只适用简单的几何边界，只适用简单的几何边界，不适合的边界条件不仅不能达到预期目的，不适合的边界条件不仅不能达到预期目的，还会导致计算失败还会导致计算失败Navier-Stokes方程及其离散化方程及其离散化21Reuuuput 0u( , , )(2 ,2 , ),( , , )(2 ,2 , )u x y tu xytv x y tv xyt 二维N-S方程的解22/Re( , , )cos()sin()k tu x y tkxky e 22/Re( , , )sin()co

14、s()k tv x y tkxky e24/Re( , , )0.25cos(2 ) cos(2 )k tpx ytkxky e 采用采用Adams-Bashforth-Crank-Nicolson(ABCN)方法得到离散化的方法得到离散化的N-S方方程程2111/211RennnnuupSt 21/20nnpS1123()()12RennnnnnnuuuuuSut 伪谱方法解离散伪谱方法解离散N-S方程方程 (1)在求解区在求解区 分解成个分解成个 网格点，网格点，格点间隔为格点间隔为 ，格点坐标，格点坐标为为(2) 确定的初始分布确定的初始分布. 初始的速度场初始的速度场必须满足连续性方程

15、必须满足连续性方程(不可不可压缩条件压缩条件)和一定的能谱和一定的能谱.0,2 0,2 xyNN2/,2/xyxNyN 2/,2/jxjyxj Nyj N 引入流函数引入流函数. 速度分量用流函数表示为速度分量用流函数表示为 在波数空间流函数的实部和虚部分别为在波数空间流函数的实部和虚部分别为 ,uvyx 121212121212( , )( , )cos(2 ( , )( , )( , )sin(2 ( , )RIk kak kbk kk kak kbk k 将能量平均地分配给每一个模式将能量平均地分配给每一个模式,则则 在波数空间在波数空间212( ,)( )/()ka k kE kk N

16、2212kkk121212( ,)( ,)( ,)RIk kk kik k121211()1200( , )( ,)yxNNi k x k ykkx yk k e (3)将函数将函数u(x,y,t),v(x,y,t),p(x,y,t)等场量展开等场量展开成成Fourier级数级数. 离散点离散点(j,k)(在物理空间点在物理空间点)上函数上函数u(x,y,t)的的值分别为值分别为 /2 1/2 1(),/2/2/2 1/2 12(/),/2/2yxjkxyyxxyxyNNi lxmyj kl mlNmNNNi jl Nkm Nl mlNmNuueue0,1,.,1;0,1,.,1xyjNkN

17、函数函数u(x,y,t),在谱空间的分量为在谱空间的分量为11(),00112( /),0011yxjkyxxyNNi lxmyl mj kjkxyNNi lj Nmk Nj kjkxyuueN NueN N /2,/2 1,.,/2 1/2,/2 1,.,/2 1xxxyyylNNNmNNN （4）求场量的空间导数求场量的空间导数 /2 1/2 12(/),/2/2,( )yxxyxyNNni jl Nkm Nnl mnlNmNj kd uiluedx/2 1/2 12(/),/2/2,()yxxyxyNNni jl Nkm Nnl mnlNmNj kd uim uedy(4)混淆误差修正混

18、淆误差修正. 对于对于线性微分方程，谱方法的精度取决于线性微分方程，谱方法的精度取决于谱展开的精度，即谱截断误差谱展开的精度，即谱截断误差。对于。对于非线非线性微分方程性微分方程，有限项谱展开的，有限项谱展开的非线性项会非线性项会产生附加的误差产生附加的误差，这种误差在谱方法中称，这种误差在谱方法中称为为混淆误差混淆误差（aliasing error） 假设，在离散点上的函数值分别等于 ( )( )( )u xw xu xxjx,()jjuux/2 1/2jNilxjllNuu e/2 1/2()( )jNilxjllNuil u ex0,1,.,1jN 他们的乘积他们的乘积 是在物理是在物理

19、空间中离散点上的值，这时就有高波数的空间中离散点上的值，这时就有高波数的成分混在乘积中。令在谱空间的分量成分混在乘积中。令在谱空间的分量为为()jjjuwux( )w xjxjw101jNilxljjww eN1/2 1/2 10/2/2/2 1/2 11()/2/201()( )1( ) ()jjjjNNNipxiqxilxlpqjpNqNNNNi p q l xpqpNqNjwueiqueeNu iqueN ()() ()lpqpqp q lp q l Nlpqp q l Nwiq u uiq u uwiq u u 常用的一种常用的一种消除混淆误差的方法消除混淆误差的方法成为成为“3/2规

20、则规则”，具体做法是将函数在，具体做法是将函数在Fourier空间空间中的展开的系数进行延拓到原来的中的展开的系数进行延拓到原来的3/2倍，倍，延拓后的函数的延拓后的函数的 Fourier系数有系数有M=3N/2项，项，延拓按以下的规则延拓按以下的规则-赋值赋值 llUu/2/2 1NlN 0lU /2/2MlN /2/2 1NlM 以延拓的谱系数求出延拓后函数在物理空以延拓的谱系数求出延拓后函数在物理空间中的离散值间中的离散值/2 1/2jMilxjllMUU e/2 1/2()( )jMilxjllMUil U ex0,1,.,1jM 然后在物理空间求乘积然后在物理空间求乘积在谱空间得到类

21、似的乘积公式在谱空间得到类似的乘积公式 可以证明当时可以证明当时()jjjUWUx ( )( )lpqpqp q lp q l MWiq U Uiq U U /2/2 1NlN ()()lpqpqp q lp q lWiq U Uiq u u 当求得后当求得后,将范围内的值将范围内的值赋值给赋值给,从而求得在谱从而求得在谱空间的值空间的值,再经过逆再经过逆Fourier变换后变换后,就可以求就可以求出出 lW/2/2 1NlN lw( )( )( )u xw xu xx( )( )( )u xw xu xx（）先在谱空间解方程先在谱空间解方程求出压强，然后类似解方程求出压强，然后类似解方程求出

22、速度求出速度21/20nnpS2111/211RennnnuupSt （四）流动湍流控制的数值模拟结果四）流动湍流控制的数值模拟结果2003年年2004年年管曙光管曙光先后在先后在chaos和和PRE上发表两篇流动湍流控制的文章上发表两篇流动湍流控制的文章,他采他采用的方法是用的方法是全局钉扎和局域钉扎负反馈控全局钉扎和局域钉扎负反馈控制方法制方法,我们的工作就是在此基础上完成我们的工作就是在此基础上完成.下下面介绍我们的两个工作面介绍我们的两个工作 (1)用移动控制器控制流动湍流用移动控制器控制流动湍流 N-S方程方程 控制方案控制方案 两种目标态：两种目标态：层流态层流态 空间周期态空间周

23、期态 21ReuuupuFt 0u()TFuu5110.05,0TTuv2 /Re2( , , )0.1cos( )sin( )tTux y txy e 2 /Re2( , , )0.1sin( )cos( )tTvx y txy e 控制器放置：控制器放置：有控制信号注入的网格点称有控制信号注入的网格点称为控制器为控制器，在这个工作中，在这个工作中，控制器是均匀控制器是均匀分布在计算区内分布在计算区内，并且对两个速度分量都，并且对两个速度分量都加控制，总控制器个数为加控制，总控制器个数为 初始能量分布初始能量分布初始能量初始能量 系统自然演化系统自然演化5个时间单位后才对系统进行个时间单位后

24、才对系统进行控制控制 2(256)xyMI I20( /)( ,0)k kE kcke05.0k 0.01E 参考态（自由演化） 静态控制器控制湍流，目标态是周期态(a)(b) (a) 是 与t的关系曲线,(b)t=50是时刻的涡度等值线图.(c)(d) (c) 是 与t的关系曲线,(d)t=80是时刻的涡度等值线图8,4xyIIln( )t8,8xyIIln( )t 用运动控制器控制湍流用运动控制器控制湍流选控制器为参照选控制器为参照系系 ，N-S方程为方程为 控制器运动速度控制器运动速度211( )()ReMiTiuuupVuurr tuut (1.2,0)xyVVV移动控制器控制湍流，目

25、标态是周期态(a)(b) (a) 是 与t的关系曲线,(b)t=50是时刻的涡度等值线图.(c)(d) (c) 是 与t的关系曲线,(d)t=80是时刻的涡度等值线图8,4xyIIln( )t8,8xyIIln( )t 控制效率的分析控制效率的分析 涡度差涡度差控制差控制差控制差随时间变化控制差随时间变化暂态时间暂态时间控制开始时刻控制开始时刻是的时刻是的时刻 ( , , )( , , )Tx y tx y t1/2256221,11( )( , , )256ijti j t( )( )/(0)tt21tt15t 2t35 10 tet )(运动控制器控制效果描述目标态是层流态.(a)是随 l

26、n(M) 变化的关系曲线 . (b)是随ln(M)变化的关系曲线 ln(100) 静态控制与移动控制暂态比较(a)目标态为层流态情况下的与ln(M)关系曲线(b)目标态为周期态情况下的与ln(M)关系曲线(c) 与移动速度关系曲线( ) 8,8xyII32 32M xV（）用间歇反馈控制流动湍流（）用间歇反馈控制流动湍流 N-S方程方程 控制方案控制方案(全局钉扎全局钉扎) 目标态：空间周期态目标态：空间周期态222222221()Re1()Reuuupuuuvftxyxxyvvvpvvuvtxyyxy ()Tfuu2 /Re( , , )cos( )sin( )tTux y txy e 2

27、/Re( , , )sin( )cos( )tTvx y txy e1 初始能量分布与前一种方法相同初始能量分布与前一种方法相同,速度场的速度场的初始能量为初始能量为 控制差控制差 00.1E 1/225621/2221,11( )( , , )( , , )256ijFF x y tF x y t ( , , )( , , )( , , )TF x y tF x y tFx y t, ,Fu vvuxy 无间歇全局钉扎反馈控制的控制差随时间t变化曲线 无间歇反馈控制存在无间歇反馈控制存在 间歇反馈控制方案间歇反馈控制方案式中式中 是间歇周期是间歇周期. 是整数是整数, 是时间步长是时间步长,定义定义 工作时间比率为工作时间比率为50limFFtc1(),0Tfuutkkktfother012()ktkkt 120,k k k k0.0025t 10/kk 间歇反馈控制. 控制参数 (a)控制差随时间变化.(b)涡度等值线图0.5,0.55间歇反馈控制湍流.(a)固定控制参数 ,t=50时刻的涡度控制差与工作比率关系曲线.(b)(c)(d)t=50涡度差 等值线图,控制参数 (b) (c) (d)T0.50.350.551.0 速度分