计算流体力学第讲差分方法

1、 数学楼数学楼315； 52377625 email: 计算流体力学讲义计算流体力学讲义 9/2/20211 *21. 耗散与色散误差耗散与色散误差2.半离散分析与全离散分析半离散分析与全离散分析 -如何理论计算修正波数1. 耗散与色散误差耗散与色散误差3Copyright by Li Mingjun数值实验数值实验 boundaryperiodicxxxuxutu,2 , 0)sin()0 ,(0 时间推进时间推进： 3步步TVD型型Runge-Kutta, 且时间步长足够小且时间步长足够小 （误差忽略）（误差忽略） 空间离散空间离散： 1阶及阶及2阶迎风格式阶迎风格式 （20个网格点）个网

2、格点）实验观察到的现象实验观察到的现象 两类误差：两类误差： 振幅误差振幅误差 相位误差相位误差 （波速误差）（波速误差） 9/2/202134Copyright by Li Mingjun)sin()0 ,(0 xxuxutu精确解精确解1阶迎风阶迎风2阶迎风阶迎风xuuujjj / )(1)2/()43(21xuuuujjjj )sin(),(txtxu 9/2/202145Copyright by Li Mingjun)sin()0 ,(0 xxuxutu 9/2/202156对以上对以上“实验现象实验现象”进行理论分进行理论分析析考查问题：精确解：ikxikctctxikeeetxf

3、)(),(差分格式：0 jxjfctf （1）xfffjjjx / )(1 )2/()43(21xffffjjjjx 其他格式 jikxjfe假设对于：假设对于：有有jikxxjkfex隐含假设隐含假设: 线性差分格式，线性差分格式， 非线性非线性 系统作用于单波，系统作用于单波， 会产会产 生多个谐波生多个谐波（2） 差分没有误差差分没有误差xikk Copyright by Li Mingjun)0(0 cxfctfikxexf )0 ,(0,2 ,xperiodicboundary1阶迎风格式：2阶迎风格式： 9/2/202167半离散分析：半离散分析： 假设时间推进是精确的，仅分析空间

4、 离散带来的误差（难度小、常用） 全离散分析：全离散分析： 同时分析时、空离散的误差 （难度大） 注意：注意：实际上就是普通三角函数，采用复数形式仅仅是实际上就是普通三角函数，采用复数形式仅仅是为了理论推导方便。若用实数形式为了理论推导方便。若用实数形式 sin(kx), cos(kx)推导推导形式上略显繁琐。形式上略显繁琐。Copyright by Li Mingjunikxikctctxikeeetxf )(),(精确解：精确解：2.半离散分析与全离散分析半离散分析与全离散分析 -如何理论计算修正波数如何理论计算修正波数 9/2/202178令：jikxjetutxf)( ),( jikx

5、jxxektutxf /)( ),( （1）式化为：0/)( )( jjikxikxxektucedttud0/)( )( xktucdttud“半离散化半离散化”： 空间导数差分 计算，时间方程 （常微）精确计算xctkeutu /)0( )( xctkikxikxjjjeuetutxf /)0( )( ),( ixikk 如果如果 ， 无误差！无误差！Copyright by Li Mingjun设：()( , )ik x ctf x te精确解： 9/2/202189分析分析 （修正波数）与误差的关系：（修正波数）与误差的关系：kirikkk )/(/)0( ),(xkctkxikxct

6、kjijreeutxf 的误差导致解的幅值误差 耗散误差耗散误差 的误差导致解传播速度的误差 色散误差色散误差xkkkir , 0rkikCopyright by Li Mingjun (理想情况理想情况： )设设则则 的实部：的实部：耗散耗散误差误差 的虚部：的虚部：色散误差色散误差k 9/2/2021910jikxjexu )(jikxjxexku 含义：含义： 反应波数（谱）空间内差分的误差jikxkkjefxf )(任意函数：任意函数：jikxkkjefikxf )(定义：定义：设设 ，kkfikf 求导数，精确解差分解jikxkkjxjefxk ifF / kkfxkF Copyri

7、ght by Li Mingjun修正波数修正波数 定义 为：k 9/2/202110 反映了一个波内的点数。反映了一个波内的点数。PPW （波内的点数）（波内的点数）= 11 计算出计算出 ，并考差其与，并考差其与 的逼近程度。的逼近程度。k ixik 考察格式分辨率（考察格式分辨率（resolution）的重要指标）的重要指标0 xxk 优秀的差分格式，优秀的差分格式，1个波长里面个波长里面6个点个点 即可即可 /2Copyright by Li Mingjun精度精度 分辨率分辨率Fourier 分析的任务：分析的任务： 精度：精度： 反映 时的情况 分辨率：分辨率：网格点数很少（例如波

8、里面只有6个点）时的 性能， 对于多尺度问题，分辨率更重要。 （牺牲精度， 提高分辨率） 9/2/20211112如何理论计算如何理论计算修正波数？修正波数？方法方法：根据差分具体表达式及定义计算例例1：xuuujjjx 1 令,jikxjeu 则：jjjjikxiikxxxikikxjjjxexkexeeexxuuu )1()(1)(1 xk 其中 。于是： sin)cos1(1ieki 1阶迎风格式： sin,cos1 irkkCopyright by Li Mingjun 9/2/20211213例例2：)2/()43(21xuuuujjjjx 2阶迎风格式：)43(2243221 ii

9、ikxjjjjxeexexuuuuj 2/ )43(2 iieek 2/ )2sinsin4(,2/ )2coscos43( irkk令jikxjeu 则：Copyright by Li Mingjunxk 其中 。 9/2/202113 *144.1 Jacobian 系数矩阵及其性质系数矩阵及其性质4.2 对流项的分裂对流项的分裂 (2) 严格特征分裂严格特征分裂 基架点上冻结系数基架点上冻结系数(1) 两类典型的逐点分裂两类典型的逐点分裂0 xuatu格式F+格式F-visf(U)U xt2( ,) ,(, ()TTUu Efuup u Ep x f(U)xua 对流项：对流项：信息（波

10、）从上游传至下游 ，上游更重要，迎风差分扩散项：扩散项： 信息从中心向周围扩散，不区分上、下游，中心差分迎风差分优点：迎风差分优点： 有效利用信息传播的方向，增强稳定性N-S方程：单波方程：单波方程：单波方程： 一个波，容易判断波传播方向N-S对流项（对流项（Euler）：）：方程组，多波问题， 复杂15?0?,0 aaCopyright by Li Mingjun4.1 Jacobian 系数矩阵及其性质系数矩阵及其性质 9/2/2021150 xuatu单波方程单波方程：16Copyright by Li Mingjun对于正数对于正数a,说明波沿说明波沿x正方向传播，如下图所示：正方向传

11、播，如下图所示：10iiuuuatx采用空间向后差分格式采用空间向后差分格式：何为迎风格式？何为迎风格式？11()nniiiituuauux 9/2/202116微分与差分方程的影响域双曲方程组的原则：双曲方程组的原则： 特征分解，找到独立传播的波00VV0USSU0UAU1 xvtvxtxtxtkkk 常系数矩阵常系数矩阵A的情况：的情况：完全解耦，独立求解变系数矩阵变系数矩阵A的情况：的情况： 局部讨论17Copyright by Li Mingjun(1) 常系数矩阵常系数矩阵A两类典型系数矩阵两类典型系数矩阵A： 9/2/202117(2) 一般一般Jacobian 系数矩阵系数矩阵A

12、,UAf(U)xx ,Uf(u)A f(U)U( f f(U)U)U()U)f( 1 AUUU)f(U)(f(U) xxx UA)AU(f(U)重要性质 uucucuuu 222322311221)3(23010Attt UA)AU(f(U)特点：特点： A 可以像常数一样，和求导运算交换visf(U)U xtTTpEupuufEuU)(,(,),(2 18Copyright by Li MingjunN-S方程方程： 9/2/2021184.2 对流项的分裂对流项的分裂 方法：方法： (1) 逐点分裂， (2) 严格特征分裂 (1) 两类典型的逐点分裂两类典型的逐点分裂（利用性质)AU(UA

13、f(U)xxx ,SSA1 ,kdiag , kdiag ,22kkkk=+优点：优点：耗散小缺点：缺点：导数间断SSA1,1fA US SU2)(2/122 kkkA: Steger-Warming 分裂分裂19Copyright by Li Mingjun目的：目的： 确定波传播方向，便于使用迎风差分或者：或者： 9/2/202119 Steger-Warming 具体步骤具体步骤 （以一维为例）（以一维为例）已知已知TEuU),( 1） 计算计算2） 计算计算3） 计算计算4） 代入（代入（1）式得到）式得到5） 利用不同的迎风格式，分别计算利用不同的迎风格式，分别计算pu, cucuu

14、 321, 2)(2/122 kkk)3 , 2 , 1(, kkk wcucuucucuu232221321321)(2)(2)1()()()1(2)1(22)(f )1(2)(3(232 cw)(ff),(ff （1） f ,f fffxx f,fxx ff0 a0 xuatu0 a（后差，前差）（后差，前差） 6） xxx fff计算计算7） 时间推进时间推进0fU xt20Copyright by Li Mingjun 9/2/202120二维问题的二维问题的steger-Warming 分裂分裂0ffU21 yxt,)(,f21TpEuuvpuu ,)(,f22TpEvpvuvv ,

15、UTEvu 21fff令：则： WvuvuVvvvuuu2222)(f2222421213202413024130430 )1(2)(3(243 cW10)1(2 11kcuu 12kcuu /1 k /2 k21kcvv 22kcvv 22 具体使用步骤，具体使用步骤， 以计算以计算 为例为例x 1f令令 计算特征值计算特征值 分裂特征值，计算分裂特征值，计算 带入左式，计算正、负流通矢量带入左式，计算正、负流通矢量 计算计算0, 1 cucuu 4321, 2)(2/122 kkk)4 , 3 , 2 , 1(, kkk )(ff1 xxx 111fffy 2f计算计算 设置设置 ，并注意

16、，并注意1, 0 cvcvv 4321, 对于曲线坐标系对于曲线坐标系0ffU21 t21111fffyxJJ 仅需令yxJJ 11, 21Copyright by Li Mingjun 9/2/202121B: L-F分裂分裂AUf(U) 2/ ) IA(A, 2/ ) IA(A* 特点：特点： A 正特征值 负特征值 A2/ )Uf (f* UA=+缺点：缺点：耗散偏大局部局部L-F分裂分裂，每个点上计算，每个点上计算 全局全局L-F分裂分裂，全局（一维）上计算，全局（一维）上计算 cu * * 足够大cu * 数学性质（光滑性）数学性质（光滑性）最好，但耗散偏大最好，但耗散偏大)(max

17、*cux 常数例如，可取例如，可取22Copyright by Li Mingjun 9/2/202122UfU)(2f )(21*0*xxxxxxx 0f(U)U xt22*Uf(U)Uxxt 23Copyright by Li Mingjun)Uf ()Uf (21fff* xxxxx（人工粘性人工粘性）)Uf ()Uf (21* xx 与迎风格式结合，等价于人工粘性与迎风格式结合，等价于人工粘性 9/2/2021232,2kkkkkk kkk 方式很多，典型的有方式很多，典型的有3种种=+S-W:L-F:=+2,2* kkkk0, 0 kk Van Leer:=+ 221112)1()1

18、(22)1(cufcufff 2121 Macf ff 0 f1 Ma1 Maff cuMa/ 0 f1 Ma24Copyright by Li Mingjun 9/2/20212425Copyright by Li Mingjun优点：优点： 无需矩阵运算，计算量小缺点：缺点： 分裂后改变了特征方向， 耗散大xxx (U)f(U)ff(U)SUSUA(U)f1 利用性质 )AU(UAf(U)xxx 逐点分裂的特征逐点分裂的特征 9/2/202125 分裂后分裂后 失去了失去了A的性质（可以像常数一样与求导交换）的性质（可以像常数一样与求导交换）xxx UAUA(U)f变系数，变系数， 不能与

19、导数交换不能与导数交换 A A,UxB AB实质：实质： 没有做到解耦；没有做到解耦； 只是把原变量重新组合，组合后波的传播方向的保证 f+ 向正向传播，f-向负向传播 缺点：缺点： 由于未解耦，各变量的误差会相互传递由于未解耦，各变量的误差会相互传递 26Copyright by Li Mingjun一般情况下一般情况下 9/2/202126 概念澄清：概念澄清： 流通流通矢量分裂本身不带来耗散矢量分裂本身不带来耗散， 但其会但其会影响到差分的耗散影响到差分的耗散;举例：举例：0)( xUftU, fff分裂过程：),(21),(21UffUff fffxxx UaftUxxx00 （耗散）

20、如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。ffffffxxxxx0000)( 分裂差分格式耗散 fff精确满足，不引入误差！27Copyright by Li Mingjun 9/2/202127=+向上平移向下平移分裂后的流场越偏离原先流场，则总体耗散越大分裂后的流场越偏离原先流场，则总体耗散越大如使用高精度格式（低耗散），则对分裂形式不敏如使用高精度格式（低耗散），则对分裂形式不敏感感 （可使用（可使用逐点分裂逐点分裂）28Copyright by Li Mingjun如使用低阶精度差分格式，如使用低阶精度差分格式， 则对分裂形式敏感则对分裂形式敏感 （推荐使用下面将介

21、绍的（推荐使用下面将介绍的特征分裂特征分裂）结论：结论： 9/2/202128(2) 严格特征分裂严格特征分裂 基架点上冻结系数常系数方程组：常系数方程组：00VV0USSU0UAU1 xvtvxtxtxtkkk 完全解耦29Copyright by Li Mingjun优点：优点： 严格保证（局部）特征方向，数值解质量好；缺点：缺点： 大量矩阵运算，计算量大。 9/2/2021290UAU0f(U)U xtxt(a) 变系数情况变系数情况: 局部冻结系数0UAU jjjxt j-2 j-1 j j+1 jjx UAjjx UA计算：在差分基架点上在差分基架点上Aj 不变，不变， 可按常矩阵处

22、理可按常矩阵处理jj1jSSA j)VV(SVS)US(SUSSUA1j1j1jj1jjjjjjjjjxxxxxx 局部冻结系数局部冻结系数分别采用后差和前差,USVj 30,USVj USVjCopyright by Li Mingjun 9/2/202130通常写成守恒型差分，计算x)(Ufxxjjj / )ff ()U( f2/12/1 j-2 j-1 j j+1 在基架点上系数 不变jjx UAjA)VV(f2/12/112/12/1 jjjjSUSV2/12/1 jj31,2kkk 2* kk均可均可.Copyright by Li Mingjun(b) 通量分裂表示通量分裂表示局部

23、特征值采用局部特征值采用 或者或者 9/2/202131(c) 具体步骤具体步骤;/ )(2/12/1xvvvjjj （1）1） 计算出计算出2/12/112/1, jjjSS各变量在j+1/2的值（例如 ）可使用j, j+1 点值的算术平均 （如 ）或Roe平均 ;2/1 ju2/ )(12/1 jjjuuu32 2/1j2/1 j2kkk 2* kkCopyright by Li Mingjun假设初值 U， 且针对模型方程（线性单波方程）已构造出差分格式计算；方法很多，例如前面介绍的或由./ )(2/12/1xvvvjjj 0UAU jjjxt 9/2/2021322） 在网格基上计算在网格基上计算)1, 1.(USV2/12/1 jjjkegkjjk j-2 j-1 j j+1 计算fj+1/2用到的点注意，注意，在该网格基上（例如k=j-1, j, j+1） 保持不变2/12/1,S jj3) 利用已构造好的差分格式，计算通量利用已构造好的差分格式，计算通量 2/12/1V,Vjj)VV(Sf2/12/112/12/1 jjj