低阶格式在大涡模拟计算中的适用性

刘同新， 马宝峰

（北京航空航天大学流体力学教育部重点实验室，北京 100191）

低阶格式在大涡模拟计算中的适用性

刘同新， 马宝峰

（北京航空航天大学流体力学教育部重点实验室，北京 100191）

采用三维Taylor-Green涡作为研究对象，利用工程中常用的低阶数值格式，研究格式本身的数值误差对大涡模拟计算的影响.结果表明：三种数值格式的数值耗散行为都与亚格子模型行为类似，即在小雷诺数下，流场比较光滑时，耗散很小，当雷诺数增加，流动转捩为湍流，流场梯度增大，耗散显著增大.对于MUSCL格式和二阶有界中心格式，在高雷诺数下，亚格子尺度模型没有明显改善计算结果，但也没有使计算结果恶化.中心格式相比其它两种格式，数值耗散最小，但是在高雷诺数湍流情况下，中心格式的数值耗散仍然主导了能量的耗散，再添加亚格子模型，计算结果反而变得稍差.对于工程中的低阶格式而言，采用中心格式计算大涡模拟是比较好的选择，而且在计算不存在稳定性问题时，采用不添加亚格子模型的隐式大涡模拟效果更好.

大涡模拟；数值耗散；亚格子模型；低阶格式；

0 引言

湍流的准确计算在航空、航天器设计和研制中具有重要意义，因而能够准确地模拟湍流一直是学术界和工业界追求的目标.在工程湍流计算中以往大都采用基于雷诺平均的湍流模式（RANS），由于湍流模式是对湍流内大小不同尺度的流动结构统一建模，难以得到普适模型，计算准确度欠佳.同时雷诺平均基于时间平均，抹掉了非定常脉动细节，无法计算转捩、气动噪声等非定常问题.直接数值模拟方法（DNS）尽管计算准确，并可获得流场的全部信息，但计算资源的要求过高，在可以预见的将来难以应用到工程计算中，目前主要用来对低雷诺数、简单几何外形的流动做机理研究.

大涡模拟（LES）是目前最有希望取代雷诺平均来进行湍流工程计算的方法［1］.其在相对较低的计算资源下能够模拟较高的雷诺数运动，同时能够获得非定常流场信息.大涡模拟不对时间求平均，而是在空间上对湍流小尺度涡进行滤波，因而流动在时间方向上大涡脉动信息不丢失.滤波掉小尺度涡后由于切断了湍流从大涡向小涡的级串，因而能量会在截断的尺度上堆积，产生非物理解或数值不稳定，因而在计算时需要设计一种机制，将累积的能量耗散掉.用来耗散小涡能量的模型称为亚格子尺度模型（subgrid-scale model，SGS）.亚格子模型是针对小涡建模，而小涡运动更容易趋向均匀各向同性，因而更容易获得相对普适的模型.

对于亚格子尺度模型，除了本身的建模误差外，模型提供的物理耗散与计算格式误差之间干扰是一个需要解决的重要问题.Ghosal［2］通过理论分析指出即使是四阶精度的中心格式，其截断误差与亚格子应力也是同一量级的.所以为了准确模拟亚格子应力，大涡模拟一般要求采用高精度低耗散的数值格式进行计算或者采用显示滤波方法［3］滤除网格分辨率内的高波数段模拟不准的流场.高精度格式和显示滤波方法虽然在学术研究领域，针对简单几何外型的流场取得较好的计算结果，但出于计算效率和鲁棒性的原因，一直没有广泛应用到工程计算中.工程领域采用的仍然是二阶到三阶精度的低阶格式.基于数值格式误差对亚格子模型的干扰较大，有研究人员提出不显式建立亚格子尺度模型，而利用数值格式自身的数值耗散作为隐式模型模拟SGS应力，该大涡模拟方法称为隐式大涡模拟.隐式大涡模拟类似于超声速计算中的激波捕捉方法，即都是用数值格式自身的耗散来模拟物理问题.

