1、Copyright by Li Xinliang1知识回顾知识回顾 有限体积法有限体积法0y(U)fx(U)ftU210F(U)tU在以某节点为中心的在以某节点为中心的控制体控制体上积分上积分i,jk非结构网格的控制体非结构网格的控制体i+1,ji-1,ji,j+1i,j-1k3k1k2k4k5结构网格的控制体结构网格的控制体0dF(U)tU01dsnFtUdUU1dsdnFF(U)xyn01mmijHtU体积平均)()(1ymxmmmmnUnUss2ffnFH控制体边界垂控制体边界垂直于节点连线直于节点连线（也可选其他（也可选其他方式）方式）垂直平分线垂直平分线n1) 建立控制体建立控制体
2、mx2) 在控制体上积分在控制体上积分离散方程离散方程重构：重构： 由节点上平均值由节点上平均值 给出函数分布，最终给出函数分布，最终给出通量给出通量ijU表示第m个界面上的值1m2m3m4m1. 有限体积法的离散过程有限体积法的离散过程重构重构1） 重构重构 两种不同的重构方案，向左偏及向右偏。两种不同的重构方案，向左偏及向右偏。 给出两种结果：给出两种结果： 及及Copyright by Li Xinliang2LjiU,2/1RjiU,2/1例如：例如： 0阶重构：阶重构：jiLjiUU,2/1jiRjiUU, 1,2/1 线性重构线性重构: )(21, 1, 2/1jijijiLjiU
3、UUU)(21, 1, 2, 1, 2/1jijijiRjiUUUU)()(2)(, 1, 1, 1, 2/1jijijijijjijiLjiUUxxxxUU2. 迎风型有限体积法迎风型有限体积法注：注： 有激波的情况，需使用限制器有激波的情况，需使用限制器 （即激波捕捉格式）（即激波捕捉格式）j-1jj+1/2如不使用限制器，很可能如不使用限制器，很可能“寸寸步难行步难行” （第（第1步就发散）步就发散）例：例： Sod 问题问题 35. 02/ )(, 1 . 0, 112/11jjjLjjjpppppp11jp35. 02/1Ljp1 . 0jp压力出现负值，计算发散压力出现负值，计算发
4、散),mod(min21, 1, 1, 2/1jijijijijiLjiUUUUUU),mod(min21, 1, 2, 1, 1, 2/1jijijijijiRjiUUUUUU例如：例如：NND格式格式在在 二阶中心二阶中心 / 二阶迎风二阶迎风/1阶迎风阶迎风 三个方案中选取最合适的。三个方案中选取最合适的。1阶迎风阶迎风修正项修正项线性插值， 得到二阶迎风格式两个修正方案（二阶迎风及二阶中心），给出趋势相两个修正方案（二阶迎风及二阶中心），给出趋势相同时，选修正量最小的，趋势相反时，不修正同时，选修正量最小的，趋势相反时，不修正Copyright by Li Xinliang3i,ji+
5、1,ji-1,ji,j+1i,j-1n左重构左重构右重构LmU2） 由左右重构得到的自变量：由左右重构得到的自变量： 和和 给给出通量出通量 方案方案A: FVS 方案方案B: 解解Riemann问题问题 （常用）（常用）LjU2/1RjU2/11111,pvu2222,pvuxy看似二维看似二维Riemann问题，其实是一维问题，其实是一维的，坐标旋转一下的，坐标旋转一下就行了就行了RmU)(RmLmmU(Uf(UffA. 精确精确Riemann解解B. Roe-近似近似Riemann解解 （求解常系数线性化的（求解常系数线性化的Euler方程）方程）C. HLL型近似解型近似解 双激波近似
6、，积分平均解双激波近似，积分平均解D. HLLC型近似解型近似解 三波近似三波近似3. 粘性通量的计算粘性通量的计算 中心型中心型Copyright by Li Xinliang4Part 1 代数方程组的求解代数方程组的求解微分方程（组）微分方程（组）代数方程组代数方程组数值解（离散解）数值解（离散解） 差分差分有限体积有限体积nnnnnnnnaaaaaaaaaaaaaaaaA321333323122322211131211大部分计算量大部分计算量11.1.1 Gauss 消去法消去法bAx nnnnnaaaaaaaaaaA0000000333223221131211nnnnnnnaaaaa
7、aaaaaaaaA3233332223221131211000消元nnnnabx/iiknikikiiaxabx/ 1为了计算稳定，通常使用主元消去法为了计算稳定，通常使用主元消去法列主元消去法；列主元消去法； 全主元消去法全主元消去法计算量：计算量： 乘法：乘法： 加法：加法：11.1 代数方程组求解的直接法代数方程组求解的直接法3/3/23nnn6/52/3/23nnn) 3/(3nO优点：优点： 简单精确，缺点：计算量大简单精确，缺点：计算量大上三角矩阵Copyright by Li Xinliang511.1.2 LU分解法分解法bAx LUA nnnnnnnnnnnnnnnnuuuu
8、uuuuuullllllaaaaaaaaaaaaaaaa3332232211312113213231213213333231223222111312111111nkkjikijula1jjnkjkjkjauuula111111jnkjjjjkjkjulauuulula121122212122001111111111111/aaualululaiiinkikiki2212112222212122/ )(uulalulululainkiiiikikimjkmkmkjkjulau11Step 1Step 2Step kkkmkimkmikikuulal/ )(11nkkj,.1,nki,.1对角线上
9、不能有对角线上不能有0， 计算之前先交换矩阵计算之前先交换矩阵A的元素，将主值交换到对角线上的元素，将主值交换到对角线上Copyright by Li Xinliang6bAx bLUX 回代过程回代过程bLy yUx bL1yYUX1nnnnnbbbbyyyyllllll.11112212213213231211111,ikkikiiylbybyiinikkikiinnnnuxuyxuyx/ )(,/1计算量：计算量： 分解分解O(n3/3), 回代回代 O(n2)优点：优点： 1） 重复求解重复求解 ， 时，仅需一次时，仅需一次LU分分解，计算量小；解，计算量小； 2） LU分解不破坏带状
10、稀疏矩阵的性质，可大幅减小计算量。分解不破坏带状稀疏矩阵的性质，可大幅减小计算量。11bAX 22bAX L带宽的带状矩阵带宽的带状矩阵：LU 分解分解: O(nL)回代：回代： O(nL)Copyright by Li Xinliang711.1.3 带状矩阵求解的追赶法带状矩阵求解的追赶法追赶法：等价于带状矩阵的追赶法：等价于带状矩阵的LU分解分解例：例： 三对角矩阵三对角矩阵bAX jjjjjjjdxcxbxa11一般项：一般项：边界项：边界项：12111dxcxbnnnnndxbxa122xutu21111112xuuutuunjnjnjnjnj追赶法追赶法jjjjBxAx1令：111
11、jjjjBxAxjjjjjjjjjjjjbAaBadxbAacx1111jjjjjjjjjjjjbAaBadBbAacA111,Step 1:Step 2: Step 3: Step 4: 111111/,/bdBbcAjjjjjjjjjjjjbAaBadBbAacA111,)/()(11nnnnnnnAabBadxjjjjBxAx1 计算量：计算量： 9n次次 （乘法）（乘法） A为固定值时，为固定值时， 3n次（乘法）次（乘法）简单易用，简单易用，计算量小计算量小Copyright by Li Xinliang811.2 代数方程组求解的迭代法代数方程组求解的迭代法bAx 直接法直接法迭代
12、法迭代法优点优点算法简便，准确（未知数少算法简便，准确（未知数少时）时）计算量小，误差容易控制计算量小，误差容易控制缺点缺点计算量大计算量大O(n3)舍入误差积累，不易控制舍入误差积累，不易控制快速收敛的算法设计较为复杂快速收敛的算法设计较为复杂11.2.1 Jocabi 及及Gauss-Seidel迭代迭代ininiibxaxaxa.2211解出对角元素解出对角元素nijjijijjijiiiixaxabax1111nijkjijijkjijiiikixaxabax1)(11)()1(1nijkjijijkjijiiikixaxabax1)(11)1()1(1Jocabi迭代迭代Gauss-
13、Seidel迭代迭代“对角占优对角占优”Copyright by Li Xinliang911.2.2 松弛迭代松弛迭代超松弛（超松弛（SOR）、亚松弛）、亚松弛Step1: 采用采用Jocabian 或或 Gauss-Seidel迭代产生新的值迭代产生新的值)1()(kikixxStep 2: 进行松弛进行松弛)()()1()()1(kikikikixxxx含义：含义： 改变步长改变步长1超松弛超松弛x精确解)1 (x)2(x“步子迈大一些步子迈大一些”，加快收敛，加快收敛1亚松弛亚松弛 “步子迈小一些步子迈小一些”，稳定性好，稳定性好收敛性：收敛性： 对角占优矩阵，对角占优矩阵，Jocab
14、ian及及Gauss-Seidel迭代可收敛迭代可收敛Copyright by Li Xinliang10举例：举例： Laplace方程的求解方程的求解),(),(2222yxguyxfyuxu 1 , 01 , 0),(yx) 1/(1Nyxjijijijijijijifuuuuuu,2,1,1,2, 1, 1222,1,1, 1, 14jijijijijijifuuuuu五点格式五点格式Jacobi迭代迭代Gauss-Seidel迭代迭代缺点缺点: 每迭代一步，信息只传递到周围网格点，每迭代一步，信息只传递到周围网格点，n很大时收敛较慢很大时收敛较慢4/ )(2,1,1, 1, 11,j
15、injinjinjinjinjifuuuuu4/ )(2,11,1,1, 1, 11,jinjinjinjinjinjifuuuuun+1nnnnn+1nn+1n+1Copyright by Li Xinliang11 对称对称Gauss-Seidel迭代迭代 （SGS）2,1,1, 1, 14jijijijijijifuuuuun+1nn+1n+1n+1nnnn+1n+14/ )(2,11,1,1, 1, 11,jinjinjinjinjinjifuuuuuStep 14/ )(2,11,11,1, 11, 11,jinjinjinjinjinjifuuuuuStep 2NjNi.2 , 1
16、;.2 , 11.1,; 1.1,NNjNNi特点：特点： 两次扫描，反复迭代两次扫描，反复迭代Copyright by Li Xinliang1211.2.3 交替方向迭代（交替方向迭代（ADI）方法）方法2,1,1, 1, 14jijijijijijifuuuuu例：例：2,1,1,*,*, 1*, 14jinjinjijijijifuuuuuStep 1:认为认为 已知已知 （使用上一步的（使用上一步的值），值）， 求解三对角方程求解三对角方程,得到中间步的值得到中间步的值njinjiuu1,1,Step 2: 代入中间步的值，求解三对角方程，得到代入中间步的值，求解三对角方程，得到n+
17、1步的值步的值 2,11,11,1,*, 1*, 14jinjinjinjijijifuuuuu三对角方程采用追赶法求解，效率较高三对角方程采用追赶法求解，效率较高在每一条线上采用直接法，信息快速传递，有利于收敛在每一条线上采用直接法，信息快速传递，有利于收敛Step 3: 重复以上两个步骤，直至收敛重复以上两个步骤，直至收敛因追赶法实际上是因追赶法实际上是LU分解法，因此又称分解法，因此又称LU-ADI方法方法n+1n+1n+1nn+1n+1n+1nnnCopyright by Li Xinliang1311.2.4 近似解法近似解法 LU-SGS方法方法2,1,1, 1, 14jijiji
18、jijijifuuuuu SGS 方法信息传递速度仍较慢，需要加速方法信息传递速度仍较慢，需要加速 ULDDUDDLUDL11)()()(近似近似LU分解分解2,1, 14jijijijifuuuStep 1:Step 2:jijijijiuuuu,1,1, 144优点：优点： 不含任何迭代过程不含任何迭代过程，两步扫描即可完成，效率很高；，两步扫描即可完成，效率很高；缺点：缺点： 近似近似LU分解，结果分解，结果不够准确不够准确4/ )(2,1, 1,jijijijifuuu4/ )(1, 1,jijijijiuuuuOK， 不迭代不迭代是是LU分解的分解的SGS方法，因方法，因此成为此成为
19、LU-SGS近似解法近似解法Copyright by Li Xinliang1411.2.5 加速收敛的多重网格法加速收敛的多重网格法周期边期边界),(2222yxfyuxu2,1,1, 1, 14jijijijijijifuuuuuGauss-Seidel迭代迭代4/ )(2,11,1,1, 1, 11,jinjinjinjinjinjifuuuuu)(,21jiykxkinnjie含义：含义： 线性系统，误差线性系统，误差 满足同样的方程满足同样的方程4/ )(11,1,1, 1, 11,njinjinjinjinjiuunn定义误差定义误差：1） 收敛速度的收敛速度的Fourier分析分
20、析iiiinneeeeG4/1增长（收敛）因子增长（收敛）因子ykxk21,含义： PPW/2极端高波数情况极端高波数情况2/1:,G迭代一次，误差减小一半迭代一次，误差减小一半极端低波数情况极端低波数情况1:0,G收敛速度趋近于收敛速度趋近于0 Copyright by Li Xinliang15策略：策略： 多重网格多重网格 粗网格加速低波数扰动收敛，细网格加速高波数收敛粗网格加速低波数扰动收敛，细网格加速高波数收敛ykxk21,PPW/2细网格粗网格使用多重网格法求解方程：使用多重网格法求解方程：4/ )(2,1,1, 1, 11,jinjinjinjinjinjifuuuuu2,1,1
21、, 1, 14jijijijijijifuuuuu迭代方程迭代方程：njinjinjiuu,1,4/ )(1,1, 1, 11,njinjinjinjinji 以以Jacobian迭代为例迭代为例修正方程修正方程Step 1: 在在细网格细网格上迭代一定步数（无需收上迭代一定步数（无需收敛），得到中间步的值敛），得到中间步的值Step 2: 将将修正项修正项 插值到插值到粗网格粗网格上，并上，并迭代求解迭代求解Step 3: 将求解后的将求解后的修正项修正项插值到插值到细网格，细网格， 并计算出细网格上新的值并计算出细网格上新的值Step 4: 重复重复Step 1-3 直到收敛直到收敛mum
22、修正项修正项Copyright by Li Xinliang16常用方法：常用方法： V 型型 及及 W 型迭代型迭代细网格细网格粗网格粗网格更粗网格更粗网格细网格细网格粗网格粗网格更粗网格更粗网格V 型迭代型迭代 W型迭代型迭代Copyright by Li Xinliang17Part 2 网格生成技术网格生成技术10.3 代数网格生成法代数网格生成法基本思路：基本思路： 通过代数方程计算出网格点的位置通过代数方程计算出网格点的位置优点：优点： 灵活、计算量小灵活、计算量小缺点：缺点： 光滑性差，光滑性差， 过于依赖人工过于依赖人工)(22xyy )(11xyy 如图：如图： 叶栅通道叶栅
23、通道已知计算域上边界（红线）及下边界已知计算域上边界（红线）及下边界（蓝线）的方程为：（蓝线）的方程为： 和和)(22xyy )(11xyy 则网格为：则网格为：)()()()() 1/() 1(,1,2,1,jjijijijijigxyxyxyyNix其中其中 可控制法向的疏密分布可控制法向的疏密分布)(jg) 1/() 1(Njj)(g均匀分布；均匀分布；)tanh()tanh()(bbg在下壁面处密集在下壁面处密集分布分布11)(bbeeg)(tanh2)12(tanh1)(11bbg上下壁面两侧加密上下壁面两侧加密Copyright by Li Xinliang18),(),(yxyx
24、10.3 椭圆形方程网格生成法椭圆形方程网格生成法ABCDEFABCDEF对于如图对于如图单联通单联通的计算域，可通过坐标变换的计算域，可通过坐标变换变换到图示矩形计算域变换到图示矩形计算域),(物理空间边界物理空间边界 计算空间边界计算空间边界物理空间内点物理空间内点 计算空间内点计算空间内点物理空间计算空间通常：通常： 给定边界点的对应关系给定边界点的对应关系 （代数方法）（代数方法） 通过求解方程获得内点的对应关系通过求解方程获得内点的对应关系方程的边值问题方程的边值问题椭圆型方程椭圆型方程 边值问题边值问题抛物型方程抛物型方程 双曲型方程双曲型方程初边值初边值问题问题椭圆型方程椭圆型方
25、程Copyright by Li Xinliang19通常：0022222222yxyx或),(),(22222222yxQyxyxPyx含义：含义： 给物理空间的每个点找到计算空间的对应位置。给物理空间的每个点找到计算空间的对应位置。注注： 由于拓扑的对应性，物理空间必须是由于拓扑的对应性，物理空间必须是单联通域单联通域如果是多联通的，可通过切割，形成单联通域如果是多联通的，可通过切割，形成单联通域Copyright by Li Xinliang20),(),(22222222yxQyxyxPyx求解求解1. 形式变换，改写成以形式变换，改写成以 为自变量为自变量),(),(),(yx便于进
26、行便于进行差分求解差分求解1001yyyyxxxxyyyyxxyxyyxyxxxx xyJyJx ,xyJyJx ,Jx yx yy xxyxy ( , ) / ( , )()1 222220 xxxxxxxxxxxxxxx()() 222220 xyxxxxxyyyyyyyy()()yyyxxx22 JPxQx2() JPyQy2()xy22 ()x xy yxy22Copyright by Li Xinliang21离散化：离散化： 中心差分中心差分 1yyyxxx22jijijixxxx, 1, 12jijijixxxx,1,1,22/ )(, 1, 1jijixxx2/ )(1,1,j
27、ijixxx JPxQx2() JPyQy2()xy22 ()x xy yxy22ijjijijijijijijijijijixCxCxCxCxC,1,1,1,1, 1, 1, 1, 1ijjijijijijijijijijijiyDyDyDyDyD,1,1,1,1, 1, 1, 1, 1离散方程离散方程jijijijijixxxxxx,1,1, 1, 14迭代求解：迭代求解： Jacobi, Seidel, SOR,LU-ADI, LU-SGS ,多重网格，多重网格， Copyright by Li Xinliang22Step 1: 确定边界网格确定边界网格通常采用代数方法生成通常采用代数
28、方法生成 一维网格，容易生成一维网格，容易生成注意：注意： 1）边界对应关系容易出错）边界对应关系容易出错 2）考虑网格的疏密分布）考虑网格的疏密分布 （翼型尾缘区，激波区，（翼型尾缘区，激波区， 近壁区近壁区)(,jijiyxStep 2: 利用上页的离散方程，利用上页的离散方程， 解出全部网格坐标解出全部网格坐标 )(,jijiyx不足：不足： 内部区域的网格分布不易控制内部区域的网格分布不易控制 无法做到指定区域网格加密无法做到指定区域网格加密 无法保证网格正交无法保证网格正交边界网格可控，边界网格可控，内部网格只能内部网格只能“听天由命听天由命”方案方案1： 源项源项P,Q为为0 ，
29、求解求解Laplace方程方程Copyright by Li Xinliang23方案方案2： 设定源项设定源项P,Q 求解求解Poison方程方程),(),(22222222yxQyxyxPyx源项源项P,Q对网格的影响对网格的影响 数值实验发现，在某点处加入点源数值实验发现，在某点处加入点源P: P0 使使 方向网格线发散方向网格线发散 123123123123123123点源点源PP0P0在某点处加入点源在某点处加入点源Q, 可对可对 方向方向网格线产生同样效果网格线产生同样效果1） 网格线的汇聚网格线的汇聚启发：启发： 在某条网格线上加入负的在某条网格线上加入负的源项，可令网格汇聚源项
30、，可令网格汇聚2/12211)()(exp)()exp()(),(lllLlllIiiiiidsignbcsignaP2/122 1 1 )()(exp)()exp()(),(lllLllljjjJjjdsignbcsignaQ使网格汇使网格汇聚于聚于llCopyright by Li Xinliang242） 边界网格的正交（并指定边界网格间距）边界网格的正交（并指定边界网格间距）)()1()()1(maxmax)()(),()()(),(dbcaeseqQerepP令：源项在边界处，内部衰减利用边界处网格的利用边界处网格的正交性正交性及及网格间距网格间距要求要求确定系数确定系数P和和Q122s指定值yyyxxx22 JPxQx2() JPyQy2()xy22 ()x xy yxy22基本思路： 以壁面线 处为例 1网格线正交网格线正交指定法向网指定法向网格间距格间距计算第计算第2层层网格线网格线 上上的坐标的坐标2通过差分计算边界处的通过差分计算边界处的yyyRxxxR2221如需要利用如需要利用 的信息，可用上一迭代步的值的信息，可用上一迭代步的值312)(RQxPxJ22)(RQyPyJ计算出边界处的计算出边界处的P,Q根据指数衰减原则给出全场的根据指数衰减原则给出全场的P,Q具体公式见：具体公式见： 傅德薰傅德薰 计算空气动力学计算空气动力学284-286
