第10讲-有限体积法2.ppt

计算流体力学讲义第十讲有限体积法（2）李新亮lixl；力学所主楼219；82543801,知识点：,1,讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”讲课录像及讲义上传至网盘http:/cid-,CopyrightbyLiXinliang,近似Riemann求解器HLL/HLLC中心型有限体积法人工粘性具体问题的计算翼型绕流(边界条件、具体解法、粘性项处理),CopyrightbyLiXinliang,2,知识回顾,在以某节点为中心的控制体上积分,i,j,k,非结构网格的控制体,i+1,j,i-1,j,i,j+1,i,j-1,k3,k1,k2,k4,k5,结构网格的控制体,x,y,n,体积平均,控制体边界垂直于节点连线（也可选其他方式）,垂直平分线,n,1)建立控制体,2)在控制体上积分,离散方程,重构：由节点上平均值给出函数分布，最终给出通量,表示第m个界面上的值,1.有限体积法的离散过程重构,1）重构两种不同的重构方案，向左偏及向右偏。给出两种结果：及,CopyrightbyLiXinliang,3,i,j,i+1,j,i-1,j,i,j+1,i,j-1,n,左重构,右重构,2）由左右重构得到的自变量：和给出通量方案A:FVS方案B:解Riemann问题（常用）,例如：0阶重构：,线性重构:,用i,i-1点的值插i+1/2点的值（网格剧烈变化时，应当用实际坐标插值）,用i,i+1点的值插i+1/2点的值,x,y,看似二维Riemann问题，其实是一维的，坐标旋转一下就行了,2.迎风型有限体积法,（称为数值流通量）的含义,CopyrightbyLiXinliang,4,重要概念澄清：重构与插值,A.有限差分法：,j+1/2,切线,j-1/2,j,j-1,注意：与f在xj+1/2点的值含义不同！,用周围几个点的值计算的过程称为“重构”，不能理解为用来插值,记号确实容易混淆，让人容易联想起。记为更好些,否则，最高只能达到2阶精度了！,是控制体内的平均值,（称为数值流通量）的含义,CopyrightbyLiXinliang,5,重要概念澄清：重构与插值,B.有限体积法：,j+1/2,j-1/2,确实为f在xj+1/2点的值！,通常做法：1）用计算出2）,u在xj+1/2点的值！,关键：是用计算（称为重构），而不是用计算（是标准的插值）；否则最高也只能达到2阶精度。,6,概念：MUSCL与非MUSC类方法,j+1/2,切线,j-1/2,j-1,差分,有限体积,方法1（非MUSCL类）：直接利用周围几个点的函数值或）直接计算（或）,如何计算或？,方法2（MUSCL类）：利用周围几个点的自变量值（或）计算出（或）；然后再计算（或）,当f=f(u)是连函数时，二者精度相同,f的误差与u的误差同阶,CopyrightbyLiXinliang,7,10.1近似Riemann解,迎风型有限体积法通常需要求解Riemann问题；精确Riemann解计算量大,10.1.1HLLRiemann近似求解器（Harten,Lax&vanLeer）,思路：在控制体上积分，分析积分方程,在控制体上积分Euler方程,Ref.：E.F.Toro:RiemannSolversandNumericalMethodsforFluidDynamics,Springer,2009(ThirdEdition),控制体内质量（动量、能量）的减少等于流出控制面的通量,CopyrightbyLiXinliang,8,若控制体空间足够大（或时间跨度足够小），扰动波未达到控制体的边界（如图）,未扰动的值,边界上，未受扰动，保持常数,把积分域分成三段；中间为受扰动段,为扰动波传播的速度（激波或稀疏波的波头传播速度）,左、右波的速度假设已知,CopyrightbyLiXinliang,9,未扰动段，不影响积分，可不考虑,未扰动区，保持常值ULUR,消项,T时刻，扰动区内的平均值,T时刻的流动,以平均值代替分布值,CopyrightbyLiXinliang,10,思路：以扰动区内的平均值，代替瞬时值（扰动区很小，误差也不大）,含义：T时刻扰动区内的平均值,于是，HLL的近似Riemann解为：,HLL近似解是左右波都是激波，中心区域均匀分布的一个近似模型,双激波近似,假设：1）左右波都是激波2）中心区无接触间断（近似）,显然，假设中心区无接触间断是不合理的（造成方程数多于未知数个数），但可作为一种近似（类似最小二乘解）,CopyrightbyLiXinliang,11,利用HLL近似Riemann解计算中心线上的通量,在控制体（图中红色部分）积分Euler方程，得：,控制体积（图中绿色部分）积分,二者计算结果相同,带入表达式,CopyrightbyLiXinliang,12,最终界面的HLL通量如下：,原理：1）双激波近似，假设初始间断演化成两道激波，之间物理量为常数；2）利用积分关系式，确定了两道激波（左行激波与右行激波）之间的物理量；3）利用积分关系式，计算了穿过x=0界面的流通量,优点：计算简单缺点：没有考虑到接触间断。该现象（初始间断产生两道激波、中间没有接触间断）不可能发生。Riemann问题是三波问题，该近似为两波问题，假设过强。,CopyrightbyLiXinliang,13,10.1.2HLLCRiemann近似求解器（Toro）,发展了HLL近似解，用三波模型来近似（如图）,三波近似，左、中、右波的波速为，中间波为接触间断,HLLC近似解,T时刻的流动状态,CopyrightbyLiXinliang,14,根据积分关系，可知,红色区域积分可得,蓝色区域积分可得,假设中间波为接触间断：,利用积分关系计算接触间断的速度及其左右的物理量,如果把ZL,ZR也当成未知数，可精确求解（6方程6未知数，Riemann精确解），但计算复杂。,R-H关系式；弱解定义式,含义：控制体内质量的增加等于,控制体1控制体2,CopyrightbyLiXinliang,15,问题最终描述：求解超定方程4个未知数，6个方程,求解方程组：,其中：,假设波速ZR,ZL已知，未知数4个,与精确Riemann解的区别：本近似解假设左、右波速度已知，因此，未知数只有4个；精确Riemann解认为ZL,ZR也是未知数，因此6个方程，6个未知数,6个方程，4个未知数，怎么解？,6个方程,求解思路只使用其中的4个方程方案（1）：利用的连续性方程及动量方程（共4个），解出全部4个未知数,求解方案不能太复杂，复杂度不能超过求解原6个未知数的方程（Riemann精确解）；原则务必给出解析解,CopyrightbyLiXinliang,16,近似解：假设左、右波速已知，简化了计算（求解线性方程）精确解：左、右波速都是压力的函数，方程复杂（非线性），需迭代求解,线性方程动量积分关系式,求出,解出,利用压力相等，求出解出速度,利用质量及动量方程,CopyrightbyLiXinliang,17,最终HLLC（方案1）的表达式为：,HLLC近似解(方案1),CopyrightbyLiXinliang,18,HLLC近似解(方案2),根据方案(1)算出u*,p*,近似解方案（2）与方案（1）的区别方案（1）只使用了质量及动量积分关系（未使用能量积分关系）；方案（2）在最后一步也使用了能量积分关系,HLLC（方案2）给出的最终形式：,CopyrightbyLiXinliang,19,HLLC近似解（方案3）,假设压力时左、右两个解的平均,方案（1）令这两个解相等，以此计算Z*方案（3）不要求他们相等，而以平均值代替,最终形式（近似解方案3）：,6个方程，4个未知数，近似解方法很多，希望读者思考，开发出更好的方案,CopyrightbyLiXinliang,20,10.1.3波速的估算,HLL及HLLC均假设ZL,ZR已知，实际上它们仍需要估算,准确计算ZL,ZR，实际是计算Riemann精确解，计算量大,方法1：直接估算,1a:假设以声速传播（Davis）,竟然假设激波以声速传播，太OUT了,小常识：激波的传播速度激波相对于波前介质以超声速传播，相对于波后介质以亚声速传播；弱激波（Ma趋近于1）以声速传播。,1b:左、右两种状态声速的平均(Davis,Einfeldt),要平均吗？用Roe平均吧，效果可好了。,激波速度介于波后（相对）声速与波前（相对）声速之间，平均是个好思路,CopyrightbyLiXinliang,21,Roe平均：,1c.Roe平均的修正（Einfeldt）,方法2：基于压力的波速估算法（Toro）,已知中心区压力，容易计算波传播速度,中心区的估算：,CopyrightbyLiXinliang,22,10.1.4“扩展一维Riemann问题”的HLLC近似解,特点：切向速度v,w的处理方法与被动标量相同；被动标量场穿过激波值不变,CopyrightbyLiXinliang,23,10.1.5HLLC方法计算通量的步骤,Step1：估算左、中、右波的速度（选用压力估算，还可采用直接估算）,Step2:计算HLLC型通量（选用近似方案1，还可选用其他两种方法）,左、右激波的速度,中间波（接触间断）的速度,中心区压力的估算,CopyrightbyLiXinliang,24,在有限体积法中的具体使用方法：,控制体上积分，得到积分方程,j-1jj+1j+2,j-1/2j+1/2,Step1:重构，用计算,相当于差分方法的差分格式；方法非常丰富,例如：,Step2:求通量,把代人上页的公式，得到：,相当于差分法的通量分裂（属于FDS类），方法很多,CopyrightbyLiXinliang,25,10.2中心型有限体积法,中心型及迎风型方法的对比,为了使计算稳定，中心型格式通常添加人工粘性,CopyrightbyLiXinliang,26,j-1jj+1j+2,j-1/2j+1/2,控制体上积分,采用两侧的值重构j+1/2点的通量，例如：,非MUSCL型,MUSCL型,注意：不要把重构理解为插值，否则精度不会超过2阶,添加人工粘性,二阶人工粘性,四阶人工粘性,CopyrightbyLiXinliang,27,Jameson人工粘性,j-1jj+1j+2,j-1/2j+1/2,二阶和四阶人工粘性,光滑区：间断区,经验常数，可调整人工粘性大小,CopyrightbyLiXinliang,28,10.3具体算例讲解RAE2822翼型绕流,算例：RAE2822超临界翼型绕流的数值模拟,RAE2822翼型及压力测量装置,问题描述：二维跨声速流动吹过RAE2822翼型，试计算流场及翼型表面的压力分布。,流动参数：Mach0.729,攻角（AoA）=2.31，Re=6.5106（基于弦长及来流值）,背景：NASANPARCAllianceVerificationandValidationArchive的标准算例（Bechmark）之一，有可靠的实验数据、网格及多个数值结果。成为目前CFD程序验证的标准模型之一。,网站：/WWW/wind/valid/validation.html（用Google搜索RAE2822即可）,CopyrightbyLiXinliang,29,有限差分法的计算结果5阶及7阶迎风差分+Steger-WarmingFVS+3阶Runge-KuttaBL湍流模型,网格：NPARC网站提供，369*65,流动为湍流，二维计算需使用湍流模型；湍流模型将在后续课程介绍，目前先按层流处理。并假设流动定常,CopyrightbyLiXinliang,30,1.控制方程、网格及边界条件控制方程无量纲N-S方程,无粘和粘性项,计算网格,边界条件：外边界：无穷远来流，出流条件固壁：无滑移绝热边界,无穷远来流条件,出流条件,网格：从NASA网站下载,CopyrightbyLiXinliang,31,边界条件的处理方法：,方案1：严格按照特征方向（精确，但复杂）,外边界1,外边界2,边界,对于每个边界单元上讨论特征方向，在法方向局部一维无粘化,边界处的物理量,在边界法向的坐标系中，控制方程为近似为“扩展的1维Euler方程”,求解，计算边界值，单边差分（重构）,CopyrightbyLiXinliang,32,计算模型：4个特征波，分别以速度在该一维直线上传播，根据传播方向独立讨论。若传播方向指向边界内，则给定流动参数；否则，用内部点计算,边界,方案2：近似处理,外边界1：流场远离翼型，扰动很弱，认为（无论亚、超均按超音速处理）外边界2：按超音速出流处理,外边界1,外边界2,近似模型：,对于内流，有问题,固壁边界条件,物理条件：,补充计算条件（常用）：,CopyrightbyLiXinliang,33,2.空间离散,无粘通量,粘性通量,控制体上积分：,注：选取控制体的方法很多,本计算采用的控制体节点中心型,控制体的多种选择方式,CopyrightbyLiXinliang,34,2.空间离散,无粘通量,粘性通量,（1）无粘通量的计算,旋转到xy坐标系,（迎风）重构控制体边界处的左、右值,线性重构：平均值=中心点值,求解“扩展”Riemann问题，得到控制体边界m处的通量,建议：使用HLL或HLLC型通量（公式见本PPT12、17页）,q是x轴与x轴的夹角,CopyrightbyLiXinliang,35,具体步骤:无粘通量,讨论以i,j为中心的控制体,在所在的控制体表面（记为m）上,1）重构自变量的左、右值,2）进行坐标旋转，得到坐标系（x,y）下的自变量,q是x轴与x轴的夹角,3）求解Riemann问题，获得界面m上的通量,m,使用HLL或HLLC型通量，见12，17页,4）在其他3个控制面（图中虚线所示）上，同样求出通量,5）4个面的通量求和，得到该控制体的总无粘通量,CopyrightbyLiXinliang,36,（2）粘性通量的计算采用中心格式,可用周围角点的平均,关键问题：如何计算？,方案1：借鉴有限差分法，引入坐标变换,具体表达式见:任玉新等计算流体力学基础130-131页及127页,转化到计算坐标系,通常，令,计算坐标实际上就是下标（i,j）,（1）,CopyrightbyLiXinliang,37,于是：,坐标变换系数的计算：,控制体的体积,半点上的Jocabian系数（体积加权）插值得到,完全的有限差分法，坐标变换