有限差分法（Finite Difference Method）解方程：边界和内部结点的控制方程

最新推荐文章于 2022-12-22 12:15:42 发布

CFD_Tyro

最新推荐文章于 2022-12-22 12:15:42 发布

阅读量3.7k

点赞数 4

分类专栏：理论基础文章标签：线性代数算法矩阵

原文链接：https://www.jianshu.com/p/39f70fdce192

版权

理论基础专栏收录该内容

11 篇文章 5 订阅

订阅专栏

本文转载自本人发布在简书平台上的文章

1. FDM解常微分方程

1.1. 问题描述

$\frac{d^2\phi}{dx^2}=S_{\phi} \tag{1}$
这是二阶常微分方程（second-order Ordinary Differential Equation, ODE），考虑最简单的情况即 $S = 0$ ，积分后可得 $\phi=c_1x+c_2$ ，有两个待定系数，因此要求解该方程必须提供两个边界条件（因为方程中不包含时间项，因此无初始条件），例如

$\phi(x_L)=\phi_L \quad \phi(x_R)=\phi_R$

1.2. 场和边界

我们把上述方程转化为一个场（field）。试想存在一维（因为 $\phi$ 仅与一个变量有关）直线，线上有若干等距分布的结点（node），每个结点都有唯一的 $\phi$ ，那么 $\phi$ 和 $x$ 的关系满足上述方程。

内部结点满足方程，那么边界上的 $\phi(x_L)$ 和 $\phi(x_R)$ 呢？事实上，边界应是内部结点的延申，即边界点也应该满足上述方程。这也是为什么我们可以通过两个边界条件求解 $\phi=c_1x+c_2$ 的原因，这隐含的假设就是边界点满足常微分方程。

1.3. 结点上的值

用数值方法求解上述方程等价于求解场内每个结点上的 $\phi$ ，结点 $i$ 的上式表达为

$\frac{d^2\phi}{dx^2}\bigg|_{i}$

如何从结点值 $\phi_i$ 得到导数值？很自然想到用Taylor展开。需要注意的是，Taylor展开隐含的假设是 $\phi$ 无限可导。将 $i$ 关于周围两点做Taylor展开，即

假设等距，则

$\frac{d^2\phi}{dx^2}\bigg|_{i} =\frac{\phi_{i+1}+\phi_{i-1}-2\phi_i}{(\Delta x)^2}-\frac{(\Delta x)^2}{12} \frac{d^4\phi}{dx^4}\bigg|_{i} + ... \tag{2}$

$\varepsilon_i=-\frac{(\Delta x)^2}{12} \frac{d^4\phi}{dx^4}\bigg|_{i} + ...$

其中 $\varepsilon_i$ 为离散误差或截断误差（discretization or truncation error），该误差正比于距离的平方，因此我们称上式对导数的逼近有二阶精度。

1.4. 边界条件

