第10讲有限体积法2

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