第7讲-差分方法3

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,Copyright by Li Xinliang,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,Copyright by Li Xinliang,*,知识回顾,1,）,守恒型方程与守恒型格式,守恒型方程： 散度型,守恒型格式： 差量型,仅为记号，,与,j+1/2,点上的值无关,！,守恒型方程,+,守恒型格式,=,“解守恒”： 数值解的积分误差为,0,（如果边界准确）,保证总量（总质量、总动量、总能量）严格守恒 （无误差）,早期的,CFD:,极为重视守恒性,；,目前,CFD:,仍很重视守恒性,难点,复杂非线性系统的守恒性很难保证,中间项全部消去，只剩两端,2,）,通量分裂,便于使用迎风格式,方法（,1,）： 逐点分裂,（,Steger-Warming, Van Leer, L-F,）,原理： 利用了性质,使得,的,Jocabian,阵特征值,纯正或纯负,优点： 无需矩阵运算，计算量小，使用方便,不足： 仅重新组合，没有做到真正解耦。 原因：,具体方法：,Steger-Warming,L-F,Van Leer,or,方法,2,） 特征（投影）分裂,在网格基上冻结系数, j-2 j-1 j j+1 ,在基架点上系数 不变,优点：特征分解，（局部）解耦,耗散小，数值振荡低,缺点：大量矩阵运算，计算量大,每计算一个点的导数，要进行,m,次矩阵运算 （,m,为网格基点数）,原理描述 （非守恒，很少采用； 实际上使用下一页的方法）,守恒型方式,计算, j-2 j-1 j j+1 ,在基架点上系数 不变,具体步骤：,假设已知,U,， 且针对模型方程（线性单波方程）,已构造出差分格式,（,1,）,1,） 计算出,教材,130,页的公式,(6.1.11-6.1.13),式中用到各变量在,j+1/2,的值（例如 ）,可使用,j, j+1,点值的算术平均 （如 ）,或,Roe,平均,；,由 计算；方法很多，例如前面介绍的 或,均可,推荐使用,Roe-,平均！,2,）,在网格基上计算, j-2 j-1 j j+1 ,计算,f,j+1/2,用到的点,注意，在该网格基上（例如,k=j-1,j,j+1,） 保持不变,例如：,3),利用已构造好的差分格式，计算通量,4,） 得到总通量,5,） 计算差分 （,j,点处）,步骤的算法描述,（注意： 实际上是两重循环）,do j=1,N,do k=j-1,j+1 (,网格基，可以是更多或更少点,),enddo,enddo,do j=1,N,enddo,需要多次矩阵运算，计算量大,守恒性好，耗散小，数值解质量好,通量分裂,+,迎风差分,引入数值耗散,分裂本身不带来耗散， 但会放大（或减少）差分的耗散,举例：,分裂过程,耗散,如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。,=,+,向上平移,向下平移,分裂后的流场越偏离原先流场，则总体耗散越大,如使用低精差分度格式， 则对分裂形式敏感 （推荐使用特征分裂）,如使用高精度格式（低耗散），则对分裂形式不敏感 （可使用逐点分裂）,概念澄清,耗散放大系数, 7.,1 Roe,格式,守恒型格式的范例,破坏守恒性,后果很严重（,?,）,为了便于使用迎风格式、特征分裂解耦， 通常把守恒型方程改写为非守恒型,守恒方程,+,守恒格式,=,解守恒,方程不守恒， 即使差分方法守恒，也无法做到解守恒,由于,a,非常数，无法消去中间项 ！,不再守恒,思路： 保证特征方向，找回守恒性,守恒型方程优点： 守恒性,非守恒方程优点： 清晰的特征方向,兼顾守恒与非守恒方程,1.,单方程的,Roe,格式,1,阶迎风（直接从守恒方程出发）,（,1,）,（,2,）,体现了特征方向,只有这种表达式，,才能保证 （,2,）与（,1,）等价,（,3,）,or,都无法保证 （,2,） 与（,1,）等价。,简单的线性平均不行,(,非线性系统，中点的斜率不等于平均斜率）,关键： 构造,Roe,格式,“平均斜率”，不等于“斜率的平均值”，也不等于中点处的斜率,平均斜率,2.,方程组的情况 （,Roe,格式的意义）,需满足如下条件 （,Uniform,特性）,单方程的简单推广,1,） 连续，且,2,） 可通过相似变换对角化，即,保证双曲性,3,）对于任意,U,R,U,L,有,单方程的推广，含义为,平均增长率,标量方程向矩阵方程的简单推广，,但 构造很困难。,均不满足,Uniform,特性,平均斜率,3.,矩阵 的构造,关键：,“向量除以向量” ？,直接求平均增长率：,u,f(u),u,L,u,R,u,Roe,Roe,点的斜率为平均斜率,（根据拉格朗日中值定理，,U,L,U,R,区间内肯定存在,Roe,点）,思路,1,： 在,U,L,与,U,R,之间寻找一个点,U,Roe,该点处的增长率为平均增长率,f(u)=u,2,u,二次函数, Roe,点与中点重合,标量函数的启示：,Roe,点肯定存在（,Langrage,中值定理）,二次函数的中点即为,Roe,点,思路,2,： 进行坐标变换，得到一个二次（齐）函数,引入,如果,是二次（齐）函数，则其中点 即为,Roe,点,重要启示,更准确地讲，应当是要求 为,W,的线性函数，,即增长率为线性函数 （中点处的增长率刚好为平均增长率）,针对,Euler,方程的具体构造,引入新变量：,则：,目的： 使得,F(w),是,W,二次齐函数 （增长率为线性函数）,f(U),不是,U,的二次齐函数,二次齐函数！,中点处的斜率即为平均斜率。,Roe,点,Roe,点为：,增长率为线性函数！,最终：,其中 如下计算：,平均增长率（矩阵）,含义：,左、右两个状态点的某种平均 （称为,Roe,平均，为密度加权平均）,该状态点对应的增长率（矩阵）为平均增长率（矩阵）,实际上是一种“等效平均”。 效果优于简单的算数（或几何）平均。,三维情况下，还有,其他量（如压力、温度、音速等）用这三个量计算,（,5,）,简单易记：,Roe,格式的计算步骤 （半离散）,已知,n,时刻所有网格点上的物理量，对于,j,点：,1,） 令,U,R,=U,j+1,U,L,=U,j,（密度、压力、速度等）,2,） 采用,Roe,平均公式（,5,）计算,Roe,平均值,3,） 将,Jacobian,矩阵 进行特征分解,:,计算,4),计算,5,）计算,6,） 计算空间导数,7,）时间推进，计算下一时间步的值。,j-1 j j+1,与前文（第,3,，,4,讲）的形式相同，仅需把式中的密度、压力、速度等换成经过,Roe,平均的密度、压力、速度即可,其中：,可能出现导数不连续， 可能引起数值振荡,实际使用时 可用如下函数代替,所谓“,熵修正”,实际上是在特征值,0,点周围增加了耗散,Roe,格式的优点：,1,） 保持守恒性的同时，严格保证了特征方向,2,） 便于推广到高精度格式,特征投影分裂中使用,Roe,平均即可 （见本,PPT,第,5,页）。推广到高阶后，虽不再保证严格的特征方向， 但仍优于采用算数平均方法。,Roe,格式的不足：,本身精度只有一阶；,推广到高阶后，特征方向无法严格保证 ；,推广到二维或三维后，特征方向无法严格保证，出现振荡。,作业,7.1,对于一维,Euler,方程：,引入新变量：,推导出 及其,Jacobian,矩阵 的具体表达式（以,W,为自变量），并证明对于任意 ，有：,提示： 写出 表达式后，将向量 分别代入上式左、右两端，容易证明相等。,要求：推导过程要详细，切勿简单从书本上摘抄。,重要的,CFD,基本功练习，一定要重视！,针对如下,Sod,激波管问题,用,Roe,格式,计算其数值解，画出,t=0.14,时刻密度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。,要求：,1,） 空间网格数,100,， 时间推进格式选用,3,阶,Runge-Kutta,时间步长自选。,2,） 尝试使用熵修正与不使用熵修正两种情况（见本,PPT 15,页）,3,） 欢迎与其他数值方法得到的结果对比（最好画在同一张图上，便于比较）。,作业,7.2,7.2,非物理振荡及,TVD,格式,1.,数值解中的非物理振荡,间断附近非物理振荡的根源,理论,1,：色散误差导致各波传播速度不同 （第,4,讲）,理论,2,：物理粘性的错误计算,思路： 物理问题,有粘； 物理粘性足以克服本身振荡,数值方法错误计算了物理粘性,不足以克服振荡,物理问题本身也可能振荡。但如果错误计算物理粘性，则会错误地加剧（或衰减）振荡。,1,） 非物理振荡的原因分析,数值实验,二阶中心差分,计算域,0,1,网格点,201,（,D,x=0.005,）,时间步长,D,x=0.0005,T=0.1,时刻的,u,分布,Re=200,D,x=0.005,现象：,D,x,一定时，减小,Reynolds,数可抑制振荡,Reynolds,数一定时，减小,D,x,可抑制振荡,暗示,是某一特征量,Re=2000,D,x=0.005,Re=2000,D,x=0.0005,相同,对流,-,扩散方程的特性：,n,n+1,差分方程：,某点的值是上一时刻周围几个点上值的线性组合,物理上要求系数,a,k,均非负,含义： 某处浓度的增加对下一时刻周围浓度的影响为正。,j-2 j-1 j j+1 j+2,差分方程单调性（无振荡）条件： 差分方程 （,1,）中的系数非负,网格,Reynolds,数,2,） 重要概念： 网格,Reynolds,数,以网格尺度度量的,Reynolds,数,含义：,数值振荡,流动尺度为网格尺度,网格,Reynolds,数小，该尺度的能量被耗散掉,不发生振荡,j,j+1,j-1,过于苛刻的条件,单方向网格点数,10,6,， 三维,10,18,单纯靠物理粘性抑制振荡，网格间距必须足够小，通常难以实现,网格足够小：不会发生振荡；,网格小于激波的实际厚度，则不会振荡,网格,Reynolds,数足够小时，物理粘性发挥作用，抑制振荡,3,）,人工粘性,物理粘性 足够小才发挥作用，,Reynolds,数很高时很难做到,思路： 人为增加粘性系数 （添加人工粘性） 抑制振荡,优点：方法简便，有抑制振荡效果,缺点：改变了物理问题，带来误差,湍流、分离流等,对粘性敏感： 非物理解,分离流,对粘性敏感,转捩,对粘性敏感,很难计算对粘性敏感的问题,改进措施：,A:,局部施加人工粘性,B:,高阶人工粘性,Von Neumann,MacCormack,4),数值振荡的定量描述,总变差,对于离散函数,uj,定义总变差,:,单调函数,振荡函数,j=1,j=N,含义： 反映了振荡的剧烈程度,双曲型守恒方程,特点： 沿特征线 ，,u,不变,特征线未相交,总变差不变,特征线相交,总变差减小,结论： 单个双曲型方程，总变差不增,（,Total Variation Diminishing,：,TVD,）,2,概念： 单调格式、保单调格式与,TVD,格式,n,时刻： 单调函数,j=1,j=N,n+1,时刻： 仍是单调函数,j=1,j=N,设,n,时刻 是单调的，如果,n+1,时刻的解 仍保证单调，则称该格式为,保单调格式,。,保单调格式,基本结论：,常系数的单调格式只能是一阶 （,Godunov,）,单调格式必是保单调的；,线性格式，单调与保单调等价,格式：,如果满足 则称其为单调格式,。,单调格式：,单调格式,保单调格式：,TVD,格式,总变差不增,TVD,保单调,单调,3. TVD,格式的理论基础, Harten,定理,Harten,定理：,如果差分格式可写成如下形式：,且,则格式（,1,）是,TVD,格式,（,1,）,可验证：,Roe,格式是,TVD,格式,保证“系数非负”,含义：“单调格式必是,TVD,格式”,例,7.2.1,：,考虑线性单波方程,：,试讨论如下,Lax-Wendroff,格式,二阶中心,人工粘性,是否满足,Harten,条件,单调格式,只有一阶精度,对比条件：,不满足,Harten,条件,知识回顾：,Lax-Wendroff,格式,Taylor,展开，写出修正方程,时,-,空二阶精度,巧妙添加人工粘性，不但克服了不稳定性，而且抵消了时间误差，提高了时间精度,类似方法：,Beam-Warming,格式,人工粘性,二阶精度迎风差分,人工粘性，且提高时间精度,特点： 全离散、时刻耦合,4.,构建,TVD,格式,思路： 对现有格式进行改造，使之符合,Harten,条件,通常在,Roe,、,L-W,、,B-M,（或其组合）基础上改进,80,年代初、这些格式是主流,（,1,） 以,L-W,格式为基础改造的格式,L-W,原格式（,2,阶）,= 1,阶迎风,+,修正项,新格式,= 1,阶迎风,+,限制函数,*,修正项,引入限制函数 （限制器）,1,阶迎风部分,修正项,显然 格式为,LW,（,2,阶）,可验证： 格式为,B-M,（,2,阶）,新格式,= 1,阶迎风,+,j,*,（,LW,格式,-1,阶迎风）,新格式：,LW, BM,均为线性格式，二者组合仍为二阶,根据,Harten,定理，可知,时，可满足,TVD,性质,(2),精度条件,Beam-Warming,二阶精度区,TVD,区,二阶精度,TVD,区（二者交集）, 7.,3 WENO,格式,高精度的激波捕捉法,1.,基本思路,j-3, j-2,j-1,j,j+1,j+2,j-3,j-2,j-1,j,；,j-2,j-1,j,j+1,；,j-1,j,j+1,j+2,五个基架点被分成三个组,1,） 若高精度逼近 ， 必然利用多个基架点,2,） 如果该基架点内函数有间断，会导致振荡,3,） 间断不可能处处存在,4,） 把基架点分成多个组（模板），,每个模板,独立,计算,j,点导数的逼近。,得到,多个差分,5,）根据每个模板的光滑程度，设定权重,6,） 对多个差分结果进行加权平均 。光滑度越高，权重越大。如果某模板存在间断，则权重趋于,0,；,如果都光滑，则组合成更高阶格式。,2. WENO,格式的原理描述,考虑线性单波方程：,注： 为了简便，以非守恒型形式为例讲授其思路，实际使用时，请采用下一节介绍的守恒形式,（,1,） 确定网格基架点：,6,个点,j-3 , j-2,j-1,j,j+1,j+2,构造出该基架点上的目标差分格式,计算,这,6,个点可构造,5,阶迎风差分：,该格式为,WENO,的“目标”格式，,即， 光滑区,WENO,逼近于该格式。,利用,Taylor,展开，可唯一确定系数 （可利用小程序,coeff-schemes.f ),实际上，还可利用分辨率优化技术，可构造出新的目标格式（降低精度、提高分辨率，见第,4,讲）。目前大量,WENO,的优化版做这种工作。,将这,6,个基架点分割成,3,个组（称为模板）,每个组独立计算 的差分逼近,模板,1,模板,2,模板,3,模板,1,：,j-3,j-2,j-1,j,模板,2,：,j-2,j-1,j,j+1,模板,3,：,j-1,j,j+1,j+2,利用这,三个,模板的基架点，可构造出逼近 的,3,阶精度差分格式,计算,j,点的导数,u,竟然算出了三个不同的值，怎么办？,ENO,方法： 选择最优（最光滑）的，舍弃其余两个,WENO,的处理方法： 三个都要，加权平均它们。,利用,Taylor,展开式，可唯一确定这些系数）（可利用小程序,coeff-schemes.f ),也可运用优化技术，降低精度、提高分辨率,（,3,） 对这,3,个差分值进行加权平均，得到总的差分值,原则：,A.,模板内函数越光滑，则权重越大； 模板内有间断时，权重趋于,0,B.,三个模拟内函数都光滑时，这,三个三阶精度,的逼近式可组合成,一个五阶精度,的逼近式。,“理想权重”,（,3.1,） 确定理想权重,令：,5,阶精度,容易解出：,（,3.2,） 度量每个模板内函数的光滑程度,IS,越大，表示越不光滑。,光滑区，不同模板上的,IS,趋近同一值。,具体形式见下一节。,（,3.3,） 给出实际权重,构造,IS,方法很多， 例如：,：第,k,个模板,光滑区逼近,O,（,1,）量级,间断区 量级，很大,特点： 间断区权重很小,光滑区，趋近于理想权重,（,3.4,） 给出最终的差分逼近,3. Jiang & Shu,的五阶,WENO,格式,守恒型；目 前使用的,WENO,格式均为守恒型,针对方程：,模板,1,模板,2,模板,3,构造差分格式如下：,构造方法与前文相同 （但注意这里构造的是通量，而前文是直接构造差分格式）,针对整个网格基，构造出,5,阶精度的通量（理想情况下的通量）,并构造出每个模板上的通量，计算出理想权重。,仍利用程序,coeff-schemes.f,求系数,理想权重,光滑度量因子,实际权重,光滑度量因子的计算 （,Jiang & Shu,）,k=1 k=2 k=3,其中：,j-2 j-1 j j+1 j+2,是使用模板,k,得到的插值函数,利用,j-2,j-1,j,点上的值构造的插值函数,，,特点： 光滑区趋近同一个值,非光滑区值远大于光滑区,O(1),j,点一阶、二阶导数的差分逼近（用模板,k,计算）,代入,最终,5,阶,WENO,格式为,正通量情况（,a0,）,负通量情况 （,a0,）的方法计算,采用针对负通量（,a0,）的方法计算,可采用,Steger-Warming, L-F, Van Leer,等分裂， 见第,4,讲,1),计算,，它们是 的函数（推荐使用,Roe,平均计算 ）,效果更好,但计算量较大,计算量小,但效果略差,有轻微振荡,数值测试：,1,维,Sod,激波管问题，网格点,128,差分方法,: 7,阶精度经典,WENO (Jiang & Shu),分裂方式：,Steger-Warming FVS;,特征投影分裂（算术平均,/Roe,平均）,t=0.1,时刻密度及速度分布,Steger-Warming FVS + WENO,仍无法避免振荡,特征投影分裂,+WENO,可避免振荡,39,针对如下,Sod,激波管问题,用,5,阶,WENO,格式,计算其数值解，画出,t=0.14,时刻密度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。,要求：,1,） 空间网格数,100,， 时间推进格式选用,3,阶,Runge-Kutta,时间步长自选。,2,） 结合使用,Steger-Warming,流通矢量分裂，画出结果,3,） 结合使用特征投影分裂（,Roe,平均计算,U,j+1/2,),，画出结果,作业,7.3,针对如下,Shu-Osher,激波,-,密度扰动波干扰问题：,用,5,阶,WENO,格式,计算其数值解，画出,t=0.1,时刻密度、速度及压力的分布,要求：,1,） 空间网格数,200,， 时间推进格式选用,3,阶,Runge-Kutta,时间步长自选。,2,） 结合使用,Steger-Warming,流通矢量分裂，画出结果,3,） 结合使用特征投影分裂（,Roe,平均计算,U,j+1/2,),，画出结果,4,） 使用,2000,个网格点计算，其结果作为“精确解”，与其它结果画在一起，便于比较。,作业,7.4,初值为,:,作业,7.5,（选作题）,推导,7,阶精度的,WENO,格式，给出详细的推导过程及格式的具体表达式,提示： 与,5,阶,WENO,推导思路相同，但网格基架点扩大了，模板数目也增加了,j-3 j-2 j-1 j j+1 j+2 j+3,正通量,(a0),的,WENO,通量 使用基架点,j-3,j-2,j-1,j,j+1,j,j+2,j+3,如上图示，将其分割为,4,个组（模板），每个模板上,4,个基架点。,1,） 先构建整个网格基（,7,个点）上,7,阶迎风格式的通量表达式,2,） 对于每个模板，构造逼近,j,点导数的,4,阶差分格式的通量表达式,3,） 按照理想权重进行组合可得到理想格式（,7,阶迎风格式）,对照具体表达式，求出理想权重,C,k,4,） 构建,WENO,通量的加权表达式,5,）仿照 本,PPT,第,35,页的方法，给出光滑度量因子的具体表达式,可利用小程序求系数,coeff-schemes.f,减小工作量,