1.3.4边界条件的处理

椭圆型方程边值问题的边界条件差分处理详解

各位同学，今天我们讲解有限差分法处理一般曲边区域椭圆型方程的核心难点——边界条件的离散处理。上一讲我们介绍了正则内点的标准差分格式，但工程中绝大多数问题的求解区域都是曲边的，必然会出现非正则内点：即差分模板的相邻节点落在区域之外，无法直接使用标准五点格式。边界条件处理的本质，就是结合边界上的已知条件，为非正则内点构造合理的差分方程，同时保证格式的精度、稳定性、守恒性等核心性质。

我们以二维Poisson方程 \(-\Delta u = f\) 为模型，分第一类（Dirichlet）、第二/三类（Neumann/Robin）边界条件，系统讲解主流处理方法，分析每种方法的精度、优缺点与适用场景。

---

一、问题背景：非正则内点的定义

对于步长为 \(h_x,h_y\) 的矩形网格，我们定义：

• 非正则内点：属于内部节点集 \(J_\Omega\)，但差分模板（五点格式的相邻节点）中至少有一个节点落在闭区域 \(\overline{\Omega}\) 之外，这类节点的集合记为 \(\tilde{J}_\Omega\)。
• 核心矛盾：非正则内点的标准差分模板需要用到区域外的虚拟节点值，而这些值是未知的，必须通过边界上的已知条件，将虚拟节点值用内部节点和边界节点的已知值表示，或者直接为非正则内点构造不依赖虚拟节点的差分方程。

我们以典型非正则内点 \(P(i,j)\) 为例（如图1-1）：其北邻节点\(N\)、东邻节点\(E\)落在区域外，南邻\(S\)、西邻\(W\)在区域内；边界与网格线的交点为 \(N^*\)（\(PN\)与边界的交点）、\(E^*\)（\(PE\)与边界的交点），\(P^*\)是边界上距\(P\)最近的点。边界上的函数值（Dirichlet边界）或法向导数值（Neumann边界）是已知的，这是构造格式的核心依据。

---

二、第一类（Dirichlet）边界条件的处理方法

Dirichlet边界条件的形式为 \(u|_{\partial\Omega_D} = u_D(x,y)\)，即边界上的函数值已知。我们的目标是为非正则内点\(P\)构造差分方程，核心是利用边界上的已知值\(U_{N^*},U_{E^*},U_{P^*}\)。

方法1：直接插值法

插值法是最直观的方法，核心思想是用边界点和内部节点的函数值，插值得到非正则内点\(P\)的函数值，直接作为该点的离散方程。

（1）零次插值（常数插值）

• 格式：\(U_P = U_{P^*}\)，直接用边界上距\(P\)最近点的函数值作为\(P\)点的取值。
• 截断误差：\(O(h)\)，一阶精度。
• 特点：形式最简单，计算量极小；但精度低，仅适用于网格极密或对精度要求不高的场景，会降低整体格式的收敛阶。

（2）一次线性插值

沿x或y方向做线性插值，利用边界点和同侧内部节点值构造\(U_P\)的表达式：

• x方向插值：设\(h_x^* = |PE^*|\)（\(P\)到边界点\(E^*\)的距离），则
  \[U_P = \frac{h_x U_{E^*} + h_x^* U_W}{h_x + h_x^*} \]

• y方向插值：设\(h_y^* = |PN^*|\)（\(P\)到边界点\(N^*\)的距离），则
  \[U_P = \frac{h_y U_{N^*} + h_y^* U_S}{h_y + h_y^*} \]

• 截断误差：\(O(h^2)\)，二阶精度，与正则内点的五点格式精度匹配。
• 特点：形式简单，精度与内部格式匹配，是工程中最常用的基础方法；缺点是仅用到单方向的两个节点，没有利用Poisson方程的控制方程信息，严格来说是“赋值”而非“差分方程”，可能轻微破坏格式的对称性。

方法2：虚拟节点外推法

外推法的核心思想是：保留标准五点差分格式的形式，引入区域外的虚拟节点\(N,E\)，通过边界条件和泰勒展开，将虚拟节点的函数值用内部节点和边界点的已知值表示，再代入标准五点格式得到\(P\)点的差分方程。

二阶外推的构造

