数值方法简介

数值计算方法是将连续的偏微分方程及其边值条件通过一些数值方法将其离散为一个代数方程组,然后利用计算机进行数值求解。常用的数值方法有有限差分法,有限体积法,有限元法等。

有限差分法通过用差商代替微分,将偏微分方程离散为代数方程,这种方法简单有效,容易推导,但通常只针对结构网格。

有限体积法根据守恒定理将方程的积分形式在每一个控制体中列出一个守恒方程,这样可以保证物理量的局部守恒,而且可以适用于任意网格。

有限元方法利用变分原理或加权残值法将偏微分方程变为一种积分形式,这种形式被称为弱形式,然后利用分片多项式对未知量进行插值,带入弱形式中进行求解。有限元方法非常通用和灵活,很容易处理任意边界条件,同时适用于任意网格。

差分

迎风差分格式

$$\frac{\partial \phi}{\partial x} = \frac{\phi_i - \phi_{i-1}}{Δx}$$

表示点$i$只与点$i-1$相关

矩阵形式

$$[\frac{\partial \phi}{\partial x}Δx] = \begin{bmatrix} 1 & 0 & 0 &…&0
\ -1&1 & 0 & …&0
\ 0&-1 &1 & …&0
\ \vdots&\vdots &\vdots&\ddots & \vdots
\ 0&0 &0&-1 & 1\end{bmatrix}\begin{bmatrix} \phi_1 \ \phi_2 \ \phi_3 \ \vdots \ \phi_n \end{bmatrix}$$

矩阵中第一行和最后一行应该由边界条件决定,该矩阵是一个下三角矩阵.

中心差分格式

$$\frac{\partial \phi}{\partial x} = \frac{\phi_{i+1} - \phi_{i-1}}{2Δx}$$
表示点$i$与点$i-1$和点$i+1$相关

矩阵形式

$$[\frac{\partial \phi}{\partial x}Δx] = \frac{1}{2}\begin{bmatrix}
0 & 1 & 0& 0 &…&0
\ -1&0 & 1 & 0&…&0
\ 0&-1 &0 & 1 &…&0
\ \vdots&\vdots &\vdots&\vdots& \ddots& \vdots
\ 0&0&0&-1&0 & 1
\ 0&0 &0&0&-1 & 0
\end{bmatrix}
\begin{bmatrix}
\phi_1 \ \phi_2 \ \phi_3 \ \vdots \ \phi_{n-1} \ \phi_n
\end{bmatrix}$$

矩阵中第一行和最后一行应该由边界条件决定,该矩阵是一个反对称矩阵.

对对流占优方程使用迎风差分格式会产生拖尾效应 (smearing),使用中心差分格式会产生伪震荡现象(wiggles)

smearing_wiggles

有限元

有限元方法在方程两端同时乘以基函数,然后再进行积分,导数项的单元矩阵可以表示为
$$\mathbf{B}_{1x} = \int_{e}\vec{N}(x)\frac{\partial \vec{N}^{T}(x)}{\partial x}dx$$

如果采用线性单元,此时
$$\vec{N}(z) = [λ_1,λ_2]^{T}$$

$$\frac{\partial \vec{N}(z)}{\partial z} = \frac{[-1,1]^{T}}{l}$$

$$\int_{x}dx=l$$

$$\int_{x}λ_idx=\frac{l}{2}$$

$$\int_{x}λ_i^{2}dx=\frac{l}{3}$$

$$\int_{x}λ_iλ_jdx=\frac{l}{6},(i\ne j)$$

所以
$$\mathbf{B}_{1x} =\frac{1}{2}\begin{bmatrix}
-1&1\\
-1&1
\end{bmatrix}$$

将矩阵进行组装之后,得
$$(\varphi,\frac{\partial \phi}{\partial x}) = \frac{1}{2}\begin{bmatrix}
-1 & 1 & 0& 0 &…&0
\ -1&0 & 1 & 0&…&0
\ 0&-1 &0 & 1 &…&0
\ \vdots&\vdots &\vdots&\vdots& \ddots& \vdots
\ 0&0&0&-1&0 & 1
\ 0&0 &0&0&-1 & 1
\end{bmatrix}
\begin{bmatrix}
\phi_1 \ \phi_2 \ \phi_3 \ \vdots \ \phi_{n-1} \ \phi_n
\end{bmatrix}$$

矩阵中第一行和最后一行应该由边界条件决定,除此之外,该矩阵与中心差分是完全一样的,该矩阵是一个反对称矩阵.所以有限元方法是与中心差分类似的方法. 不能直接用于求解双曲方程.

迎风格式探索

格式1

$$\mathbf{B}_{1x} =\begin{bmatrix}
0&0\\
-1&1
\end{bmatrix}$$
那么
$$(\varphi,\frac{\partial \phi}{\partial x}) = \begin{bmatrix} 0 & 0 & 0 &…&0
\ -1&1 & 0 & …&0
\ 0&-1 &1 & …&0
\ \vdots&\vdots &\vdots&\ddots & \vdots
\ 0&0 &0&-1 & 1\end{bmatrix}\begin{bmatrix} \phi_1 \ \phi_2 \ \phi_3 \ \vdots \ \phi_n \end{bmatrix}$$

格式2

$$\mathbf{B}_{1x} =\begin{bmatrix}
1&0\\
-1&0
\end{bmatrix}$$
那么
$$(\varphi,\frac{\partial \phi}{\partial x}) = \begin{bmatrix} 1 & 0 & 0 &…&0
\ -1&1 & 0 & …&0
\ 0&-1 &1 & …&0
\ \vdots&\vdots &\vdots&\ddots & \vdots
\ 0&0 &0&-1 & 0\end{bmatrix}\begin{bmatrix} \phi_1 \ \phi_2 \ \phi_3 \ \vdots \ \phi_n \end{bmatrix}$$

格式3

$$\mathbf{B}_{1x} =\frac{1}{2}\begin{bmatrix}
1&0\\
-2&1
\end{bmatrix}$$
那么
$$(\varphi,\frac{\partial \phi}{\partial x}) = \begin{bmatrix} 1/2 & 0 & 0 &…&0
\ -1&1 & 0 & …&0
\ 0&-1 &1 & …&0
\ \vdots&\vdots &\vdots&\ddots & \vdots
\ 0&0 &0&-1 & 1/2\end{bmatrix}\begin{bmatrix} \phi_1 \ \phi_2 \ \phi_3 \ \vdots \ \phi_n \end{bmatrix}$$

数值方法的性质

数值方法应该满足以下几点性质

  1. 一致性,随着网格尺寸的减小,截断误差应该趋近于0
  2. 稳定性
  3. 收敛性
  4. 物理量守恒
  5. 物理边界