Boris等［4］首次分析了单调格式的截断误差，发现其耗散行为与SGS应力很类似，因而提出可利用格式的截断误差隐式作为SGS模型的想法，而不再建立显式SGS模型.但正如Grinstein等［5］所说的，并不是所有的数值格式不加显式SGS模型，在粗网格上计算都叫做隐式大涡模拟，只有某些特定的数值格式才能够得到较好的计算结果.多种数值格式已被验证适合做隐式大涡模拟，可参照相关综述文章［5］.Garnier等［6］采用均匀各向同性湍流，研究了一系列迎风格式在大涡模拟中的适用性，发现迎风类格式的数值耗散行为与SGS模型类似但不完全一样，即使是5阶WENO格式的数值耗散已足够大，不建议另外再加显式SGS模型.而Hickel［7］等进一步指出，利用已有的格式采用逐个尝试的方法研究隐式大涡模拟有些盲目，他们采用改型方程分析方法，针对格式的截断误差性质，专门建立了一种适合做隐式大涡模拟的新格式，针对模型流动，取得了很好的模拟结果.Kravchenko and Moin［8］对渠道流动进行大涡模拟误差分析，他们利用谱方法进行分析，发现在低阶格式下误差对于亚格子模型具有较大的影响，但仍然能够得到较好的计算结果.杨小龙和符松［9］采用pade格式，研究了低阶到高阶格式数值误差对大涡模拟的影响，结果表明亚格子模型能够缓解数值误差对计算结果的影响，采用低精度格式也可能够得到较好计算结果.谢志刚等［10］也采用二阶格式得到了较为满意的计算结果.

综上，目前对于数值耗散对大涡模拟计算结果的影响，结论尚不一致，不同研究者采用不同的计算格式和方法可能会得到不同的结论，因而需要开展更广泛的研究，对不同的数值格式特别是工程中常用的低阶格式在大涡模拟中的适用性，采用不同的方法进一步进行验证，从而为工程计算提供参考.本文采用一种后验的方法（a posteriori），以Taylor-Green涡为计算对象，对工程中常用的低阶数值格式的数值耗散对大涡模拟的影响进行评估.

1 数值方法

1.1 主控方程及亚格子尺度模型

本研究的所有计算基于FluentTM14.5代码，该代码目前在工业界应用较多.

计算基于不可压N-S方程，大涡模拟需要对方程进行滤波处理，滤波后的方程为

大涡模拟的亚格子尺度建模既是针对ij进行，本研究的亚格子尺度模型选用动态Smagorinsky模型，这是目前大涡模拟计算中的主流模型，在计算效率和准确性方面较适合于工程计算.动态Smagorinsky模型是将经典Smagorinsky模型中的系数Cs在计算过程中根据大涡信息自动计算［11－12］.

经典Smagorinsky模型基于Boussinesq涡粘假设

动态亚格子建模过程是一种方法，而不是单一的模型，该过程可应用于任何常系数亚格子模型，将其变为动态模型.应用于经典Smagorinsky模型则构成动态Smagorinsky模型.该方法的基本思路是对流动方程进行二次滤波，根据第二次滤波及大涡模拟第一次滤波信息自动计算涡粘系数，从而保证亚格子模型可适应更广泛的流动计算.

1.2 离散格式

数值离散基于非结构网格有限体积法，不可压流动方程采用分数步方法（fractional step）［13］求解.相比于SIMPLE类方法，分数步方法不将算子分裂误差通过迭代降到接近零，而是将分裂误差降到与时间截断误差为同一量级.这样能够大大减少在每个时间步内推进的迭代次数，从而在非定常计算中能够提高计算效率.该算法时间截断误差和算子分裂误差都由时间步长确定，因而时间步长不宜太大.本研究时间离散采用隐式二阶格式，具有二阶精度，时间步长取为0.001.

本计算空间对流项采用三种离散格式，包括二阶中心格式、三阶MUSCL格式和二阶有界中心格式.粘性项离散采用二阶中心格式.

Fluent程序中的有限体积法基于格心法，即计算变量位于各网格单元的几何中心.不同的空间格式体现为，由相邻单元中心变量值，计算界面上的值所用的差值或重构方法.

二阶中心格式：二阶中心格式具有较少的数值耗散常被用在大涡模拟中.该格式利用界面两侧单元中心的值计算界面变量值.对于界面流动变量，基于二阶中心格式的界面重构公式为

其中0和1代表了相邻两个单元，其中中间的公用面为f，∇φr，0和∇φr，1是位于单元0和1中心的变量的梯度，其中r为单元中心到面中心的向量.单元中心的变量计算采用Taylor展开法［14］.

因为二阶中心格式奇偶失联，计算时会出现非物理振荡，因而需要做校正（deferred correction），φf，up为低阶迎风格式计算的结果，“old”为上一时间步的结果，当计算收敛后，结果收敛到二阶中心格式.

三阶MUSCL格式：MUSCL格式是通过中心格式和二阶迎风格式混合构造的格式，最早由Van Leer提出［15］.MUSCL格式相比二阶迎风格式，具有更小的数值耗散，对网格类型适应性更强.