以y方向的虚拟节点\(N\)为例：边界点\(N^*\)在\(PN\)连线上，设\(|PN^*|=h_y^*\)，\(|N^*N|=h_y - h_y^*\)。对\(u\)在\(P\)点做泰勒展开，利用\(u_{N^*}\)、\(u_S\)、\(u_P\)做二阶外推，得到虚拟节点\(N\)的表达式：

\[U_N = \frac{h_y U_{N^*} - (h_y - h_y^*) U_P}{h_y^*} \]

同理，x方向的虚拟节点\(E\)的表达式为：

\[U_E = \frac{h_x U_{E^*} - (h_x - h_x^*) U_P}{h_x^*} \]

将\(U_N,U_E\)代入标准五点差分格式，整理后即可得到\(P\)点的差分方程，截断误差为\(O(h^2)\)，二阶精度。

• 特点：完美保留了标准五点格式的结构，精度与内部格式匹配，格式的对称性、最大值原理都能得到较好的保持；缺点是推导稍复杂，需要处理虚拟节点的外推表达式。

方法3：不等距差分法

不等距差分法的核心思想是：直接对Poisson方程中的二阶偏导数，用非等距节点构造差分逼近，完全不引入虚拟节点，直接得到\(P\)点的差分方程。

格式推导

Poisson方程为 \(-\left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) = f\)，分别对x、y方向的二阶导数做不等距差分：

• x方向：节点为\(W\)（距\(P\)距离\(h_x\)）、\(P\)、\(E^*\)（距\(P\)距离\(h_x^*\)），二阶导数的不等距差分逼近为：
  \[\frac{\partial^2 u}{\partial x^2}\bigg|_P \approx \frac{2}{h_x + h_x^*} \left( \frac{U_{E^*} - U_P}{h_x^*} - \frac{U_P - U_W}{h_x} \right) \]

• y方向：节点为\(S\)（距\(P\)距离\(h_y\)）、\(P\)、\(N^*\)（距\(P\)距离\(h_y^*\)），二阶导数的不等距差分逼近为：
  \[\frac{\partial^2 u}{\partial y^2}\bigg|_P \approx \frac{2}{h_y + h_y^*} \left( \frac{U_{N^*} - U_P}{h_y^*} - \frac{U_P - U_S}{h_y} \right) \]

将两者代入Poisson方程，得到不等距差分格式：

\[-\left[ \frac{2}{h_x + h_x^*} \left( \frac{U_{E^*} - U_P}{h_x^*} - \frac{U_P - U_W}{h_x} \right) + \frac{2}{h_y + h_y^*} \left( \frac{U_{N^*} - U_P}{h_y^*} - \frac{U_P - U_S}{h_y} \right) \right] = f_P \tag{1.3.23} \]

• 截断误差：\(O(h)\)，一阶精度。
• 特点：直接利用控制方程构造，物理意义明确；但格式失去了正则内点五点格式的对称性，最终的线性方程组系数矩阵非对称，无法使用共轭梯度法等针对对称矩阵的快速求解算法，对大规模问题不友好。

对称修正格式

为了保持格式的对称性，可将不等距格式修正为：

\[-\left[ \frac{1}{h_x} \left( \frac{U_{E^*} - U_P}{h_x^*} - \frac{U_P - U_W}{h_x} \right) + \frac{1}{h_y} \left( \frac{U_{N^*} - U_P}{h_y^*} - \frac{U_P - U_S}{h_y} \right) \right] = f_P \tag{1.3.24} \]

• 特点：格式恢复了对称性，系数矩阵对称，可使用快速求解算法；虽然局部截断误差变为\(O(1)\)，但通过最大值原理可证明，整体收敛阶仍然保持\(O(h^2)\)，是工程中兼顾对称性与精度的优选格式。

方法4：有限体积法（守恒型格式）

有限体积法的核心思想是：从Poisson方程的积分守恒形式出发，在非正则内点\(P\)对应的修正控制体\(V_P\)上积分，利用散度定理将体积分转化为边界积分，再用差商逼近法向导数，得到守恒型差分格式。

格式推导

Poisson方程的积分守恒形式为：对任意控制体\(V_P\)，有

\[-\oint_{\partial V_P} \frac{\partial u}{\partial \nu} ds = \int_{V_P} f dxdy \tag{1.3.25} \]

其中\(\nu\)是控制体边界的单位外法向量。

