计算流体力学(中科院力学所)_第9讲-有限体积法1.ppt

计算流体力学讲义第九讲有限体积法 1 李新亮lixl 力学所主楼219 82543801 知识点 1 讲义 课件上传至 流体中文网 流体论坛 CFD基础理论 讲课录像及讲义上传至网盘http cid CopyrightbyLiXinliang 有限体积法的基本概念 重构和反演迎风型有限体积法 Riemann求解器 Roe格式的新理解 近似Riemann解多维迎风型有限体积法 坐标旋转 CopyrightbyLiXinliang 2 知识回顾 1 差分方法的基本概念 差分格式 修正方程 相容性 收敛性 稳定性 LAX等价定理 2 精度分析 稳定性分析与分辨率分析 修正波数 Taylor分析 Fourier分析 修正波数 激波捕捉格式GVC NND Roe Godnov MUSCL TVD WENO Euler N S 方程的通量分裂逐点分裂 特征投影分裂 建议使用Roe平均 5 隐格式求解的LU SGS方法 要点 a 引入差量 方程线性化b 单边差分 隐式代数方程显式 推进 化 以一维为例 多维可直接推广 方法1 直接隐式离散 直接求解 非线性方程组 计算量大 方法2 差量化 线性化 已知项 线化微分方程 CopyrightbyLiXinliang 3 CopyrightbyLiXinliang 4 求解思路 如果直接离散 得到线性代数方程组 仍需求解 计算量大 多维情况 如果能单侧差分就好解了 多对角方程组 不好解 多维情况 中心 双侧 离散 如果单侧离散 单侧离散 可推进求解 免受解方程组之苦 真简单 CopyrightbyLiXinliang 5 可是 A有正有负 无法单侧差分化 还是个三对角的 奇思妙想 如果分成两个子步 各自用单侧值 就简单多了 强行单侧差分会不稳定的 近似LU分解 Step1 近似LU分解 Step2 均为递推求解 两次扫描 免受解方程组之苦 j 1 j j 1 j 以上描述适用于求解定常问题 求解非定常问题该过程可用于内迭代 迭代收敛后q趋于0 精度由右端项决定 CopyrightbyLiXinliang 6 9 1有限体积法入门 有限体积法主要优势 处理复杂网格 差分法处理复杂外形 坐标变换 坐标变换函数必须足够光滑 否则损失精度 实际问题 外形复杂 光滑的结构网格生成困难 CopyrightbyLiXinliang 7 9 1 1有限体积法的基本概念 实质 把几何信息包含于离散过程中 虽然简单 但有助于建立基本概念 j 1jj 1 j 1 2j 1 2 1 全离散型过程 含义 f在j 1 2点的值 注意与差分法的区别 在控制体上积分原方程 定义 空间平均 时间平均 精确推导 不含误差 提示 为区间内的空间及时间平均值 如果把它们理解为某点的值 会产生误差 CopyrightbyLiXinliang 8 积分 精确 重构 Reconstruction 有限差分法的离散 数值微分过程有限体积法的离散 数值积分过程 积分方程 离散化 反演 evolution 1 重构过程 A 零阶重构 假设分片常数 j 1 B 线性重构 假设分片线性函数 零阶重构与一阶重构示意图 j j 1 or or 或其他方法 C 更高阶的重构例如 分片二次函数 PPM WENO等 重构是有限体积的空间离散化过程 有多种方法 CopyrightbyLiXinliang 9 2 演化过程 以线性方程为例 需要得知时间演化信息 通常利用特征方程 若采用零阶重构 则 假设时间步长足够小 则方程为 等价于一阶迎风差分 Riemann解 CopyrightbyLiXinliang 10 若采用线性重构 若 Warming Beam Lax Wendroff 0阶重构 1阶精度线性重构 2阶精度 一维均匀网格的有限体积法等价于有限差分法 Euler方程 演化过程可通过Riemann解或近似Riemann解进行 CopyrightbyLiXinliang 11 2 半离散方法 全离散 积分方程 代数方程 守恒性好 但复杂 半离散 积分方程 常微分方程 简便 便于使用R K等成熟方法 仅空间积分 f在j 1 2点的值 仍需要使用周围点进行插值 通常无法精确计算 可采用近似值代替 等价于二阶中心差分 半离散 j 1jj 1 j 1 2j 1 2 重构 CopyrightbyLiXinliang 12 9 1 2一维Euler方程的迎风型有限体积法 j 1jj 1 j 1 2j 1 2 半离散 1 重构 控制体积 j 1 j j 1 左重构值 右重构值 选择不同的模板会得到不同的重构方案 向左偏的模板产生向右偏的模板产生 差分法 同一点的导数可使用向前差分和向后差分 根据特征方向选择之 例如 0阶重构1阶单边重构 根据特征方向 选择左通量或右通量 途径1 FVS 途径2 FDS CopyrightbyLiXinliang 13 2 分裂方法 1 FVS方法 流通矢量分裂 逐点分裂 具体方法 Steger Warming分裂Lax Friedrichs分裂VanLeer分裂 Liou Steffen分裂 压力项与其他项分开 AUSM类格式的基础 根据当地Mach数分裂 保证的Jocabian阵特征值为正 的为负 正通量 向左偏斜重构 负通量 向右偏斜重构 偏重向上游 与迎风差分法类似 网格基 或权重 偏重上游 差分 有限体积都可使用 一个参数 反映全部特征 CopyrightbyLiXinliang 14 小知识 Liou Steffen分裂 对流项 压力项 思路 决定特征的关键参数 当地Mach数 超音速 x 方向 超音速 x 方向 因此 对Mach数进行分裂更为简洁 显然 参考文献 Toro RiemannSolversandNumericalMethodsforFluidDynamics section8 4 4Liou TenYearsinthemakingAUSMfamily NASATM 2001 210977 类似VanLeer分裂 但压力单独处理 M 保证光滑过渡 M 1 CopyrightbyLiXinliang 15 3 FDS方法 通量差分分裂 特征投影分裂 1 利用精确Riemann解 Godnov格式 目的 j 1jj 1 j 1 2j 1 2 控制体积 j 1 j j 1 左重构值 右重构值 1 精确求解Riemann问题 2 精度 取决于重构的精度 原则上可任意阶 差分法 Godnov格式使用分片常数 精度1阶有限体积法 先重构 再解Riemann问题 可高阶 精确Riemann解 见本讲座第2讲 需迭代求解 计算量大 近似Riemann解 整体思路 先重构自变量 两种方案得到 再求解Riemann问题 或用FVS 得到通量的方法通常称为MUSCL方法 CopyrightbyLiXinliang 16 差分法与有限体积法区别与联系 二阶迎风 FVS为例 差分 有限体积 差分 通常做法 直接插值通量fi 1 2 有限体积 先插值自变量U 然后计算通量f 先插值自变量 再计算通量的方法 称为MUSCL类方法 是有限体积法的常用方法 差分法也可以用 单侧重构 以避免跨过激波 还可使用FDS方法 重构后求解Riemann问题 当f f U 连续时 对f插值与对U插值精度相同 称为数值流通量 的含义 CopyrightbyLiXinliang 17 重要概念澄清 重构与插值 A 有限差分法 j 1 2 切线 j 1 2 j j 1 注意 与f在xj 1 2点的值含义不同 用周围几个点的值计算的过程称为 重构 不能理解为用来插值 记号确实容易混淆 让人容易联想起 记为更好些 否则 最高只能达到2阶精度了 是控制体内的平均值 称为数值流通量 的含义 CopyrightbyLiXinliang 18 重要概念澄清 重构与插值 B 有限体积法 j 1 2 j 1 2 确实为f在xj 1 2点的值 通常做法 1 用计算出2 u在xj 1 2点的值 关键 是用计算 称为重构 而不是用计算 是标准的插值 否则最高也只能达到2阶精度 19 概念 MUSCL与非MUSC类方法 j 1 2 切线 j 1 2 j 1 差分 有限体积 方法1 非MUSCL类 直接利用周围几个点的函数值或 直接计算 或 如何计算或 方法2 MUSCL类 利用周围几个点的自变量值 或 计算出 或 然后再计算 或 当f f u 是连函数时 二者精度相同 f的误差与u的误差同阶 CopyrightbyLiXinliang 20 2 近似Riemann解 例 Roe格式 与差分法的Roe格式形式相同 理解 近似Riemann解 Euler方程常系数线性化解 u f u uL uR uRoe 利用Roe平均 刚好是左右两点间的平均增长率 实现了常系数线性化 常系数双曲方程组 易解 思路 用平均增长率矩阵取代瞬时增长率矩阵A 不但实现了线性化 而且实现了常系数化 利用二次齐函数的性质 可找到了Roe点 即Roe平均点 该点处的增长率刚好等于平均增长率 Roe平均 常系数化 线性化 常系数线性单波方程的Riemann问题 太简单了 21 常系数方程组的Riemann问题 解耦了的单波方程 有精确解 初值 CopyrightbyLiXinliang 22 解为 线性化条件 并利用齐函数性质 与差分法的Roe格式相同 还有各种其他类型的近似Riemann解 今后介绍 CopyrightbyLiXinliang 23 9 1 3多维问题的有限体积法 二维问题 一维Riemann问题 坐标选取不当 变为 二维 Riemann问题 x y 差分法 独立计算只考虑各自的特征方向 由于非线性 实际 二维 特征方向并非x y方向特征量的线性组合 特征方向计算不严格 带来误差 差分方法 多维情况 特征理论复杂 通常x y方向独立计算 转化为x方向与y方向的两个一维问题 逐点分裂 特征投影分裂 完全按照一维情况独立处理 局部坐标旋转 差分算法设计造成局部旋转困难 差分法的多维处理方法 1 小知识 差分方法如何处理高维问题的 优缺点 优点 简单缺点 特征方向计算不准 CopyrightbyLiXinliang 24 2 二维有限体积方法的离散过程 在以某节点为中心的控制体上积分 i j k 非结构网格的控制体 i 1 j i 1 j i j 1 i j 1 k3 k1 k2 k4 k5 结构网格的控制体 x y n 体积平均 控制体边界垂直于节点连线 也可选其他方式 垂直平分线 n 1 建立控制体 2 在控制体上积分 离散方程 重构 由节点上平均值给出函数分布 最终给出通量 表示第m个界面上的值 1 重构 两种不同的重构方案 向左偏及向右偏 给出两种结果 及 CopyrightbyLiXinliang 25 i j i 1 j i 1 j i j 1 i j 1 n 左重构 右重构 2 由左右重构得到的自变量 和给出通量方案A FVS方案B 解Riemann问题 常用 3 二维迎风型有限体积法 例如 0阶重构 线性重构 用i i 1点的值插i 1 2点的值 网格剧烈变化时 应当用实际坐标插值 用i i 1点的值插i 1 2点的值 x y 看似二维Riemann问题 其实是一维的 坐标旋转一下就行了 CopyrightbyLiXinliang 26 x y x y 通常 进行坐标旋转 旋转q角后的坐标系 x y 性质 Euler方程的旋转不变性 形式上完全不变 仅需把u v x y换成u v x y 即可 其中 旋转矩阵 旋转q角 矩阵表示 CopyrightbyLiXinliang 27 i j i 1 j i 1 j i j 1 i j 1 左重构 右重构 局部坐标系 x y x y 旋转q角后的坐标系 x y 习题 设n为平行x 轴的向量 试证明 证明 坐标旋转 标量不变向量的模不变 CopyrightbyLiXinliang 28 i j i 1 j i 1 j i j 1 i j 1 左重构 右重构 于是 x y x y 旋转q角后的坐标系 x y 其中下标m表示控制体第m个面 线 表示该面的面积 长度 于是 问题转化为求控制面上的 这个量有两个重构方案 方法1 FVS 方法2 需要求解Riemann问题 旋转后 转化为 扩展的 1维Riemann问题 CopyrightbyLiXinliang 29 解释 扩展的 一维Riemann问题 x y x y 旋转q角后的坐标系 x y 问题本身是一维的 所有变量都只沿着x 方向分布 沿y 方向均匀允许有y 方向的速度v 比纯一维问题多一个变量 v 的存在对流动的一维性质无任何影响 举例 Sod激波管问题 一维 如果在沿y方向匀速运动的坐标系中观察 则方程为 扩展的一维问题 但不影响其一维性质 坐标系沿y方向匀速运动 x y 可用精确Riemann解 也可用Roe等近似解 CopyrightbyLiXinliang 30 二维迎风型有限体积法求解步骤1 对n时刻的平均量进行重构 给出控制面上的左 右重构值 2 将以上值旋转到 每个 控制面法向的局部坐标系下 3 求解上述 扩展的 一维Riemann问题 给出后续时刻控制面上的值4 利用积分型方程 计算下一时刻的平均量 i j i 1 j i 1 j i j 1 i j 1 左重构 右重构 0阶重构 1阶重构 线性重构 更复杂的重构 WENO等 下标m指的是第m个控制面上的值 CopyrightbyLiXinliang 31 知识回顾 Riemann问题精确解 Riemann问题 问题描述 初始时刻 物理量分布存在单个间断 间断两侧物理量为常数 求解思路 采用积分方程单个间断 且间断两侧物理量为常数情况下 积分方程转化为代数方程 代数方程 质量 动量 能量守恒 计算出 将与这三个值进行比较 判断会产生的情况 具体见下图 CopyrightbyLiXinliang 32 Riemann问题的具体计算步骤 全流场 1 判断可能会出现的情况 五种情形之一 a 定义函数 b 进行判断 情况5 情况4 情况3 情况1 情况5 情况4 情况2 情况1 单调增函数 性质很好 CopyrightbyLiXinliang 33 2 求解中心区的压力和速度 单未知数的代数方程 迭代求解 例如Newton法 F p 性质好 求解不困难 3 确定中心区接触间断两侧的密度以及左 右波传播的速度a 左波为激波的情况 情况1 3 b 左波为稀疏波的情况 情况2 4 5 中心区接触间断左侧的物理量 膨胀波的波头及波尾速度 激波的传播速度 对于情况 5 波尾速度为 中心区为真空 音速无定义 改由该式计算 CopyrightbyLiXinliang 34 c 右波为激波的情况 情况1 2 中心区接触间断右侧的物理量 b 右波为稀疏波的情况 情况2 4 5 4 计算稀疏波区域的值 如果有稀疏波的话 a 左稀疏波b 右稀疏波 情况2 4 情况5 注意 教科书32页c的公式有误 CopyrightbyLiXinliang 35 有限体积法 扩展的 Riemann问题的计算方法 中心线x 0处 迎风型有限体积法 需要求解 扩展的 一维Riemann问题 x y 物理问题分析 所有物理量均沿x方向一维分布 沿y方向均匀分布 仅需计算t时刻x 0