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

Wang,Kui

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

阅读量2.6k

点赞数 16

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

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

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

版权

本文通过阅读李新亮老师的代码，探讨了一维差分格式在Sod激波管问题中的应用，重点介绍了三阶迎风格式的构造和通量分裂处理，以及如何利用Steger-Warming方法和Runge-Kutta法进行高精度时间和空间求解，避免了激波问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

#最近看了李新亮老师的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}$