计算流体力学(中科院力学所)-第7讲-差分方法3

计算流体力学讲义第七讲差分方法（3）李新亮[e-mail]；力学所主楼219；82543801知识回顾：守恒型方法通量分裂——逐点分裂与特征投影分裂知识点：精确特征方向的守恒型格式——Roe格式数值解的震荡机制及对策——TVD格式1讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”CopyrightbyLiXinliangCopyrightbyLiXinliang2知识回顾1）守恒型方程与守恒型格式守恒型方程：散度型0)(xuftu0)(UFUt守恒型格式：差量型1jjjxf习惯写为)(12/12/1jjjffxxf仅为记号，与j+1/2点上的值无关！守恒型方程+守恒型格式=解守恒“解守恒”：数值解的积分误差为0（如果边界准确）保证总量（总质量、总动量、总能量）严格守恒（无误差）0)(12/12/1jjjffxtu0)(10)(12/12/12/12/1Njjjjjjjffxutffxtu边界点早期的CFD:极为重视守恒性；目前CFD:仍很重视守恒性难点——复杂非线性系统的守恒性很难保证......)()()()(2/52/72/32/52/12/32/12/1ffffffffjjj中间项全部消去，只剩两端CopyrightbyLiXinliang2）通量分裂——便于使用迎风格式0xtf(U)U方法（1）：逐点分裂（Steger-Warming,VanLeer,L-F）)(0UfAAUxUtSUΛSUAffff1,原理：利用了性质xxUAAU)(使得UAf的Jocabian阵特征值纯正或纯负优点：无需矩阵运算，计算量小，使用方便不足：仅重新组合，没有做到真正解耦。原因：xxxUAUAf具体方法：Steger-WarmingL-FVanLeer2)(2/122kkk2*kk2kkkorCopyrightbyLiXinliang4方法2）特征（投影）分裂)()(jjjjjjjxxxxxVVSUΛSSUΛSSUA1j1jj1j00xtxtUAUf(U)U0jjjxtUAU在网格基上冻结系数…j-2j-1jj+1…在基架点上系数不变jjxUAjAUSΛVj)(12/12/1jjjxxVVV优点：特征分解，（局部）解耦——耗散小，数值振荡低缺点：大量矩阵运算，计算量大每计算一个点的导数，要进行m次矩阵运算（m为网格基点数）原理描述（非守恒，很少采用；实际上使用下一页的方法）守恒型方式计算x)(Ufxxjjj/)()(2/12/1ffUf…j-2j-1jj+1…在基架点上系数不变jjxUAjA)(2/12/112/12/1jjjjSVVfUSΛV2/1j具体步骤：假设已知U，且针对模型方程（线性单波方程）已构造出差分格式xvvvxvvvjjjjjj/)(;/)(2/12/12/12/10xvatv（1）1）计算出2/12/112/1,,jjjSS教材130页的公式(6.1.11-6.1.13),式中用到各变量在j+1/2的值（例如）可使用j,j+1点值的算术平均（如）或Roe平均；由计算；方法很多，例如前面介绍的或2/1ju2/)(12/1jjjuuu5CopyrightbyLiXinliang2kkk2*kk均可2/1j2/1j2kkk2*kk推荐使用Roe-平均！2）在网格基上计算)1,,1.(2/12/1jjjkegkjjkUSΛV…j-2j-1jj+1…计算fj+1/2用到的点注意，在该网格基上（例如k=j-1,j,j+1）保持不变2/12/1,jjΛSjjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuu1111112/12/)(2/)3(例如：1211121212/12/)(2/)3(jjjjjjjjjjjjjuuuuwhenuuuuuuwhenuuu3)利用已构造好的差分格式，计算通量4）得到总通量2/12/1,jjVV)(2/12/112/12/1jjjjVVSf5）计算差分（j点处）xxjjj/)()(2/12/1ffUf步骤的算法描述（注意：实际上是两重循环）doj=1,Ndok=j-1,j+1(网格基，可以是更多或更少点)enddoenddodoj=1,Nenddo),,(),,(112/1112/1jjjjjjjjffVV,VVVV,VV)(2/12/112/12/1jjjjSVVfxxjjj/)()(2/12/1ffUf需要多次矩阵运算，计算量大守恒性好，耗散小，数值解质量好6CopyrightbyLiXinliangkjjkUSΛV2/12/1通量分裂+迎风差分引入数值耗散分裂本身不带来耗散，但会放大（或减少）差分的耗散举例：0)(xUftUfff分裂过程),(21),(21UffUfffffxxxUftUxxx00耗散如果差分格式无耗散（例如都用中心差分），则通量分裂不带来耗散。ffffffxxxxx0000)(=+向上平移向下平移分裂后的流场越偏离原先流场，则总体耗散越大如使用低精差分度格式，则对分裂形式敏感（推荐使用特征分裂）如使用高精度格式（低耗散），则对分裂形式不敏感（可使用逐点分裂）7CopyrightbyLiXinliang概念澄清耗散放大系数CopyrightbyLiXinliang8§7.1Roe格式——守恒型格式的范例0)(xuftu))((,0)(ufuaxuuatu破坏守恒性后果很严重（?）为了便于使用迎风格式、特征分裂解耦，通常把守恒型方程改写为非守恒型守恒方程+守恒格式=解守恒方程不守恒，即使差分方法守恒，也无法做到解守恒)(12/12/1jjjuuxxu02/12/1xuuatujjjj02/12/1xuuautjjjjjj由于a非常数，无法消去中间项！......)()()()(2/532/732/322/522/112/312/12/1uauauauauauauuajjjj不再守恒思路：保证特征方向，找回守恒性守恒型方程优点：守恒性非守恒方程优点：清晰的特征方向兼顾守恒与非守恒方程CopyrightbyLiXinliang91.单方程的Roe格式0)(xuftu0)(假设,0)(ufuaxuuatu0/)(0/)(11axffaxffxfnjnjnjnjnj1阶迎风（直接从守恒方程出发）)(21][2112/12/11njnjnjnnjnjuuafffj)(12/12/1njnjnjffxxf0当)(0当,)()(),(111112/1jjnjnjnjnjnjnjnjjjnjuuuauuuuufufuuaa（1）（2）体现了特征方向只有这种表达式，才能保证（2）与（1）等价（3）2/)(12/1jjjaaa)2/)((12/1jjjuuaaor都无法保证（2）与（1）等价。简单的线性平均不行(非线性系统，中点的斜率不等于平均斜率）关键：构造2/1jaRoe格式“平均斜率”，不等于“斜率的平均值”，也不等于中点处的斜率njnjnjnjuuufuf11)()(平均斜率CopyrightbyLiXinliang102.方程组的情况（Roe格式的意义）0)(xtUfUUf(U)Af(U)AU,0xt)(12/12/1jjjxxfff)(ˆ12/1jjjU,Uff)U(U)U,(UA21)]f(U)[f(U21)U,(UfLRLRLRLR~ˆSΛS)U,(UA1LR~ΛSS)U,(UA1LR~需满足如下条件（Uniform特性）单方程的简单推广1）连续，且2）可通过相似变换对角化，即)(~UAU)(U,AΛSSA1~保证双曲性3）对于任意UR,UL有))(,()()()()(),(111111jjjjjjjjjjjjuuuuaufufuuufufuua)U)(UU,(UA)f(U)f(ULRLRLR~单方程的推广，含义为平均增长率标量方程向矩阵方程的简单推广，但构造很困难。)U,(UALR~)](21[),(~)],()([21),(~VUAVUAVAUAVUA均不满足Uniform特性)(ˆ)(ˆ1LRjjU,UfU,Uf经常记为)U,(UALR~)U,(UALR~)U,(UALR~)U,(UALR~平均斜率CopyrightbyLiXinliang113.矩阵的构造)U,(UALR~)U)(UU,(UA)f(U)f(ULRLRLR~关键：错误！/)(~)U(U)f(U)f(U)U,(UALRLRLR“向量除以向量”？直接求平均增长率：uf(u)uLuRuRoeRoe点的斜率为平均斜率（根据拉格朗日中值定理，[UL,UR]区间内肯定存在Roe点）思路1：在UL与UR之间寻找一个点URoe,该点处的增长率为平均增长率f(u)=u2u二次函数——Roe点与中点重合标量函数的启示：Roe点肯定存在（Langrage中值定理）二次函数的中点即为Roe点思路2：进行坐标变换，得到一个二次（齐）函数F(W)F(U(W))F(U)U(W)U引入如果是二次（齐）函数，则其中点即为Roe点重要启示F(W))/2W(WWRL更准确地讲，应当是要求为W的线性函数，即增长率为线性函数（中点处的增长率刚好为平均增长率）WF(W)CopyrightbyLiXinliang1222312121211U(W)UHuWpEH针对Euler方程的具体构造引入新变量：则：目的：使得F(w)是W二次齐函数（增长率为线性函数）0xtf(U)UTTpEupuuEu))(,,()(,),,(2UfUf(U)不是U的二次齐函数23223121211)(Wf二次齐函数！)W)(WWC(ffLRLR)/2W(WWLR中点处的斜率即为平均斜率。Roe点Roe点为：)WU(U23122312011210wf(W)C(W)增长率为线性函数！CopyrightbyLiXinliang13最终：)UA()U,(UALR~)2)(1(22uHc其中如下计算：U平均增长率（矩阵）)/()()/()(]2/)[(2RLRRLLRLRRLLRLHHHuuu含义：左、右两个状态点的某种平均（称为Roe平均，为密度加权平均）该状态点对应的增长率（矩阵）为平均增长率（矩阵）实际上是一种“