1、计算流体力学讲义计算流体力学讲义2020 第五讲第五讲 差分方法差分方法（3）李新亮李新亮知识点：知识点： 通量技术简介通量技术简介Steger-Warming, RoeSteger-Warming, Roe 常用的隐式处理方法常用的隐式处理方法LU-SGSLU-SGS 1Copyright by Li XinliangCopyright by Li Xinliang2知识回顾：知识回顾： 单调、保单调和单调、保单调和TVD概念：概念： 网格网格Reynolds数数 单调格式、保单调格式及单调格式、保单调格式及TVD格式格式 Harten定理：定理： 正系数原则正系数原则xxReReTVD保单
2、调保单调单调单调TVD格式格式= 1阶迎风阶迎风+ j j *（修正项）（修正项）二阶精度区二阶精度区TVD区区二阶精度二阶精度TVD区（二区（二者交集）者交集）2/ )(12/1jjjjuuruujjjjjjjjuuuurr11 ),(jjCopyright by Li Xinliang30 xuatu0a0a0 xtf(U)UStep 1 针对模型方程构造差分格式针对模型方程构造差分格式Step 2 将格式推广到将格式推广到Euler方程方程方法方法1： 流通矢量分裂流通矢量分裂(FVS) (U)f(U)ff(U)通量分裂技术：通量分裂技术： 模型方程模型方程 NS/ Euler 方程方程
3、1/21/2jjuuuxx1/2=.ju1/2=.ju格式1格式2（利用对称性， 将格式1 中的 j+k替换成j-k即可）fffxxxCopyright by Li Xinliang4方法方法2： 通量差分分裂通量差分分裂 （FDS）特点：特点： 对流通矢量对流通矢量 f 的导数进行分裂的导数进行分裂 Godnov, Roe, HLL,HLLC 方法方法3： AUSM 类方法类方法van Leer分裂法分裂法+压力项单独处理压力项单独处理u 利用利用 Riemann解解1/21/2jjjffxxf0 xuatu利用差分表达式，计算利用差分表达式，计算1/21/2,LRjjUU 求解求解Riem
4、ann问题，获得通量问题，获得通量1/2jf1. Jacobian 系数矩阵及其性质系数矩阵及其性质xxUAf(U)Uf(u)Af(U)Uf)(f(U)UUU)f()()(1AUUUf(U)f(U)(重要性质uucucuuu222322311221)3(23010Avisf(U)Uxt,),(TEuU5Copyright by Li Xinliang21222231332322113(1)2()12uufufupuuu Epfu uuuuf(U)f(U)AU2222320103(3)1223()2112uupuuuuccuuuEp uE容易验证：容易验证：验证方法2： 利用齐函数性质 5.1
5、流通矢量分裂：流通矢量分裂： Steger-Warming 、L-F2. 流通流通矢量分裂矢量分裂(FVS)kdiagkdiag2,2kk=+优点：耗散小缺点：导数间断fffA: Steger-Warming 分裂分裂6Copyright by Li Xinliangf(U)AUSUSUAf1向量分裂向量分裂 特征值标量分裂特征值标量分裂任意函数分解为：任意函数分解为：非负函数非负函数+非正函数非正函数2,2kk=+优点：耗散小缺点：导数间断SUSUAf12)(2/122kkk) 1(2)(3(232cw特点：特点： 不必进行矩阵运算，计算量小不必进行矩阵运算，计算量小Steger-Warmi
6、ng 分裂分裂wcucuucucuu232221321321)(2)(2) 1()()() 1(2) 1(22)(f)(ff7Copyright by Li Xinliang2kkk改进版 Steger-Warming 具体步骤具体步骤 （以一维为例）（以一维为例）已知已知TEuU),(1） 计算计算2） 计算计算3） 计算计算4） 带入（带入（1）式得到）式得到5） 利用不同的迎风格式，分别计算利用不同的迎风格式，分别计算 pu,cucuu321,2)(2/122kkk)3 , 2 , 1(,kkkwcucuucucuu232221321321)(2)(2) 1()()() 1(2) 1(2
7、2)(f) 1(2)(3(232cw)(),(ffff（1）ff ,fffxxff,xxff0a0 xuatu0a（后差，前差）（后差，前差） 6） xxxfff计算计算7） 时间推进时间推进0 xtfU8Copyright by Li Xinliang二维问题的二维问题的steger-Warming 分裂分裂021yxtffUTpEuuvpuu)(,21fTpEvpvuvv)(,22fTEvu,U21fff令：则：WvuvuVvvvuuu2222)(2222421213202413024130430f) 1(2)(3(243cW10) 1(211kcuu12kcuu/1k/2k21kcvv2
8、2kcvv22具体使用步骤，具体使用步骤， 以计算以计算 为例为例x1f1) 令令 2) 计算特征值计算特征值3) 分裂特征值，计算分裂特征值，计算4) 带入左式，计算正、负流通矢带入左式，计算正、负流通矢量量5) 计算计算0, 1cucuu4321,2)(2/122kkk)4 , 3 , 2 , 1(,kkk)(1ffxxx111fffy2f计算计算 设置设置 ，并注意，并注意1, 0cvcvv4321,对于曲线坐标系对于曲线坐标系021ffUt21111fffyxJJ仅需令yxJJ11, 三维问题同样处理三维问题同样处理 二维、三维具体二维、三维具体 公式见傅德薰等公式见傅德薰等计算空气动
9、力学计算空气动力学 4.7节节 （158-162）书中公式有一定的排版错误，使用前务必书中公式有一定的排版错误，使用前务必重新仔细重新仔细推导！推导！9Copyright by Li XinliangB: Lax-Friedrichs (L-F)分裂分裂AUf(U) 2/ )(, 2/ )(*IAAIAA特点： A正特征值 负特征值A2/ )(*UffUA=+缺点：耗散偏大局部局部L-F分裂分裂，每个点上计算，每个点上计算 全局全局L-F分裂分裂，全局（一维）上计算，全局（一维）上计算 cu *足够大cu *数学性质（光滑性）数学性质（光滑性）最好，但耗散偏大最好，但耗散偏大)(max*cux
10、常数与迎风格式结合，等价于人工粘性与迎风格式结合，等价于人工粘性*0*11()()()()221() () 22xxxxxxxxxxxxxxxffffUfUfUfUfUfU例如，可取例如，可取0 xtf(U)U2*2xtxxUf(U)U10Copyright by Li Xinliang2,2kkkkkkkkk方式很多=+S-W:L-F:=+2,2*kkkk0, 0kkVan Leer:=+221112) 1() 1(22) 1(cufcufff2121Macfff0f1Ma1MaffcuMa/0f1Ma11Copyright by Li Xinliang 分裂后分裂后 失去了失去了A的性质（
11、可以的性质（可以像常数一样与求导交换）像常数一样与求导交换）FVS分裂：分裂： 优点：优点： 无需矩阵运算，计算量小无需矩阵运算，计算量小 缺点：缺点： 分裂后改变了特征方向，分裂后改变了特征方向， 耗散大耗散大xxx(U)f(U)ff(U)SUSUAf1利用了性质)(AUUAf(U)xxx一般情况下：一般情况下：xxxUAUAf变系数，变系数， 不能与导数交换不能与导数交换AAxBxUf AB实质：实质： 没有做到解耦；没有做到解耦； 只是把原变量重新组合，组合后波的传播方向的保证只是把原变量重新组合，组合后波的传播方向的保证 f+ 向正向正向传播，向传播，f-向负向传播向负向传播 缺点：缺
12、点： 由于未解耦，各变量的误差会相互传递由于未解耦，各变量的误差会相互传递 12Copyright by Li Xinliang概念澄清：概念澄清： 流通流通矢量分裂本身不带来耗散矢量分裂本身不带来耗散， 但其会但其会影响到差分的耗散影响到差分的耗散;举例：0)(xUftUfff分裂过程),(21),(21UffUfffffxxxUaftUxxx00耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。ffffffxxxxx0000)(=+向上平移向下平移分裂差分格式耗散分裂后的流场越偏离原先流场，则分裂后的流场越偏离原先流场，则总体耗散越大总体耗散越大fff精确满足，不引入误差！
13、如使用低精差分度格式，如使用低精差分度格式， 则对分裂形式敏感则对分裂形式敏感 （推荐使用特征分裂）（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感如使用高精度格式（低耗散），则对分裂形式不敏感 （可使用逐点分裂）（可使用逐点分裂）13Copyright by Li XinliangCopyright by Li Xinliang14 FVS差分方法一般流程差分方法一般流程Step 1. 构造差分格式构造差分格式0uuatx1/21/2iiuuuxx1/2111/21/211(.,.)0=(.,.)0iiiiiiiiiuuu uifauuuu uifa差分格式Step 2.
14、推广到方程组推广到方程组0txUf(U)+-=+xxxfff+1/21/2=jjxxfff1/21/211=,.)jjjjjfffff（.,具体步骤具体步骤 （逐点分裂）：（逐点分裂）：,jjff通量分裂通量分裂jf差分格式差分格式+1/2+1/2,jjff+1/21/2=jjxxfff+-=+xxxfffj=1,2,3N3. 特征重构方法特征重构方法常系数常系数方程组：方程组：00 xtxtUAUf(U)U0000 xvtvxtxtxtkkkVVUSSUUAU1完全解耦变系数情况变系数情况 局部冻结系数局部冻结系数0jjjxtUAU j-2 j-1 j j+1 在基架点上系数 不变jjxUA
15、jjxUAjA计算：在差分基架点上在差分基架点上Aj 不变，不变， 可按常矩阵处理可按常矩阵处理jj1jSSAj)()(jjjjjjjjxxxxxxVVSVSUSSUSSUA1j1j1jj1j局部冻结系数分别采用后差和前差USVj优点： 严格保证（局部）特征方向，数值解质量好；缺点： 大量矩阵运算，计算量大。15Copyright by Li XinliangUSVjUSVj16具体步骤具体步骤 （特征分裂）：（特征分裂）：通量分裂通量分裂kf j-2 j-1 j j+1 j+1/2,kkffk=1,2.N1/2jU2）针对）针对j点：计算特征矩阵点：计算特征矩阵11/21/21/2=jjjA
16、SS1） 逐点分裂逐点分裂3）（将网格基架点）变换到特征空间）（将网格基架点）变换到特征空间1/2=kjkfSf1/2=kjkfSfk=, j-2, j-1, j,网格基架点4）在特征空间计算）在特征空间计算j+1/2值值1/21/211=,.)jjjjjfffff（.,差分格式， 例如：5）变换回物理空间）变换回物理空间11/21/21/21/21/21/2()()jjjjjjfffSff1/21=)/ 2jjjUUU（简单平均即可（Roe平均效果更好些）1/211=+1/2minmod(,)jjjjjjffffff1/21/211=,.)jjjjjfffff（.,6） 时间推进时间推进1/
17、21/20jjjtxUff1/2jS11/2jSCopyright by Li Xinliang17 5.2 Roe格式格式( )0UF Utx通量分裂技术：通量分裂技术：流通矢量分裂（流通矢量分裂（FVS）：）：( )( )( )F UFUFUFFFxxx0uuatx0a 0a 1/21/2jjFFx通量差分分裂通量差分分裂 （FDS）：）：Step 1: 运用差分格式，计算运用差分格式，计算1/21/2,LRjjUUStep 2: 运用运用Riemann解，解， 计算计算1/21/21/2(,)LRjjjFF UUStep 3:1/21/2jjFFFxxjj+11/21/2,LRjjUU精
18、确精确Riemann解计算量大，且效果未必比近似解计算量大，且效果未必比近似Riemann解好解好Roe 格式是目前应用最为广泛的近似格式是目前应用最为广泛的近似Riemann解解耗散低、分辨率高稳定性好、低振荡Copyright by Li Xinliang181. 单方程的单方程的Roe格式格式0)(xuftu0)(xuuatu11/21/2111()22jjjjjjfffauu)(12/12/1jjjffxxf线性化，用线性化，用平均变化率平均变化率代替代替(j,j+1)之间的变化率之间的变化率a(u)“平均斜率平均斜率”，不等于，不等于“斜率的平斜率的平均值均值”，也不等于中点处的斜率
19、，也不等于中点处的斜率非线性情况非线性情况jjjjjuuufufa112/1)()(1/21/2+11/200jjjjjfaffa2/ 1ja)(21)(12/ 12/ 1jjjjaauaa根据Langrage中值定理，uL,uR之间必有一点uRoe, 该点处的斜率为平均斜率;二次函数f(u)=u2中点处的斜率=平均斜率Copyright by Li Xinliang192. 方程组的情况方程组的情况 0)(xtUfUUf(U)AUAU, 0)(xUt)()(1LRjjU,UfU,Uf经常记为)U,(UALR平均斜率平均斜率线性化，以平均增长率代替瞬时增长率)(12/12/1jjjffxxf0
20、 xtUAUj,j+1区间内)U)(UU,(UA)f(U)f(ULRLRLR连续连续，且，且可可通过相似变换通过相似变换对角化对角化)(UAU)(U,A)U,(UALR)U,(UALR 应当具有的性质应当具有的性质)U,(UALR)U(U)U,(UA21)f(U)f(U21fLRLRLR2/1j常系数方程的Riemann解SS)U,(UA1LRRULUCopyright by Li Xinliang20平均斜率平均斜率0 xtUAU2/12/1jRjLxxUxxUU1SSAVSU0 xVtVkkk2/12/1jRkjLkkxxVxxVVxkVj+1/2LkVRkV)(2sgn2002/1LkR
21、kkRkLkkRkkLkjkVVVVVVV常系数单波方程的Riemann解)(sgn(21)(212/12/12/12/12/1LjRjRjLjjVVVVVVSU1)()sgn(21)(2112/1LRRLjUUSSUUU2/12/12/1jjjUAf)sgn()U(U21)f(U)f(U21fLjRjRjLjjSS2/12/112/12/12/1Roe 格式：格式：微分型近似微分型近似Riemann解解Copyright by Li Xinliang213. 矩阵矩阵 的构造的构造)U,(UALR)U)(UU,(UA)f(U)f(ULRLRLR关键：关键：错误！/ )()U(U)f(U)f(
22、U)U,(UALRLRLR“向量除以向量向量除以向量” ？直接求平均增长率：uf(u)uLuRuRoeRoe点的斜率为平均斜率点的斜率为平均斜率（根据拉格朗日中值定理，（根据拉格朗日中值定理，UL,UR区区间内肯定存在间内肯定存在Roe点）点）思路思路1： 在在UL与与UR之间寻找一个点之间寻找一个点URoe, 该点该点处的增长率为平均增长率处的增长率为平均增长率f(u)=u2u二次函数二次函数 Roe点与点与中点重合中点重合标量函数的启示：标量函数的启示： Roe点肯定存在（点肯定存在（Langrage 中值定理）中值定理） 二次函数的中点即为二次函数的中点即为Roe点点思路思路2： 进行坐
23、标变换，得到一进行坐标变换，得到一个二次（齐）函数个二次（齐）函数F(W)F(U(W)F(U)U(W)U引入如果如果 是二次（齐）函数，则其是二次（齐）函数，则其中点中点 即为即为Roe点点重要启示F(W)/2W(WWRL更准确地讲，应当是要求更准确地讲，应当是要求 为为W的线性函数，的线性函数， 即增长率为线性函数即增长率为线性函数 （中点处的增长率刚好为平均增长率）（中点处的增长率刚好为平均增长率）WF(W)Copyright by Li Xinliang22Huwww12/1321WpEH针对针对Euler方程的具体构造方程的具体构造引入新变量：则：目的：目的： 使得使得F(w)是是W二
24、次齐函数二次齐函数 （增长（增长率为线性函数）率为线性函数）0 xtf(U)UTTpEupuuEu)(,()(,),(2UfUf(U)不是U的二次齐函数1221323211()2w ww www wf W二次齐函数！二次齐函数！)W)(WWC(ffLRLR)/2W(WWLR中点处的斜率即为平均斜率。中点处的斜率即为平均斜率。Roe点Roe点为：)WU(U 23122312011210wwwwwwwwwf(W)C(W)增长率为线增长率为线性函数！性函数！Copyright by Li Xinliang23最终：)UA()U,(UALR)2)(1(22uHc其中其中 如下计算：如下计算：U平均增长
25、率（矩阵）平均增长率（矩阵）)/()()/()(2/ )(2RLRRLLRLRRLLRLHHHuuu含义：含义：左、右两个状态点的某种平均左、右两个状态点的某种平均 （称为（称为Roe平均，为密度加权平均）平均，为密度加权平均） 该状态点对应的增长率（矩阵）为平均增长率（矩阵）该状态点对应的增长率（矩阵）为平均增长率（矩阵） 实际上是一种实际上是一种“等效平均等效平均”。 效果优于简单的算数（或几何）平均。效果优于简单的算数（或几何）平均。 )/()()/()(RLRRLLRLRRLLwwwvvv三维情况下，还有其他量（如压力、温度、音速等）用这三个量计算其他量（如压力、温度、音速等）用这三个
26、量计算pEH)21(12uHp（5）() / 2() / 2() / 2LRLLRRLLRRuuuHHH简单易记：Copyright by Li Xinliang24 Roe 格式的计算步骤格式的计算步骤 （半离散）（半离散）0)(xtUfU已知已知n时刻所有网格点上的物理量，对于时刻所有网格点上的物理量，对于j点：点： 1） 利用差分格式计算利用差分格式计算UR,UL 2） 采用采用Roe平均公式（平均公式（5）计算）计算Roe平均值平均值 3） 将将Jacobian矩阵矩阵 进行特征分解进行特征分解: 计算计算 4) 计算计算 5）计算）计算 6） 计算空间导数计算空间导数 7）时间推进，
27、计算下一时间步的值。）时间推进，计算下一时间步的值。 j-1 j j+1U)UA(SS)UA(1与前文（第与前文（第3，4讲）的形式相同，讲）的形式相同，仅需把式中的密度、压力、速度等换仅需把式中的密度、压力、速度等换成经过成经过Roe平均的密度、压力、速度平均的密度、压力、速度即可即可SS1,SS)U,(UA1LR)U(U)U,(UA21)f(U)f(U21)U,(UffLRLRLRLR2/1j)(1)(1/2j1/2jffUfxxj其中：),(321diag推荐：推荐： 特征投影特征投影Copyright by Li Xinliang25),(321diag可能出现导数不连续，可能出现导数
28、不连续， 可能引起数值振荡可能引起数值振荡实际使用时实际使用时 可用如下函数代替可用如下函数代替 所谓所谓“熵修正熵修正”当当2/ )()(22fk)(f实际上是在特征值实际上是在特征值0点周围增加了耗散点周围增加了耗散Roe 格式的优点：格式的优点： 1） 保持守恒性的同时，严格保证了特征方向保持守恒性的同时，严格保证了特征方向 2） 便于推广到高精度格式便于推广到高精度格式 特征投影分裂中使用特征投影分裂中使用Roe平均即可平均即可 （见（见本本PPT 第第5页）。推广到高阶后，虽不再保证严格的特征方向，页）。推广到高阶后，虽不再保证严格的特征方向， 但仍优于采但仍优于采用算数平均方法。用
29、算数平均方法。 Roe 格式的不足：格式的不足： 原始的原始的Roe格式本身格式本身精度只有一精度只有一阶阶 （可通过插值推广到高阶）；（可通过插值推广到高阶）； 推广到高阶后，特征方向无法严格保证推广到高阶后，特征方向无法严格保证 ； 推广到二维或三维后，特征方向无法严格保证，出现振荡。推广到二维或三维后，特征方向无法严格保证，出现振荡。Copyright by Li Xinliang26 关于关于 f(U) 与与 f(W)Euuuu311U)(2311pEupuuffff深入讨论深入讨论Huwww12/1321W新变量新变量213213231222) 1(21) 1()3(21uuuuuu
30、uuuf(U)23223121211)(wwwwwwwWf虽然是虽然是“一次一次齐函数齐函数”但有但有变量在分母上变量在分母上干净的二次干净的二次齐函数齐函数23122312011210wwwwwwwwwf(W)C(W)自变量自变量W的线性函数！的线性函数！实质区别：实质区别： 是自变量是自变量 的线性函数，的线性函数， 而而 是自变量是自变量 的非线性函的非线性函数！数！ Uf(U)UUf(W)W)()(1(23)(1(1)(3()(3(210101221213212213212212uuuuuuuuuuuuuuuUf(U)A(U)自变量自变量U的非线性函数的非线性函数使用使用U做自变量的优
31、点：做自变量的优点： 物理意义鲜明物理意义鲜明 （质量密度、动量密度和能量密度），守恒性好（质量密度、动量密度和能量密度），守恒性好使用使用W做自变量的优点：做自变量的优点： Jocabian矩阵为线性矩阵矩阵为线性矩阵思考：思考： 如果在如果在CFD计算中，使用计算中，使用W替换替换U做自变量会怎样？做自变量会怎样？4 阶3 阶阶 （TVD型）型）2阶Copyright by Li Xinliang27 5.3 时间推进方法时间推进方法 1. 显格式显格式z(U)gy(U)gx(U)gz(U)fy(U)fx(U)ftU321321Q(U)tU推荐方法：推荐方法： Runge-Kutta法法)
32、(3/23/23/1)(4/14/14/3)()2()2(1)1()1()2()1(UtQUUUUtQUUUUtQUUnnnnn)(2/12/12/1)()1()1(1)1(UtQUUUUtQUUnnnn)(6/1)2(3/1)()(2/1)(2/1)3()3()2()1(1)2()3()1()2()1(UtQUUUUUUtQUUUtQUUUtQUUnnnnnn更高阶 Copyright by Li Xinliang282. 隐格式算法简介隐格式算法简介 （以二维（以二维Euler方程为例）方程为例） 1） 原理介绍原理介绍0yxt(U)f(U)fU21011yxtnn)(Uf)(UfUU21
33、n1n（1）时间离散后时间离散后方案方案 A： 直接将（直接将（1）进行空间离散，得到）进行空间离散，得到 Un+1 的代数方程组的代数方程组困难：困难： 大型大型非线性非线性方程组，求解困难方程组，求解困难方案方案B: 设计一种迭代算法设计一种迭代算法令n1nnUUQ01211yxtnnffQn)()(21212111yxtyxtnnnnnnffffffQn两端同时添加显式项两端同时添加显式项右端项，已知右端项，已知21yxttnnffRHS是已知项，可采用某种差分方法是已知项，可采用某种差分方法显式显式计算得到计算得到（对算法无限制）（对算法无限制）Copyright by Li Xinl
34、iang29tRHSyxtnnnn)()(212111ffffQnUfUfA1)()(,)()(2UUBUU)(),(22122111tOtOnnnnnnUBffUAffRHSBQAQQntyxtnn)B(UB),(UAAnn(2)已知，（已知，（2）为线性方程）为线性方程。（1） (2) 是一个线性化过程是一个线性化过程 含义：含义： 先用显格式计算，再先用显格式计算，再用隐格式计算修正量用隐格式计算修正量线性方程，离散求解线性方程，离散求解离散后为大型带状方程组，求解计算量大离散后为大型带状方程组，求解计算量大LU-SGS原理：原理： 将矩阵分解为上、下三角阵，避免矩阵求逆运算将矩阵分解为
35、上、下三角阵，避免矩阵求逆运算隐式问题显式化隐式问题显式化UQ 未知量：nn1nnUUUQ显式部分隐式修正若隐式修正为若隐式修正为0，则为显格式，则为显格式Copyright by Li Xinliang301,11,1ijijijijiji jijijiji jijabcdeQQQQQRStep 1: 求解1,1ijijijijiji jijadeQQQR1,1()/ijijijijiji jijdeaQRQQ111111/ aQRStep 2: 求解1,1ijijijijiji jijijabcaQQQQ1,1()/ijijijijijiji jijabcaQQQQ/mnmnmnaQRLU
36、-SGSi,ji-1,ji+1,ji,j+1i,ji-1,ji+1,ji,j+1i,j+1i,j+1ULDLUD1)()(1UDDLD2） LU-SGS方法方法原理简介：忽略右、上未知量忽略左、下未知量i,ji-1,ji+1,ji,j+1i,j+1原则上：原则上： （1） （2） 两步反复迭代，直至收敛两步反复迭代，直至收敛Copyright by Li Xinliang31LU-SGS 方法求解方法求解Euler方程方程a) 将矩阵将矩阵A,B分裂分裂BBB,AAA为了简化计算，通常采用为了简化计算，通常采用L-F分裂分裂RHSBQAQQntyxtnnRHSQBQBQAQAQntyyxxtn
37、nnn矩阵分裂矩阵分裂1阶迎风格式离散阶迎风格式离散RHSQQQQQAQAQAQAQtBBBBytxtjijijijijijijijijijijijijijijijiij,1,1,1,1, 1, 1, 1, 1,)(max(, 2/ )(*AkAAIAA)(max(, 2/ )(*BBBkBBIRHSQQQAQAAAQtBBytxtBBytxtIjijijijijijijijijijijijiij)()(1,1,1,1, 1, 1, 1, 1,IBBIBAjijijiji*,*,AA整理整理迭代收敛迭代收敛后，不影响精度后，不影响精度Copyright by Li Xinliang32RHSQ
38、QAQQAQtBytxtBytxtytxtjijijijijijijijiBAij1 1,1, 1, 11,1, 1, 1*RHS)LU(DtUULDLUD1)()(1UDDLDUDU,LDLRHSUQLD1t)(RHSDLUQ11t1 1,1, 1, 1*jijijijiijBABytxtytxtQQAQLQ其中：其中：ijBAytxtQDQ1 *1 1,1, 1, 1*jijijijiBAijBytxtytxtQQAQUQ特点：特点： 严格对角占优，收敛性好严格对角占优，收敛性好 稳定性好稳定性好可采用局部时间步长等加速收敛措施可采用局部时间步长等加速收敛措施近似近似LU分解分解 do j
39、=1,ny do i=1,nxEnddoEnddoCopyright by Li Xinliang33具体求解方法具体求解方法0yxt(U)f(U)fU21（1） 计算右端项计算右端项yxnn21ffRHS显式，可采用前面构造的各种方法计算显式，可采用前面构造的各种方法计算 （2） 计算计算)(max(, 2/ )(*AkAAIAA)(max(, 2/ )(*BkBBIBB（3） 计算计算 自变量采用自变量采用n时刻的值时刻的值RHSDLQ1)1(由于由于L矩阵的下三角特性，可采用推进方法矩阵的下三角特性，可采用推进方法显式显式计算计算1,1, 1, 1)1(jijijijiijijytxtR
40、HSQBQAQ(4)计算计算(1)1QUQt1/*)1(1,)1(, 11, 1BAjijiijytxtytxttjijiQBQAQ显式计算显式计算（5）QUUn1n 时间推进时间推进1阶精阶精度，空间精度由度，空间精度由RHS的算法决定。的算法决定。 可用于定常及非可用于定常及非定常问题；定常问题； 非定常问题通常非定常问题通常采用采用内迭代内迭代 （保（保证精度）证精度）向上扫描向上扫描向下扫描向下扫描LU分解的对称分解的对称Gauss-Seidel迭代算法迭代算法Copyright by Li Xinliang341） 可简化计算（不出现矩阵运算）：可简化计算（不出现矩阵运算）：*1,1
41、,1,1,*1,1,1,1,1,()ijijijijijijijijij AQAIUAUUFQ1,1, 1, 1)1(jijijijiijijytxtRHSQBQAQ流通量，利用特性流通量，利用特性,ddFAUFAUdtdt已知量i,j-1i,j-1i,j-1F(U+Q)-F(U)UQ 2） 一个时间步仅进行一次一个时间步仅进行一次 “上扫上扫”和一次和一次“下扫下扫”过程过程 （不迭代）（不迭代） 注解：注解：3）对于定常问题，）对于定常问题， 精度由右端项决定，精度由右端项决定，LU-SGS方法并不降低精度方法并不降低精度Copyright by Li Xinliang353. 非定常问题
42、的双时间步法非定常问题的双时间步法目的：定常问题的加速收敛技术目的：定常问题的加速收敛技术 非定常问题非定常问题()WR Wt1） 构造隐格式构造隐格式例：例：11134()2nnnnWWWR Wt2阶精度隐格式2） 构造发展方程构造发展方程()0()WR WR Wt定常问题，非定常处理，推进到收敛*1*34()2nnWWWWR Wtt3）对）对t* （伪时间）（伪时间） 进行时间推进，直到收敛进行时间推进，直到收敛*1*340()2nnWWWR Wt推进到收敛推进到收敛可使用定常问题的加速收敛手段：可使用定常问题的加速收敛手段： 局部时间步局部时间步长法、隐格式、多重网格法长法、隐格式、多重
43、网格法1*nWW针对如下针对如下Sod 激波管问题激波管问题0)()(0)()(0)(2xpuEutExputuxut5 . 01 . 0 ,125. 0 , 05 . 01 , 1 , 0),(:0 xxput 1 , 0 x 用用 激波捕捉格式激波捕捉格式计算计算其数值解，画出其数值解，画出t=0.14时刻密度、速度及压力的分布；并时刻密度、速度及压力的分布；并与精确解进行比较与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。（要求数值解与精确解画在同一张图上，便于比较）。 要求要求： 空间网格数空间网格数100， 时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长时间步长自选；自选； 激波捕捉格式可自选（如激波捕捉格式可自选（如TVD, NND, GVC, MUSCL等）；等）； 36Copyright by Li Xinliang作业作业 5.1 初值提示：可采用基于流通矢量分裂方法， 具体步骤见 p8;也可以采用 流通矢量特征分裂的方法， 具体步骤见 p16 （推荐）；也可以采用 Roe 格式计算, 具体步骤见 p24 (重构时，推荐使用特征变量 ）
