采用有限差分法求解一维欧拉方程：三阶R-K时间推进——三阶迎风差分（Steger-Warming）

Wang,Kui

已于 2024-01-05 13:38:33 修改

阅读量1.3k

点赞数 10

文章标签：学习线性代数制造笔记

于 2024-01-05 13:34:23 首次发布

本文链接：https://blog.csdn.net/ocean1702/article/details/135396559

版权

#最近看了李新亮老师的OpenCFD一维差分格式验证的求解器代码，对通量的求法有了一些新的认识，过去自己使用的要么是一阶的，精度差，要么就没考虑到激波问题。在实际求解器中，会使用一些方法来构造高精度的差分格式并避免激波问题，这就需要学习各种差分格式的构造方程，于是有此文章，全当记录。#

一、问题描述

针对一维Sod激波管问题：

其控制方程可以写为：

$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F(U})}{\partial x}=0$

$\mathbf{U}=\begin{bmatrix} \rho \\ \rho u \\ E \end{bmatrix}$ $\mathbf{F(U)}=\begin{bmatrix} \rho \\ \rho u^{2}+p \\ (E+p)u \end{bmatrix}$

二、求解思路

Step 1：初始化。划分网格，初始化 $u,p ,\rho$ 。

Step 2：计算通量。采用三阶迎风格式构造构造通量。由于迎风格式需要知道方向，即正负，这里已经把通量分为了正值和负值，在第三步会处理。

$\phi _{j}=\frac{(F_{j+1/2}^{+}+F_{j+1/2}^{-})-(F_{j-1/2}^{+}+F_{j-1/2}^{-})}{\Delta x}$

对于三阶迎风格式：

$F_{j+1/2}^{+}=-\frac{1}{6}F_{j-1}^{+}+\frac{5}{6}F_{j}^{+}+\frac{2}{6}F_{j+1}^{+} ; F_{j+1/2}^{-}=\frac{2}{6}F_{j-1}^{-}+\frac{5}{6}F_{j}^{-}-\frac{1}{6}F_{j+1}^{-}$

$F_{j-1/2}^{+}=-\frac{1}{6}F_{j-2}^{+}+\frac{5}{6}F_{j-1}^{+}+\frac{2}{6}F_{j}^{+} ; F_{j-1/2}^{-}=\frac{2}{6}F_{j-2}^{-}+\frac{5}{6}F_{j-1}^{-}-\frac{1}{6}F_{j}^{-}$

Step 3：通量分裂处理。为了采用高阶的迎风格式，需要对通量进行处理，这里采用Steger-Warming方法分裂通量，避免出现间断。其思路是将通量分为正值和负值，然后将二者相加，即得到了总通量，即：

$F_{j+1/2}=F_{j+1/2}^{+}+F_{j+1/2}^{ -}$

由上一步可知：

$F_{j+1/2}^{+}=-\frac{1}{6}F_{j-1}^{+}+\frac{5}{6}F_{j}^{+}+\frac{2}{6}F_{j+1}^{+} ; F_{j+1/2}^{-}=\frac{2}{6}F_{j-1}^{-}+\frac{5}{6}F_{j}^{-}-\frac{1}{6}F_{j+1}^{-}$

$F_{j-1/2}^{+}=-\frac{1}{6}F_{j-2}^{+}+\frac{5}{6}F_{j-1}^{+}+\frac{2}{6}F_{j}^{+} ; F_{j-1/2}^{-}=\frac{2}{6}F_{j-2}^{-}+\frac{5}{6}F_{j-1}^{-}-\frac{1}{6}F_{j}^{-}$

因此我们需要求出 $F_{j}^{\pm }$ ，过程网上很多人都有介绍如（https://zhuanlan.zhihu.com/p/148377847），这里我们不做推导，只给出结果：

$F_{j}^{\pm }=\frac{\rho_{j}}{2\gamma }\begin{bmatrix} 2(\gamma-1)\lambda _{j,1}^{\pm }+\lambda _{j,2}^{\pm }+\lambda _{j,3}^{\pm }\\ 2(\gamma-1)\lambda _{j,1}^{\pm }u_{j}+\lambda _{j,2}^{\pm }(u_{j}-c_{j})+\lambda _{j,3}^{\pm }(u_{j}+c_{j}) \\ (\gamma-1)\lambda _{j,1}^{\pm }u_{j}^{2}+\lambda _{j,2}^{\pm }(u_{j}-c_{j})^{2}/2+\lambda _{j,3}^{\pm }(u_{j}+c_{j})/2+w_{j}^{\pm} \end{bmatrix}$

$w_{j}^{\pm}=\frac{(3-\gamma)(\lambda _{j,2}^{\pm }+\lambda _{j,3}^{\pm })c_{j}^{2}}{2(\gamma-1)}$

式中：

$\lambda _{j,k}^{\pm }=(\lambda _{j,k}\pm \left | \lambda _{j,k} \right |)/2\; \; \; (k=1,2,3\; ;j=1,2,3,4,5,...,N)$

为了避免导数在间断处不连续，一般会写成这种形式：

$\lambda _{j,k}^{\pm }=(\lambda _{j,k}\pm \sqrt{\lambda_{j,k}^{2}+\varepsilon ^{2}}\,\,)/2 \, \,\, \, \,\, \, \,\varepsilon =10^{-3} \, \, \, to\, \, \,10^{-6}$

其中 $\lambda _{j,1}=u_{j}\; \; \; \; \lambda _{j,2}=u_{j}-c_{j}\; \; \; \;\lambda _{j,3}=u_{j}+c_{j}$

Step 4 :时间推进求解。这里采用三阶Runge-Kutta法进行时间推进。

经过我们上述的处理，我们可以得到常微分方程如下：

$\frac{\partial U_{j}}{\partial t}+\phi _{j}=0$

进一步地，将其展开处理：

$\frac{U_{j}^{n+1}-U_{j}^{n}}{\Delta t}=-\phi _{j}^{n}=L(U_{j}^{n})$

采用三阶RK方法进行时间推进，可分为三个子步：

$U_{j}^{(1)}=U_{j}^{n}+\Delta tL(U_{j}^{n})$

$U_{j}^{(2)}=3/4U_{j}^{n}+1/4[U_{j}^{(1)}+\Delta tL(U_{j}^{(1)})]$

$U_{j}^{n+1}=1/3U_{j}^{n}+2/3[U_{j}^{(2)}+\Delta tL(U_{j}^{(2)})]$

至此，我们已经通过当前时间步 $U_{j}^{n}$ ，获得了下一时间步的 $U_{j}^{n+1}$ ，遍历整个网格，即可得到所有网格节点上的值 $U_{N}^{n+1}$ 。

重复Step2——Step4，直到达到停止时刻。

Wang,Kui

关注

10
点赞
踩
12

收藏

觉得还不错? 一键收藏
3
评论
采用有限差分法求解一维欧拉方程：三阶R-K时间推进——三阶迎风差分（Steger-Warming）

采用有限差分法求解一维欧拉方程，并求解Sod激波管问题，采用三阶R-K时间推进——三阶迎风差分（Steger-Warming）来离散控制方程。
复制链接

扫一扫