为非正则内点\(P\)构造修正控制体\(V_P\)（如图1-1的虚线区域），控制体的边界一部分是网格线，一部分是区域的边界\(\partial\Omega_D\)。对控制体的每条边，用中点积分公式和差商逼近法向导数，边界上的积分直接用Dirichlet边界条件处理，最终得到守恒型格式：

\[-\left( \frac{U_W - U_P}{h_x} + \frac{U_{E^*} - U_P}{h_x^*} \right) \frac{h_y + \phi h_y^*}{2} - \left( \frac{U_S - U_P}{h_y} + \frac{U_{N^*} - U_P}{h_y^*} \right) \frac{h_x + \theta h_x^*}{2} = f_P \frac{(h_x + \theta h_x^*)(h_y + \phi h_y^*)}{4} \tag{1.3.26} \]

其中\(\theta,\phi\)是控制体边界的修正系数，当\(e,n\)取\(PE^*,PN^*\)的中点时，\(\theta=\phi=1\)，格式与不等距格式(1.3.23)一致。

• 特点：天然满足物理守恒律，这是最大的优势，尤其适用于多介质、间断系数的问题；格式的稳定性、最大值原理都能得到保证；可直接推广到非结构化网格，适配复杂几何区域；缺点是格式推导稍复杂，控制体的构造需要贴合边界。

---

三、第二、三类边界条件的处理方法

第二类（Neumann）边界条件的形式为 \(\frac{\partial u}{\partial \nu}\bigg|_{\partial\Omega_N} = g(x,y)\)，第三类（Robin）边界条件为 \(\frac{\partial u}{\partial \nu} + \alpha u\bigg|_{\partial\Omega_N} = g(x,y)\)，核心难点是边界上法向导数的离散逼近，因为法向通常不与网格线平行，需要结合梯度的投影、法向差分来构造格式。

我们以第二类边界条件为例，第三类边界条件仅需在格式中加入\(\alpha u\)项即可直接推广。

方法1：零阶外推法（梯度投影法）

这是最简单的处理方法，核心思想是：用\(P\)点的梯度在边界外法向上的投影，近似边界点\(P^*\)的法向导数。

设边界点\(P^*\)的单位外法向量\(\nu\)与x轴夹角为\(\alpha\)，则法向导数为：

\[\frac{\partial u}{\partial \nu}\bigg|_{P^*} = \frac{\partial u}{\partial x}\bigg|_P \cos\alpha + \frac{\partial u}{\partial y}\bigg|_P \sin\alpha \]

用向后差商逼近\(P\)点的偏导数：\(\frac{\partial u}{\partial x}\approx \frac{U_P - U_W}{h_x}\)，\(\frac{\partial u}{\partial y}\approx \frac{U_P - U_S}{h_y}\)，代入边界条件\(\frac{\partial u}{\partial \nu}=g(P^*)\)，得到差分方程：

\[\frac{U_P - U_W}{h_x} \cos\alpha + \frac{U_P - U_S}{h_y} \sin\alpha = g(P^*) \tag{1.3.27} \]

• 截断误差：\(O(h)\)，一阶精度。
• 特点：形式最简单，易于实现；但精度低，仅适用于对精度要求不高的场景，法向与网格线夹角较大时，误差会显著增大。

方法2：法向差分法

核心思想是：沿边界的外法线方向构造差分格式，结合线性插值逼近法向导数。

如图1-2，过边界点\(N^*\)的外法线交内部网格线\(PW\)于\(\xi\)点，过\(E^*\)的外法线交\(PS\)于\(\eta\)点。法向导数是沿法线方向的方向导数，因此可以直接构造沿法线的一阶差分：

\[\frac{U_{N^*} - U_\xi}{|\xi N^*|} = g(N^*), \quad \frac{U_{E^*} - U_\eta}{|\eta E^*|} = g(E^*) \tag{1.3.28} \]

其中\(U_\xi\)是\(U_P\)和\(U_W\)的线性插值，\(U_\eta\)是\(U_P\)和\(U_S\)的线性插值。

将插值表达式代入，即可得到关于\(U_P,U_{N^*},U_{E^*},U_W,U_S\)的方程，再结合非正则内点的不等距差分格式，即可封闭方程组。

• 截断误差：\(O(h)\)，一阶精度。
• 特点：直接沿法向构造差分，物理意义明确，误差比梯度投影法更均匀；但需要计算法线与网格线的交点，实现稍复杂。

方法3：有限体积法（守恒型格式）

