中科院计算流体力学最新讲义CFD第讲求代数方程组及网格生成.ppt

计算流体力学讲义 第九讲 代数方程的求解及网格生成 李新亮 ；力学所主楼219； 82543801,知识点：,1,讲义、课件上传至 (流体中文网） - “流体论坛” -“ CFD基础理论 ” 讲课录像及讲义上传至网盘 /browse.aspx/.Public,Copyright by Li Xinliang,代数方程组的求解 网格生成,Copyright by Li Xinliang,2,知识回顾： 有限体积法基本流程,无粘项常用方法 （流过AB边的通量）： a. 利用周围点的值，计算出(I+1/2,J) 点处的物理量； 直接利用“差分格式” b. 利用该处的物理量，计算出流过AB边的流通量 迎风型方法需利用“通量分裂技术”,FVS类：,FDS类：,利用Riemann解,Reimann解： Godunov, Roe, HLL, HLLC,利用坐标变换，转化为一维Riemann问题,i,j,i-1,j,i+1,j,i,j+1,i,j+1,知识回顾： 隐格式,代数方程组的求解,Copyright by Li Xinliang,4,微分方程（组）,代数方程组,数值解（离散解）,差分 有限体积,大部分计算量,9.1.1 Gauss 消去法,消元,为了计算稳定，通常使用主元消去法 列主元消去法； 全主元消去法 计算量： 乘法： 加法：,9.1 代数方程组求解的直接法,优点： 简单精确，缺点：计算量大,上三角 矩阵,Copyright by Li Xinliang,5,9.1.2 LU分解法,0,0,Step 1,Step 2,Step k,对角线上不能有0， 计算之前先交换矩阵A的元素，将主值交换到对角线上,Copyright by Li Xinliang,6,回代过程,计算量： 分解O(n3/3), 回代 O(n2),优点： 1） 重复求解 ， 时，仅需一次LU分解，计算量小； 2） LU分解不破坏带状稀疏矩阵的性质，可大幅减小计算量。,L带宽的带状矩阵：,LU 分解: O(nL) 回代： O(nL),Copyright by Li Xinliang,7,9.1.3 带状矩阵求解的追赶法,追赶法：等价于带状矩阵的LU分解,例： 三对角矩阵,一般项：,边界项：,追赶法,令：,Step 1: Step 2: Step 3: Step 4:,计算量： 9n次 （乘法） A为固定值时， 3n次（乘法）,简单易用，计算量小,Copyright by Li Xinliang,8,9.2 代数方程组求解的迭代法,9.2.1 Jocabi 及Gauss-Seidel迭代,解出对角元素,Jocabi迭代,Gauss-Seidel迭代,“对角占优”,Copyright by Li Xinliang,9,9.2.2 松弛迭代超松弛（SOR）、亚松弛,Step1: 采用Jocabian 或 Gauss-Seidel迭代产生新的值,Step 2: 进行松弛,含义： 改变步长,超松弛,精确解,“步子迈大一些”，加快收敛,亚松弛,“步子迈小一些”，稳定性好,收敛性： 对角占优矩阵，Jocabian及Gauss-Seidel迭代可收敛,Copyright by Li Xinliang,10,举例： Laplace方程的求解,五点格式,Jacobi迭代,Gauss-Seidel迭代,缺点: 每迭代一步，信息只传递到周围网格点，n很大时收敛较慢,n+1,n,n,n,n,n+1,n,n+1,n+1,Copyright by Li Xinliang,11,对称Gauss-Seidel迭代 （SGS）,n+1,n,n+1,n+1,n+1,n,n,n,n+1,n+1,Step 1,Step 2,特点： 两次扫描，反复迭代,Copyright by Li Xinliang,12,9.2.3 交替方向迭代（ADI）方法,例：,Step 1:认为 已知 （使用上一步的值）， 求解三对角方程,得到中间步的值,Step 2: 代入中间步的值，求解三对角方程，得到n+1步的值,三对角方程采用追赶法求解，效率较高 在每一条线上采用直接法，信息快速传递，有利于收敛,Step 3: 重复以上两个步骤，直至收敛,因追赶法实际上是LU分解法，因此又称LU-ADI方法,n+1,n+1,n+1,n,n+1,n+1,n+1,n,n,n,Copyright by Li Xinliang,13,9.2.4 近似解法 LU-SGS方法,SGS 方法信息传递速度仍较慢，需要加速,近似LU分解,Step 1:,Step 2:,优点： 不含任何迭代过程，两步扫描即可完成，效率很高； 缺点： 近似LU分解，结果不够准确,OK， 不迭代,是LU分解的SGS方法，因此成为LU-SGS近似解法,Copyright by Li Xinliang,14,9.2.5 加速收敛的多重网格法,Gauss-Seidel迭代,含义： 线性系统，误差 满足同样的方程,定义误差：,1） 收敛速度的Fourier分析,增长（收敛）因子,含义：,极端高波数情况,迭代一次，误差减小一半,极端低波数情况,收敛速度趋近于0,Copyright by Li Xinliang,15,策略： 多重网格 粗网格加速低波数扰动收敛，细网格加速高波数收敛,细网格,粗网格,使用多重网格法求解方程：,迭代方程：,以Jacobian迭代为例,修正方程,Step 1: 在细网格上迭代一定步数（无需收敛），得到中间步的值 Step 2: 将修正项 插值到粗网格上，并迭代求解 Step 3: 将求解后的修正项插值到细网格， 并计算出细网格上新的值 Step 4: 重复Step 1-3 直到收敛,修正项,Copyright by Li Xinliang,16,常用方法： V 型 及 W 型迭代,细网格,粗网格,更粗网格,细网格,粗网格,更粗网格,V 型迭代 W型迭代,Copyright by Li Xinliang,17,9.2.6 共轭梯度法,1. 求解对称正定矩阵的共轭梯度法,化代数方程组问题为极值问题 （设A为对称正定矩阵）,的最小值问题,可用最速下降法之类方法求解,例：,的分布,最小值点（0,1）,Copyright by Li Xinliang,18,求极值问题的最速下降法,思路： 沿当地梯度方向前进,1）根据当前位置 ， 计算当地梯度方向：,2） 沿该方向前进，使得 达到极小值,方向,步长,“残余向量”,特点： 沿当地梯度方向前进，直到不能前进为止，然后以按照新的梯度方向前进；相邻路径方向正交； 缺点：局部最速下降路线并非全局最速下降路线，因而收敛速度并非最优。,“最速下降法”示意图,Copyright by Li Xinliang,19,最速下降法的改进：共轭梯度法,最速下降法：每步在一维空间求最优解：,改进：在二维空间寻求最优解 （不再沿当地梯度方向）,旧路径方向及当前梯度方向所张成的二维空间,寻找该平面内的极小值,解出,新线路：,修正,最终得：,特点： 相邻两次方向关于矩阵A正交 （称为共轭）。,Copyright by Li Xinliang,20,2. 求解一般非奇异方程组的共轭梯度法,A为一般满秩阵（非对称正定阵）的情况,正则化方法：,对称正定阵,Step 1: 设置初值,Step 2: 迭代推进,直到残余向量足够小为止,Copyright by Li Xinliang,21,9.3 代数网格生成法,基本思路： 通过代数方程计算出网格点的位置 优点： 灵活、计算量小 缺点： 光滑性差， 过于依赖人工,如图： 叶栅通道,已知计算域上边界（红线）及下边界（蓝线）的方程为： 和,则网格为：,其中 可控制法向的疏密分布,均匀分布；,在下壁面处密集分布,上下壁面两侧加密,Copyright by Li Xinliang,22,9.4 椭圆形方程网格生成法,A,B,C,D,E,F,A,B,C,D,E,F,对于如图单联通的计算域，可通过坐标变换,变换到图示矩形计算域,物理空间边界 计算空间边界 物理空间内点 计算空间内点,物理空间,计算空间,通常： 给定边界点的对应关系 （代数方法） 通过求解方程获得内点的对应关系,方程的边值问题,椭圆型方程 边值问题 抛物型方程 双曲型方程,初边值问题,椭圆型方程,Copyright by Li Xinliang,23,通常：,或,含义： 给物理空间的每个点找到计算空间的对应位置。,注： 由于拓扑的对应性，物理空间必须是单联通域,如果是多联通的，可通过切割，形成单联通域,Copyright by Li Xinliang,24,求解,1. 形式变换，改写成以 为自变量,便于进行差分求解,Copyright by Li Xinliang,25,离散化： 中心差分,离散方程,迭代求解： Jacobi, Seidel, SOR,LU-ADI, LU-SGS ,多重网格， ,Copyright by Li Xinliang,26,Step 1: 确定边界网格,通常采用代数方法生成 一维网格，容易生成 注意： 1）边界对应关系容易出错 2）考虑网格的疏密分布 （翼型尾缘区，激波区， 近壁区),Step 2: 利用上页的离散方程， 解出全部网格坐标,不足： 内部区域的网格分布不易控制 无法做到指定区域网格加密 无法保证网格正交,边界网格可控，内部网格只能“听天由命”,方案1： 源项P,Q为0 ， 求解Laplace方程,Copyright by Li Xinliang,27,方案2： 设定源项P,Q 求解Poison方程,源项P,Q对网格的影响 数值实验发现，在某点处加入点源P: P0 使 方向网格线发散,点源P,P0,P0,在某点处加入点源Q, 可对 方向网格线产生同样效果,1） 网格线的汇聚,启发： 在某条网格线上加入负的源项，可令网格汇聚,使网格汇聚于,Copyright by Li Xinliang,28,2） 边界网格的正交（并指定边界网格间距）,令：,源项在边界处， 内部衰减,利用边界处网格的正交性及网格间距要求确定系数P和Q,指定值,基本思路： 以壁面线 处为例,网格线正交,指定法向网格间距,计算第2层网格线 上的坐标,通过差分计算边界处的,如需要利用 的信息，可用上一迭代步的值,计算出边界处的P,Q,根据指数衰减原则给出全场的P,Q,具体公式见： 傅德薰 计算空气动力学284-286,Copyright by Li Xinliang,29,习题 9.1 网格生成,通过解椭圆型方程生成NACA0012翼型的网格 要求： 推荐采用图示的C型网格， 网格点不限； 外边界的位置不限； 可求解无源