CFD第9讲求代数方程组及网格生成

计算流体力学讲义

第九讲代数方程的求解及网格生成李新亮lixl@；力学所主楼219；82543801知识点：

1讲义、课件上传至(流体中文网）->“流体论坛”->“CFD基础理论

”讲课录像及讲义上传至网盘/browse.aspx/.PublicCopyrightbyLiXinliang代数方程组的求解网格生成CopyrightbyLiXinliang2知识回顾：有限体积法基本流程无粘项常用方法（流过AB边的通量）：a.利用周围点的值，计算出(I+1/2,J)点处的物理量；

直接利用“差分格式”b.利用该处的物理量，计算出流过AB边的流通量

迎风型方法需利用“通量分裂技术”FVS类：FDS类：利用Riemann解Reimann解：Godunov,Roe,HLL,HLLC利用坐标变换，转化为一维Riemann问题i,ji-1,ji+1,ji,j+1i,j+1知识回顾：隐格式代数方程组的求解CopyrightbyLiXinliang4微分方程（组）代数方程组数值解（离散解）差分有限体积大部分计算量9.1.1Gauss消去法消元为了计算稳定，通常使用主元消去法列主元消去法；全主元消去法计算量：乘法：

加法：9.1代数方程组求解的直接法优点：简单精确，缺点：计算量大上三角矩阵CopyrightbyLiXinliang59.1.2LU分解法……00Step1Step2Stepk对角线上不能有0，计算之前先交换矩阵A的元素，将主值交换到对角线上CopyrightbyLiXinliang6回代过程计算量：分解O(n3/3),回代O(n2)优点：1）重复求解，……时，仅需一次LU分解，计算量小；2）LU分解不破坏带状稀疏矩阵的性质，可大幅减小计算量。L带宽的带状矩阵：LU分解:O(nL)回代：O(nL)CopyrightbyLiXinliang79.1.3带状矩阵求解的追赶法追赶法：等价于带状矩阵的LU分解例：三对角矩阵一般项：边界项：追赶法令：Step1:Step2:

Step3:Step4:

计算量：9n次（乘法）A为固定值时，3n次（乘法）简单易用，计算量小CopyrightbyLiXinliang89.2代数方程组求解的迭代法直接法迭代法优点算法简便，准确（未知数少时）计算量小，误差容易控制缺点计算量大O(n3)舍入误差积累，不易控制快速收敛的算法设计较为复杂9.2.1Jocabi及Gauss-Seidel迭代解出对角元素Jocabi迭代Gauss-Seidel迭代“对角占优”CopyrightbyLiXinliang99.2.2松弛迭代——超松弛（SOR）、亚松弛Step1:采用Jocabian或Gauss-Seidel迭代产生新的值Step2:进行松弛含义：改变步长超松弛精确解“步子迈大一些”，加快收敛亚松弛“步子迈小一些”，稳定性好收敛性：

对角占优矩阵，Jocabian及Gauss-Seidel迭代可收敛CopyrightbyLiXinliang10举例：Laplace方程的求解五点格式Jacobi迭代Gauss-Seidel迭代缺点:每迭代一步，信息只传递到周围网格点，n很大时收敛较慢n+1nnnnn+1nn+1n+1CopyrightbyLiXinliang11对称Gauss-Seidel迭代（SGS）n+1nn+1n+1n+1nnnn+1n+1Step1Step2特点：两次扫描，反复迭代CopyrightbyLiXinliang129.2.3交替方向迭代（ADI）方法例：Step1:认为已知（使用上一步的值），求解三对角方程,得到中间步的值Step2:代入中间步的值，求解三对角方程，得到n+1步的值

三对角方程采用追赶法求解，效率较高在每一条线上采用直接法，信息快速传递，有利于收敛Step3:重复以上两个步骤，直至收敛因追赶法实际上是LU分解法，因此又称LU-ADI方法n+1n+1n+1nn+1n+1n+1nnnCopyrightbyLiXinliang139.2.4近似解法——LU-SGS方法

SGS方法信息传递速度仍较慢，需要加速

近似LU分解Step1:Step2:优点：不含任何迭代过程，两步扫描即可完成，效率很高；缺点：近似LU分解，结果不够准确OK，

不迭代是LU分解的SGS方法，因此成为LU-SGS近似解法CopyrightbyLiXinliang149.2.5加速收敛的多重网格法Gauss-Seidel迭代含义：线性系统，误差满足同样的方程定义误差：1）收敛速度的Fourier分析增长（收敛）因子含义：极端高波数情况迭代一次，误差减小一半极端低波数情况收敛速度趋近于0CopyrightbyLiXinliang15策略：多重网格——粗网格加速低波数扰动收敛，细网格加速高波数收敛细网格粗网格使用多重网格法求解方程：迭代方程：以Jacobian迭代为例修正方程Step1:在细网格上迭代一定步数（无需收敛），得到中间步的值Step2:将修正项插值到粗网格上，并迭代求解Step3:将求解后的修正项插值到细网格，并计算出细网格上新的值Step4:重复Step1-3直到收敛修正项CopyrightbyLiXinliang16常用方法：V型及W型迭代细网格粗网格更粗网格细网格粗网格更粗网格V型迭代W型迭代CopyrightbyLiXinliang179.2.6共轭梯度法1.求解对称正定矩阵的共轭梯度法化代数方程组问题为极值问题（设A为对称正定矩阵）的最小值问题可用最速下降法之类方法求解例：的分布最小值点（0,1）CopyrightbyLiXinliang18求极值问题的最速下降法思路：沿当地梯度方向前进1）根据当前位置，计算当地梯度方向：

2）沿该方向前进，使得达到极小值方向步长“残余向量”特点：沿当地梯度方向前进，直到不能前进为止，然后以按照新的梯度方向前进；相邻路径方向正交；缺点：局部最速下降路线并非全局最速下降路线，因而收敛速度并非最优。“最速下降法”示意图CopyrightbyLiXinliang19最速下降法的改进：共轭梯度法最速下降法：每步在一维空间求最优解：改进：在二维空间寻求最优解（不再沿当地梯度方向）旧路径方向及当前梯度方向所张成的二维空间寻找该平面内的极小值解出新线路：修正最终得：特点：相邻两次方向关于矩阵A正交（称为共轭）。CopyrightbyLiXinliang202.求解一般非奇异方程组的共轭梯度法A为一般满秩阵（非对称正定阵）的情况正则化方法：对称正定阵Step1:设置初值

Step2:迭代推进直到残余向量足够小为止CopyrightbyLiXinliang219.3代数网格生成法基本思路：通过代数方程计算出网格点的位置优点：灵活、计算量小缺点：光滑性差，过于依赖人工如图：叶栅通道已知计算域上边界（红线）及下边界（蓝线）的方程为：和则网格为：其中可控制法向的疏密分布均匀分布；在下壁面处密集分布上下壁面两侧加密CopyrightbyLiXinliang229.4椭圆形方程网格生成法ABCDEFA’B’C’D’E’F’对于如图单联通的计算域，可通过坐标变换变换到图示矩形计算域物理空间边界

计算空间边界物理空间内点计算空间内点物理空间计算空间通常：给定边界点的对应关系（代数方法）

通过求解方程获得内点的对应关系方程的边值问题椭圆型方程

边值问题抛物型方程双曲型方程初边值问题椭圆型方程CopyrightbyLiXinliang23通常：或含义：给物理空间的每个点找到计算空间的对应位置。注：

由于拓扑的对应性，物理空间必须是单联通域如果是多联通的，可通过切割，形成单联通域CopyrightbyLiXinliang24求解1.形式变换，改写成以为自变量便于进行差分求解

CopyrightbyLiXinliang25离散化：中心差分离散方程迭代求解：Jacobi,Seidel,SOR,LU-ADI,LU-SGS,多重网格，……CopyrightbyLiXinliang26Step1:确定边界网格通常采用代数方法生成——一维网格，容易生成注意：1）边界对应关系容易出错2）考虑网格的疏密分布（翼型尾缘区，激波区，近壁区……)Step2:利用上页的离散方程，解出全部网格坐标

不足：内部区域的网格分布不易控制

无法做到指定区域网格加密

无法保证网格正交边界网格可控，内部网格只能“听天由命”方案1：源项P,Q为0，求解Laplace方程CopyrightbyLiXinliang27方案2：设定源项P,Q求解Poison方程源项P,Q对网格的影响

数值实验发现，在某点处加入点源P:P<0使方向网格线汇聚P>0使方向网格线发散

点源PP>0P<0在某点处加入点源Q,可对方向网格线产生同样效果1）网格线的汇聚启发：在某条网格线上加入负的源项，可令网格汇聚使网格汇聚于CopyrightbyLiXinliang282）边界网格的正交（并指定边界网格间距）令：源项在边界处，内部衰减利用边界处网格的正交性及网格间距要求确定系数P和Q指定值基本思路：以壁面线处为例

网格线正交指定法向网格间距计算第2层网格线上的坐标通过差分计算边界处的如需要利用的信息，可用上一迭代步的值计算出边界处的P,Q根据指数衰减原则给出全场的P,Q具体公式见：傅德薰《计算空气动力学》284-286CopyrightbyLiXinliang29习题9.1网格生成

通过解椭圆型方程生成NACA0012翼型的网格要求：推荐采用图示的C型网格，