cfd模拟中网格尺度的误差分析

cfd模拟中网格尺度的误差分析

0工程模拟问题的数值模拟研究随着cfd数值模拟软件成为学术研究、分析和设计中不可或缺的工具，cfd模拟结果的可靠性也越来越受到重视。用CFD软件对一个工程问题进行数值模拟时,模拟的结果与真实的结果之间总会存在各种误差。对于所有这些误差,若能合适地估算出其值大小,就可对计算结果进行相应的修正,获得更为精确的计算结果。对此,文献进行了比较全面的研究。本文首先回顾了CFD模拟中可能会产生的各种误差和误差估计,然后采用Richardson外推法建立了估算网格误差的方法和准则,并以NACA4412翼型为例介绍了如何应用误差估算方法及其对CFD模拟结果的修正,更多应用见文献。1偏微分方程及其初边值问题的建模误差任何流体力学问题都要受一些物理学定律所支配,比如质量守恒定律、动量守恒定律和能量守恒定律等。对于非定常流动,这些物理学定律形成了初边值问题(IBVP);而对于定常流动则形成了边值问题(BVP)。一个初边值问题可由基于这些物理学定律所建立的一组连续偏微分方程及初边值条件所定义。若采用算子L表示这个方程组,一个初边值问题可写为:其中第一个式子是控制方程,第二、三个式子是初边值条件。T是真实值。x和t分别是几何坐标和时间。按照定义,该式中没有包含任何误差。求解流体力学问题有三种方法:实验法,即实验流体动力学(EFD);解析法,即理论流体动力学(TFD);模拟法,即计算流体动力学(CFD)。任何一种方法的应用都会伴随相应的误差产生。比如:实验流体动力学是通过实验测量系统得到测量值D,在这个过程中会产生偏差和精度误差,导致实验测量值D中的实验误差δD(δD=D-T)产生。而分析和模拟方法,或理论流体动力学和计算流体动力学,首先要选择所研究物理问题的数学模型,如N-S方程、欧拉方程、流函数方程等等,通过偏微分方程及其初边值条件建立该物理现象的模型,形成与式(1)对应的初边值问题:所有的定义与式(1)类似,但算子LM、GM、HM和xB都可能包含了建模假设。所涉及的假设可能会有:有/无黏性、定常/非定常、压缩性、湍流、非牛顿流体、化学反应、多相流和稀薄气体等等。所以M≠T。从而将产生初边值问题的建模误差:δSM=M-T。建模误差是在建立模型的过程中,或者由于人们对物理现象的认识不准确或不全面,或者由于计算或分析的条件所限制而必须做出一些简化,从而使所建立的模型未能准确和全面地描述实际物理现象而产生,因此,这部分误差难以估计。本文将不对这部分误差和不确定性进行深入研究,而假设所要面对的流体力学问题均很好地服从于NavierStokes方程。对于湍流问题,所采用的湍流模型也完美地描述了真实的湍流流动(虽然这个假设明显与很多实际问题不符),而着重研究在数值模拟的过程中所产生的误差。理论流体力学是在式(2)基础上,通过数学解析的方式获得其解析值。虽然解析的过程可能是精确的,但为了能获得解析解,不仅要对式(2)进一步简化,而且有可能要对几何模型和边界条件进行大量简化,从而必然会产生与真实值的差异,导致很大的建模误差。模拟的方法是通过采用数值技术,把连续的初边值问题转化为离散的初边值问题,即通过使用空间和时间离散技术(比如有限差分、有限体积和有限元等),近似地求解式(2)。从而导致由于空间间隔Δx和时间间隔Δt等产生的数值误差(即:当这些间隔趋于零时,数值误差是零)。离散的初边值问题可以通过一个离散算子LN和离散的初边值来定义为:其中源项ΓN与把连续方程变成离散方程过程中所产生的截断误差以及数值求解过程中的残差相关。通过计算机的算法语言编制计算机的执行程序,即CFD计算程序,求解式(3),就得到了模拟预估值S。其结果可能会包含多种误差,这些误差称为模拟误差,即模拟值与真实值之差,δs=S-T。对于数值模拟误差δs,在很多情况下可以分解为如下两部分:其中第一部分是对δs的符号和量值的估计,而εS是在该估计过程中所产生的误差。按照式(4)估计出数值模拟误差以后,就可对数值模拟结果S进行误差修正。如果认为修正后的模拟值为模拟基准值并用Sc表示(在最理想条件下,Sc就是精确解),即Sc为数值模拟结果减去可估计的数值误差。在对由偏微分方程组构成的数学模型进行数值求解时,可能会产生的数值误差包括:网格尺度误差(δG)、时间步长误差(δT)、迭代误差(δI)、以及一些输入参数所引起的误差(δP),等等。于是,如果同时考虑式(4)和式(5),修正后的模拟值Sc可表示为:可见,修正后的模拟值中包含了真实值T和不可估计部分εS(注:这里没有考虑建模误差)。式(6)适合于计及网格尺度误差、时间步长误差、迭代误差和输入参数误差等数值模拟结果的修正。本文只考虑由于网格尺度所产生的误差和不确定性。这样为了标记简单,就可在式(6)中略去下标,而把由网格尺度引起的误差用δ*表示,于是,对于任意一个网格,其模拟值S与模拟修正值SC间的关系可表示为2古龙的外推法和应用标准2.1网格大小对误差估计的影响Richardson假设一个方程的离散解等于其精确解与离散间隔(步长)的幂级数之合,即其中,gi与方程离散的方式无关。显然,对于式(7)中的误差部分,可以采用Richardson外推法估计。采用类似的方法,式(7)中由网格尺度引起的误差δ*也可以用级数表述,然后可得到广义的Richardson外推法。若只保留级数前n项,其中Δx表示网格的间隔,p(i)表示精度等级,g(i)是该级数的系数。由式(9)可知,对于某个计算,S是其计算结果,为已知,Δx也是已知的。对于求和号内的每一项有两个未知数p(i)和g(i),再加上Sc,式(9)总共有m=2n+1个未知数,需要建立m个方程才可求解,即需要进行m个不同尺度网格的计算。本文仅对n=1进行分析,忽略高阶项的影响,此时需要对(9)的误差项δ*进行修正,修正后的误差项可以表示为,详细内容可参阅文献对于n=1,式(9)可简化为对于m=3个计算,并忽略高阶项,有,其中,r为不同网格尺度间的比值,S1和S3分别为最细和最粗网格的计算结果。在这三次模拟计算中,已经假设了其它参数对每次模拟结果的影响是相同的,即在这个过程中三次所产生的误差只是由网格尺度引起的。联立求解可得:δ*所表示的是通过Richardson外推法估计的由于网格的变化在n=1的精度下所产生的误差。由于它忽略高阶项,为了提高误差估计的精度,需要把高阶项的影响也考虑进去,需引进一个对式(11)进行修正的因子:其中的pest是对高阶项指数的估计。所以修正以后,由于网格的变化所产生的误差可以表示为:其中ε21=S2-S1;ε32=S3-S2。由式(13)和式(14)就可进行网格尺度引起的误差的估算,或对CFD结果的修正。2.2多重网格技术为了进行上述误差的估算,需要进行若干次不同尺度(或密度)网格的计算。为了便于分析,最好使不同网格尺度间的比值r保持一致,即:为了要确保最粗网格的解是定性正确的,r值宜不能过大。本文采用多重网格技术形成不同尺度的网格,即只制作最细网格,通过依次取消三个网格结点中间的结点而获得其它各个层次的网格。如图1所示。总体上,可近似认为全计算域内r近似为2。2.3收敛性条件设定对工程问题的CFD计算,估计由网格尺度所带来的误差,需要对若干个不同网格数目的网格计算结果的收敛性进行分析。为了进行收敛性的分析,至少要进行三个不同网格的计算,获得S1,S2和S3三个解。它们之间的差及其比率R为于是收敛性条件为:0<R<1(单调收敛);R<0(振荡收敛);1<R(发散)。对于单调收敛,可采用广义Richardson外推法估算误差;如果是振荡收敛,只能由振幅获得误差的范围,即不确定度;对于发散情况误差和不确定度均不能进行估算。2.4最粗网格的计算在采用上面的方法时,为了能得到高可信度的误差估计,需要注意以下几点应用准则:(1)在进行外推求解时,要确保各网格的模拟值单调收敛,即0<R<1。(2)确保每个网格的迭代计算很好收敛,否则会带来较大的数值计算残差。(3)最粗网格的CFD结果要与精确解尽可能地接近,或至少要确保最粗网格的结果定性上正确。(4)如果采用低雷诺数湍流模型,要使得最粗网格能保证离壁面第一层网格处的y+值满足模型的要求。(5)为了使各个网格的计算均具有高的收敛速度,在最粗网格计算时也能够采用多重网格加速技术,在制作最细网格时要设置足够的网格数目和网格重数。(6)原则上讲,这个方法不仅可用于总体特性参数(升力、阻力、压比、效率、功率等)误差的估计,也可有用任意空间点流动参数(速度、压力和温度等)计算误差的估计,只不过对于后者,为了得到一个空间坐标点的数据,需要对各个网格密度的计算结果进行插值,从而可能会带来新的误差。3fd模型及计算对NACA4412翼型做了多个网格尺度的计算及误差分析。CFD计算采用NUMECA的Fine/Turbo软件包,其求解时间相关的雷诺平均N-S方程,采用二阶中心格式离散。选用低雷诺数SA湍流模型。采用多重网格加速。低速流动采用预处理方法。3.1细网格派出不同阶次的粗网格总共进行了6个不同网格数的计算。为了保证各个网格分布规律一致,采用多重网格技术由细网格派生出不同阶次的粗网格,各个计算的网格数目及其所对应壁面的y+如表1所列。对于该二维计算,网格指数由两个数字分别表示在I和J方向的指数。其中“0”表示对应方向的最细网格,“1”为次细网格,依次类推。详细的计算分析参见文献。3.2网格冲角条件下不同网格组合升力系数的收敛性分析图2和图3分别给出了不同网格计算的翼型升力系数和-6°冲角的翼型压力系数与实验数据的比较。由表1可见,前三个网格,满足低雷诺数湍流模型的要求,但后两个网格的y+过大,不满足准则3;从图3的压力分布可见,网格G44的结果明显不正确或定性上不正确,不满足准则2。总之,网格G00,G11,G22和G33的计算结果可用于计算误差的估计。下面分别对冲角为-6°,0°,6°,14°和18°的结果进行分析。为了展示网格粗密程度对误差估计精度的影响,将分别对G00、G11和G22以及G11、G22和G33两个组合进行分析对比。所分析的数据分别是模拟计算的Cl和Cd、实验值D和修正值Sc。表2列出了各个冲角条件下不同网格组合计算的升阻力系数的收敛性判断。可见,对小于最大升力冲角的各个计算,细网格组合(G00-G22)均为单调收敛;对于小冲角情况(0°和6°),粗网格组合也为单调收敛,可以采用Richardson外推法估计升力系数的计算误差。但对于14°和18°冲角的计算,两个网格组合均为发散或振荡收敛,不能进行升力系数的计算误差估计。对于各个冲角的阻力系数的计算,两个网格组合均为单调收敛,均可进行误差估计。图4和图5分别给出了12°冲角条件下升阻力系数计算值随网格密度的变化曲线及其与计算修正值和实验数据的比较。可以看到,采用细网格组要比粗网格组所得到的误差修正模拟值更精