和Dirichlet边界的处理一致，有限体积法是处理Neumann/Robin边界最自然、最常用的方法，核心优势是能直接将边界的法向导数条件融入积分守恒形式，无需额外处理。

为非正则内点\(P\)构造修正控制体\(V_P\)（如图1-2的虚线区域），对控制体应用积分守恒形式：

\[-\oint_{\partial V_P} \frac{\partial u}{\partial \nu} ds = \int_{V_P} f dxdy \]

控制体的边界分为两部分：内部网格线的边界，和区域边界\(\partial\Omega_N\)的部分。

• 内部网格线的边界：用相邻节点的差商逼近法向导数，和正则内点的有限体积格式一致；
• 区域边界的部分：法向导数\(\frac{\partial u}{\partial \nu}\)直接用边界条件\(g\)代替，积分直接用中点公式近似。

最终得到守恒型格式：

\[-\frac{U_{N_W} - U_P}{|\overline{N_W P}|} |\overline{an_w}| - \frac{U_W - U_P}{h_x} h_y - \frac{U_S - U_P}{h_y} |\overline{s_w b}| - g(P^*) |\tilde{ab}| = f(P) |V_P| \tag{1.3.29} \]

• 截断误差：\(O(h)\)，一阶精度。
• 特点：天然满足积分守恒律，边界条件的处理完全自然融入格式，无需额外构造差分方程；格式的稳定性、收敛性都有严格的理论保证；可直接推广到非结构化网格，是工程中处理复杂边界Neumann/Robin条件的首选方法。

---

四、边界处理的核心原则与方法总结

1. 边界处理的核心优先级

在选择边界处理方法时，格式的整体性质（对称性、守恒性、最大值原理、稳定性）的优先级，高于局部截断误差的阶数，原因如下：

• 局部截断误差仅描述单个节点的逼近精度，而整体性质决定了格式的收敛性、解的物理合理性；
• 很多时候，局部一阶精度的格式，只要保持了守恒性和最大值原理，整体收敛阶仍然能达到二阶；
• 非对称的格式会大幅增加大规模线性方程组的求解成本，甚至导致迭代法不收敛，而对称格式可以使用高效的快速算法。

2. 各类方法核心信息对比表

处理方法	适用边界类型	局部截断误差	格式对称性	守恒性	实现难度	核心优势
零次插值	Dirichlet	\(O(h)\)	无（直接赋值）	不保证	极低	形式最简单，计算量极小
一次线性插值	Dirichlet	\(O(h^2)\)	无（直接赋值）	不保证	低	精度与内部格式匹配，易于实现
虚拟节点外推法	Dirichlet	\(O(h^2)\)	好	不保证	中	保留标准五点格式结构，精度高
不等距差分法	Dirichlet/Neumann	\(O(h)\)	差	不保证	中	直接利用控制方程，物理意义明确
对称修正不等距格式	Dirichlet	局部\(O(1)\)，整体\(O(h^2)\)	好	不保证	中	保持格式对称，适配快速求解算法，整体精度高
有限体积法	所有边界类型	\(O(h)\)	较好	完美保证	中高	天然守恒，适配复杂区域与所有边界类型，稳定性强

3. 不同边界的方法推荐

边界类型	首选方法	备选方法	不推荐场景
第一类（Dirichlet）曲边边界	对称修正不等距格式、有限体积法	一次线性插值、虚拟节点外推法	零次插值（精度过低）
第一类规则矩形边界	直接赋值（边界节点直接取\(u_D\)）	无	所有复杂格式（过度冗余）
第二/三类（Neumann/Robin）边界	有限体积法	法向差分法	零阶外推法（精度低、误差不均匀）

---

五、补充说明

1. 边界处理是有限差分法处理复杂区域的核心难点，也是导致有限差分法在复杂几何区域中应用受限的主要原因。相比之下，有限元法在处理复杂区域的边界条件时，尤其是第二、三类边界条件，更加自然、简单，这也是有限元法在工程中广泛应用的核心原因之一。
2. 对于所有边界处理方法，处理完成后，需要将用到的边界点（如\(N^*,E^*\)）加入边界节点集\(J_D\)，形成扩展的边界节点集，保证线性方程组的未知数与方程数匹配，解存在唯一。
3. 边界条件的处理必须与内部格式的精度、稳定性匹配，避免边界处理的低精度拉低整个格式的收敛阶，同时保证格式的整体稳定性。