其中φf，up＝φ＋∇φ·r，为二阶迎风格式.

二阶有界中心格式：由于纯中心格式数值耗散较小，即使采用校正，在梯度大的地方计算结果中仍会出现非物理振荡.而有界中心格式的数值耗散大于纯中心格式，但比中心格式鲁棒性好，一般不会出现非物理振荡.有界中心差分格式是由变量归一化方法（NVD）［16］并在对流有界准则（CBC）［17］约束下构造的复合格式，该复合格式由中心格式、中心和二阶迎风混合格式、一阶迎风格式复合而成.

1.3 Taylor-Green涡

Taylor-Green涡（TGV）是研究湍流转捩、级串和耗散的一种重要的模型流动［18］.被广泛用来验证大涡模拟的计算正确性.相比均匀各向同性湍流，TG涡包含的流动机制更为丰富，计算难度也更大.

TG涡初始为一光滑的自由周期涡，可视为层流，开始演化后，在非线性作用下，光滑涡结构失稳，转捩为湍流，并不断级串为更小的涡，同时随着时间的演化，能量不断耗散.

本研究选用三维Taylor-Green涡（TGV）进行计算，分析湍流发展过程中的能量耗散.其初始流场为

初始涡量场可通过ω＝∇×u得到.

计算域为边长为2π的正方体，网格点为64×64×64，正方体对应的两个面形成周期性边界条件.根据TG涡直接数值模拟的结果，该网格尺度可保证在所计算的雷诺数下，网格对涡系的截断波数处于惯性区，从而能够保证大涡模拟的网格尺度是足够的.取特征速度为1，特征尺度为1，定义雷诺数为1／ν.计算过程中，雷诺数通过分子粘性系数来改变.

TG涡系统是一个自耗散系统，在演化过程中能量不断耗散.在数值模拟的情况下，能量耗散的来源分为流体分子粘性耗散，亚格子模型耗散和数值耗散.

分子粘性耗散率 εmol＝2ν＜SijSij＞，＜·＞代表体积平均；

亚格子模型耗散率 εSGS＝2νSGS＜SijSij＞，

其中ν，νSGS，Sij分别分子粘性系数，亚格子湍流涡粘系数和应变率张量，νSGS＝μSGS／ρ．分子耗散率和亚格子耗散率一起称为物理耗散率εphy＝εmol＋εSGS．

2 计算结果及分析

图1显示TG涡在不同时刻的演化结构，初始为规则的周期多涡结构，随时间演化，逐渐变为小尺度旋涡.图中涡面是利用涡的λ2判别法［19］得到，其灰度代表了动能的大小.

图2给出了不同时间步长对计算结果的影响，可见所取时间步长已足够小，改变时间步长后总能量随时间的衰减曲线完全重合，本文的以下计算时间步长都取为0.001.图3给出了不同雷诺数下总平均动能随时间的演化.在初始动能一样的情况下，雷诺数越大，耗散越慢，湍流动能在惯性作用下会传递到更小尺度的涡上去.

通过与DNS结果的对比，可以比较各种差分格式以及亚格子模型的计算效果，本研究选用BRACHETT［20］的TG涡直接数值模拟结果作为对比数据.BRACHETT［20］的关于TG涡的DNS数据是该方向的经典数据，由谱方法计算，最初计算时采用2563网格，约十年后又基于8643网格做了校核［21］，但未发现偏差.长期以来一直作为TG涡模拟的对比数据，即使最新DNS计算也与其复合得很好［22］.

图4给出了Re＝100情况下，三种数值格式计算的TG涡系统能量耗散率随时间的演化.给出了物理耗散率和总耗散率在有无亚格子模型情况下的变化规律，并给出的DNS数据作为对比.可见5条曲线基本上重合到一起，计算结果与DNS数据符合很好，无论加不加亚格子模型，流动的物理耗散率都等于总耗散率，即数值耗散很小.实际上由于在该雷诺数下，流动演化可看成层流演化，耗散过程主要是分子粘性耗散，亚格子耗散可以忽略不计.