上一部分是对场内部结点的推导。边界上结点的值为边界条件，有以下三种形式：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JFA8eir9-1659235581146)(https://upload-images.jianshu.io/upload_images/27641975-f329c2341566c1ba.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)]

1.4.1. Dirichlet

$\phi_{i=1}=\phi_L$

1.4.2. Neumann

$\frac{d\phi}{dx}\bigg|_{i=1}=J_L \tag{3}$

将结点2关于结点1做Taylor展开，可以得到一阶精度的梯度表达式

$\frac{d\phi}{dx}\bigg|_{i=1}=\frac{\phi_2-\phi_1}{\Delta x}+\varepsilon(\Delta x) \tag{4}$

将结点3关于结点1做Taylor展开，结合上式可以进一步得到二阶精度

$\frac{d\phi}{dx}\bigg|_{i=1}=\frac{4\phi_2-\phi_3-3\phi_1}{2\Delta x}+\varepsilon(\Delta x^2) \tag{5}$

对于非Dirichlet的边界条件，应尽量让边界点也满足内部结点的控制方程。
对 $\phi_2$ 和 $\phi_3$ 关于结点1展开可得

$\phi_2=\phi_1+(\Delta x)\frac{d\phi}{dx}\bigg|_{i=1}+(\Delta x^2)\frac{d^2\phi}{dx^2}\bigg|_{i=1}... \\ \phi_3=\phi_1+(2\Delta x)\frac{d\phi}{dx}\bigg|_{i=1}+(4\Delta x^2)\frac{d^2\phi}{dx^2}\bigg|_{i=1}...$

将边界条件公式(3)代入上式可得

$\phi_2=\phi_1+(\Delta x)J_L+... \\ \phi_3=\phi_1+(2\Delta x)J_L+...$

将上式整理并消去三阶导数项（ $d^3\phi/dx^3$ ）可得结点1满足的二阶导数项的离散格式：

$\frac{d^2\phi}{dx^2}\bigg|_{i=1}=\frac{8\phi_2-\phi_3-7\phi_1-(6\Delta x)J_L}{2(\Delta x)^2}+\varepsilon(\Delta x^2) \tag{6}$

上式为结点1满足的控制方程，并用上了梯度边界条件。需要说明的是，不用公式(6)而用公式(4,5)同样可以求解问题，但是前者的精度更高。

上述处理方式的意义：FDM的操作流程是对微分方程逐项离散、在每个结点上离散。本节离散的ODE只包含一个二阶导数项，需要在每个结点上对其离散。前文对内部结点采用二阶中心差分离散，我们同样希望在边界点对二阶导数项做离散，因此有了方程（6）.

除了上述采用内结点对边界结点Taylor展开的方式，《数值传热学》4.3.2小节给出了另一种方法，即设置边界外的辅助结点，此处不再展开。

1.4.3. Robin

$\alpha \phi_1 + \beta \frac{d\phi}{dx}\bigg|_{i=1}=\gamma$

类似地，令边界结点1既满足控制方程又满足边界条件。对上式移项可得

$\frac{d\phi}{dx}\bigg|_{i=1}=\frac{\gamma-\alpha \phi_1}{\beta}$

同样代入 $\phi_2$ 和 $\phi_3$ 的Taylor展开，

$\phi_2=\phi_1+(\Delta x) \left( \frac{\gamma-\alpha \phi_1}{\beta} \right) +... \\ \phi_3=\phi_1+(2\Delta x) \left( \frac{\gamma-\alpha \phi_1}{\beta} \right) +...$

消去三阶导数项并移项可得

$\frac{d^2\phi}{dx^2}\bigg|_{i=1}=\frac{8\phi_2-\phi_3-7\phi_1-(6\Delta x)J_L}{2(\Delta x)^2}+\varepsilon(\Delta x^2) \tag{7}$

其中

$J_L=\frac{d\phi}{dx}\bigg|_{i=1}=\frac{\gamma-\alpha \phi_1}{\beta}$

1.5. 组建矩阵

有了内部结点的离散格式和边界条件，我们便可对每个结点列方程。由于结点 $i$ 离散格式中势必会包含其他结点信息，例如对结点 $i$

$A_{i,1}\phi_1 + ... + A_{i,i-1}\phi_{i-1} + A_{i,i}\phi_i + A_{i,i+1}\phi_{i+1} + ... + A_{i,N}\phi_N =Q_i$

将所有结点的方程组合起来，借助矩阵运算求解。

2. FDM求解泊松方程：二维问题

二维Poisson方程：

$\nabla^2 \phi=\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = S_{\phi} \tag{8}$

如果 $S_{\phi}=0$ ，则退化为Laplace方程。二阶偏导数项同样使用Taylor展开，只不过针对该问题应使用二维展开。如果使用正交网格结点，则相当于在每一维独自做Taylor展开。

按下图所示，网格结点可分为内部结点和边界结点。

2.1. 内部结点控制方程

由上文公式(2)可知，两个二阶偏导数项可分别写为

$\frac{\partial^2\phi}{\partial x^2}\bigg|_{i,j} =\frac{\phi_{i+1,j}+\phi_{i-1,j}-2\phi_{i,j}}{(\Delta x)^2} + \varepsilon(\Delta x^2) \tag{9}$

$\frac{\partial^2\phi}{\partial y^2}\bigg|_{i,j} =\frac{\phi_{i,j+1}+\phi_{i,j-1}-2\phi_{i,j}}{(\Delta y)^2} + \varepsilon(\Delta y^2) \tag{10}$

去掉高阶项，可得内部结点控制方程

$\frac{\phi_{i+1,j}+\phi_{i-1,j}-2\phi_{i,j}}{(\Delta x)^2} + \frac{\phi_{i,j+1}+\phi_{i,j-1}-2\phi_{i,j}}{(\Delta y)^2} \approx S_{i,j} \tag{11}$

2.2. 下边界

Robin边界条件

$\alpha \phi (x,0) + \beta \frac{\partial \phi}{\partial y}\bigg|_{x,0} = \gamma$

对于下边界的某个非边角结点 $(i, 1)$ ，其中 $i\neq 1, i\neq N$ ，根据公式(7)来构建既满足边界条件又满足内部结点控制方程的离散格式，即公式(7)

$\frac{\partial^2\phi}{\partial y^2}\bigg|_{i,1} = \frac{8\phi_{i,2}-\phi_{i,3}-(7-6\Delta y \frac{\alpha}{\beta}) \phi_{i,1} - 6\Delta y (\gamma / \beta)}{2(\Delta y)^2} \tag{12}$