图5给出了Re＝400情况下的计算结果.对于每一种计算格式的计算结果而言，此时物理耗散率和总耗散率都出现差异，特别是在耗散率峰值附近（t＝6 s），差异显著，表明数值耗散此时对计算结果有明显的影响.对比三种数值格式，MUSCL格式数值耗散的影响最大，t＝4 s以后，能量总耗散率与物理耗散曲线之间明显分散（图5（a）），但是加不加SGS模型对总耗散率和物理耗散率影响不大.中心差分格式和有界中心格式数值耗散的影响相对较小.考虑施加亚格子尺度模型和不加模型两种情况，采用中心格式计算的物理耗散值与总耗散率的差值是最小的，表明中心格式的数值耗散相对较小，与DNS的结果也比较接近.

当Re＝3 000时，数值耗散的影响显著增加，如图6所示.对于MUSCL格式和有界中心格式而言，施加亚格子模型与否，总能量耗散率相差不大.对比DNS数据，亚格子模型仅在局部时间段对耗散率有所改善，但并不明显.而对于中心格式，不加SGS模型时，总能量耗散率与DNS结果反而更接近，尽管此时耗散仍然偏大，耗散率峰值出现的时间提前，但是峰值的大小与DNS结果很接近.上述结果表明对于三种数值格式而言，在该雷诺数下，数值耗散基本主控了耗散率的量值，且大于亚格子尺度模型提供的耗散.基于这一点可以看出，对于所采用三种低阶格式而言，利用格式本身的数值耗散来作为亚格子模型，构成隐式大涡模拟是更为实用的技术，而再显式添加亚格子模型的效果并不大，甚至会使计算结果变差（比如中心格式）.

对于三种数值格式而言，中心格式的总耗散率与DNS结果最为接近，特别是不加SGS模型时.更为重要的，采用中心格式时，在添加亚格子模型的情况下，中心格式的物理耗散占总耗散率的比重最大（图6（c）），有界中心格式其次，MUSCL格式最小.这表明中心格式的数值耗散对大涡模拟计算的影响最小.不加亚格子模型时，由于计算物理耗散率时少了亚格子耗散率，所以物理耗散率量值很小.

需要注意，无论加不加亚格子模型，以及采用三种之中的任一数值格式，计算得到的耗散率相比DNS结果在大部分区域都偏大，而在峰值附近又偏小.说明在高雷诺数下还没有做到准确模拟湍流的耗散.根据已有的研究结果［7］，本文用到动态Smagorinsky模型本身的耗散率相比DNS结果本来就偏大，即该亚格子模型本身也存在建模误差.因而即使采用高精度低耗散数值格式，显示添加动态Smagorinsky亚格子模型，计算的耗散率与图6（c）中计算的总耗散率在趋势上也是类似的，可参见文献［7］中四阶中心格式的计算结果.

3 结论

以Taylor-Green涡的时间演化为计算对象，研究了三阶MUSCL格式、二阶有界中心格式和二阶中心格式的数值误差对大涡模拟计算结果的影响.得到如下结论：

1）三种数值格式的数值耗散行为都与亚格子模型行为类似，即在小雷诺数下，流场比较光滑时，耗散很小，当雷诺数增加，流动转捩为湍流，流场梯度增大，耗散显著增大.与直接数值模拟结果相比较，对于添加和不添加亚格子模型两种情况能量耗散都偏大.

2）对于MUSCL格式和二阶有界中心格式，在高雷诺数下，亚格子尺度模型没有明显改善计算结果，但也没有使计算结果恶化.中心格式相比其它两种格式，数值耗散最小，但是在高雷诺数湍流情况下，中心格式的数值耗散仍然主导了能量的耗散，再添加亚格子模型，计算结果反而变得稍差.

3）对于工程中的低阶格式而言，采用中心格式计算大涡模拟是比较好的选择，而且在计算不存在稳定性问题的前提下，采用不添加亚格子模型的隐式大涡模拟效果更好.

［1］张兆顺，崔桂香，许春晓.湍流大涡数值模拟的理论和应用［M］.北京：清华大学出版社，2008.

［2］Ghosal S.An analysis of numerical errors in large-eddy simulations of turbulence［J］.Journal of Computational Physics，1996，125（1）：187－206.

［3］Gullbrand J，Chow F K.The effect of numerical errors and turbulence models in large-eddy simulations of channel flow，with and without explicit filtering［J］.Journal of Fluid Mechanics，2003，495：323－341.

［4］Boris J P，Grinstein F F，Oran E S，Kolbe R L.New insights into large eddy simulation［J］.Fluid Dynamics Research，1992，10：199－228.

［5］Grinstein F，Margolin L，Rider W.Implicit large eddy simulation：Computing turbulent flow dynamics［M］.Cambridge，UK：Cambridge University Press，2007.

［6］Garnier E，Mossi M，Sagaut P，Comte P，Deville M.On the use of shock capturing schemes for large-eddy simulation［J］. Journal of Computational Physics，1999，153：273－311.

［7］Hickel S，Adams N A，Domaradzki J A.An adaptive local deconvolution method for implicit LES［J］.Journal of Computational Physics，2006，213：413－436.

［8］Kravchenko A G，Moin P.On the effect of numerical errors in large eddy simulations of turbulent flows［J］.Journal of Computational Physics，1997，131（2）：310－322.

［9］杨小龙，符松.直接数值模拟／大涡模拟中数值误差影响的研究［J］.应用数学和力学，2008，29（7）.

［10］谢志刚，许春晓，崔桂香，王志石.方柱绕流大涡模拟［J］.计算物理，2007，24（2）：171－180.

［11］Germane M，Piomelli U，Moin P，Cabot W H.A dynamic subgrid-scale eddy viscosity model［J］.Physics of Fluids，1991，（4）：633－635.

［12］Lilly D K.A proposed modification of the germano subgrid-scale closure model［J］.Physics of Fluids，1992，（4）：633－635.

［13］Glaz H M，Bell J B，Colella P.An analysis of the fractional-step method［J］.Journal of Computational Physics，1993，108：51－58.

［14］Barth T J，Jespersen D.The design and application of upwind schemes on unstructured meshes［C］.AIAA 27th Aerospace Sciences Meeting，AIAA－89－0366，Reno，Nevada，1989.

［15］Van Leer B.Toward the ultimate conservative difference scheme.IV.A second order sequel to Godunov's method［J］. Journal of Computational Physics，1979，32：101－136.

［16］Leonard B P.The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection［J］.Computer Methods in Applied Mechanics and Engineering，1991，88：17－74.

［17］Gaskell P H，Lau A K C.Curvative-compensated convective transport：SMART，a new boundedness-preserving transport algorithm［J］.International Journal for Numerical Methods in Fluids，1988，8：617－641.

［18］Taylor G I，Green A E.Mechanism of the production of small eddies from large ones［C］∥Proceedings of the Royal Society of London，Series A，Mathematical and Physical Sciences，1937，158：499－521.

［19］Jeong J，Hussain F M.On the identification of a vortex［J］.Journal of Fluid Mechanics，1991，285：69－94.

［20］Brachet M E，Meiron D I，Orszag S A，Nickel B G，Morf R H，Frisch U.Small-scale structure of the Taylor-Green vortex ［J］.Journal of Fluid Mechanics，1983，130：411－452.

［21］Brachet M，Meneguzzi M，Vincent A.Politano H，Sulem P-L.Numerical evidence of smooth self-similar dynamics and possibility of subsequent collapse for three-dimensional ideal flow［J］.Physics of Fluids，1992，（4）：2845－2854.

［22］Chapelier J-B，Plata M D L L.Renac F.Inviscid and viscous simulations of the Taylor-Green vortex flow using a modal discontinuous Galerkin approach［A］.AIAA－2012－3073，2012.

Suitability of Low-order Numerical Schemes in Large Eddy Simulations

LIU Tongxin，MA Baofeng
（Ministry-of-Education Key Laboratory of Fluid Mechanics，Beihang University，Beijing 100191，China）

Several low-order numerical schemes were evaluated for their suitability in large-eddy simulations based on 3D Taylor-Green vortex，with and without a subgrid-scale model.It shows that dissipation characteristics of three numerical schemes used are similar to a subgrid-scale model.At lower Reynolds numbers，flow fields are relatively smooth，numerical dissipation is lower；As Reynolds number increasing，transition to turbulence occurs and numerical dissipation grows greatly.For MUSCL and bounded centered schemes，subgrid-scale model has a little influence on results at high Reynolds numbers.The second-order central scheme exhibits lower dissipation，but at higher Reynolds numbers numerical dissipation still dominates total energy dissipation of flow.With an explicit subgrid-scale model the results even become worse.Therefore，for large eddy simulations in engineering，the second order central scheme is suitable，particularly without an explicit subgrid-scale model.

large eddy simulation；numerical dissipation；subgrid-scale model；low-order numerical schemes

date： 2013－07－01；Revised date： 2013－10－02

O368

A

1001-246X（2014）03-0307-07

2013－07－01；

2013－10－02

国家自然科学基金（11272033）资助项目

刘同新（1987－），男，山东，硕士生，研究方向：非定常流动的大涡模拟研究，E-mail：bf-ma＠buaa.edu.cn