24、球形催化剂颗粒与偏微分方程数值解法

关注

原创最新推荐文章于 2026-07-04 14:49:58 发布 63 阅读 AI 写同款· GEO 优化 ›

球形催化剂颗粒与偏微分方程数值解法

1. 球形催化剂颗粒的相关计算

在处理球形催化剂颗粒问题时，我们会涉及到雅可比矩阵元素的计算。对于内部行 $j = 2, 3, \cdots, N - 1$，雅可比矩阵的元素如下：
- $J_{j, j - 1} = A_{j, j - 1}$
- $J_{j, j + 1} = A_{j, j + 1}$
- $J_{j j} = A_{j j} - \xi_{j}^{2}\Phi^{2}g(\phi_{j})$，其中 $g(\phi) = \frac{d}{d\phi}\exp\left(\frac{\gamma\beta(1 - \phi)}{1 + \beta(1 - \phi)}\right)$

在第一行和最后一行，非零元素分别为：
- $J_{11} = A_{11} + a_{1}A_{10} - \xi_{1}^{2}\Phi^{2}g(\phi_{1})$
- $J_{12} = A_{12} + a_{2}A_{10}$
- $J_{N, N - 1} = A_{N, N - 1}$
- $J_{N N} = A_{N N} - \xi_{N}^{2}\Phi^{2}g(\phi_{N})$

由这些元素构成的线性系统是三对角的，即使对于较大的 $N$，消元过程也能快速进行。

接下来是有效因子的定义。通过求解边界值问题（BVP）得到的浓度和温度场，我们可以计算催化剂颗粒内的总反应速率 $R_{tot}$：
$R_{tot} = \int_{0}^{R}k(T)c_{A}(r)(4\pi r^{2})dr$

将其用无量纲量表示为：
$R_{tot} = \left(\frac{4}{3}\pi R^{3}\right)(k_{s}c_{A, s})\left(3\int_{0}^{1}\phi_{A}(\xi)\exp\left(\frac{\gamma\beta(1 - \phi_{A}(\xi))}{1 + \beta(1 - \phi_{A}(\xi))}\right)\xi^{2}d\xi\right)$

其中，第一个方括号内的乘积是催化剂颗粒内无浓度（或温度）梯度时的总反应速率。我们定义有效因子 $\eta$ 为：
$R_{tot} = \left(\frac{4}{3}\pi R^{3}\right)(k_{s}c_{A, s})\eta$

所以，$\eta = 3\int_{0}^{1}\phi_{A}(\xi)\exp\left(\frac{\gamma\beta(1 - \phi_{A}(\xi))}{1 + \beta(1 - \phi_{A}(\xi))}\right)\xi^{2}d\xi$

在 MATLAB 中， catalyst nonisothermal scan.m 用于绘制固定 $\gamma$ 和不同 $\beta$ 值下 $\eta$ 与 $\Phi$ 的关系图。当 $\gamma = 1$，内部 $\Delta\xi = 0.01$，表面附近 $\Delta\xi = 0.001$ 时，结果如图 6.5 所示。

对于零或适度的反应加热（$\beta \approx 0$），增加 $\Phi$ 会降低 $\eta$。这是因为缓慢的扩散会在颗粒中心形成低 $A$ 浓度的耗尽区，导致反应速率降低。相反，如果相对于传导有显著的加热（$\beta \geq 1$），当 $\Phi \approx 1$ 时，$\eta$ 会大于 1，即反应速率比没有内部传输阻力时更高。但在实际中，由于固体传热比传质更好，这种情况很少观察到。因此，当内部传质阻力较强时，$\eta < 1$。

2. 对流/扩散方程的有限差分法

之前我们考虑的是仅存在扩散传输的问题，现在我们来处理引入新数值困难的对流项。考虑简单的偏微分方程：
$\frac{\partial\phi}{\partial t} = -v\frac{\partial\phi}{\partial z} + \kappa\frac{\partial^{2}\phi}{\partial z^{2}} + s(z, t, \phi)$

在稳态且无源项的情况下，该方程变为：
$-v\frac{d\phi}{dz} + \kappa\frac{d^{2}\phi}{dz^{2}} = 0$

加上狄利克雷边界条件：
$\phi(0) = \phi_{0}$
$\phi(L) = \phi_{L}$

其解为：
$\phi(z) = \phi_{0} + \frac{e^{z(Pe)/L} - 1}{e^{(Pe)} - 1}(\phi_{L} - \phi_{0})$

其中，无量纲的佩克莱特数定义为：
$Pe = \frac{vL}{\kappa} = \frac{对流强度}{扩散强度}$

这个解在图 6.6 中展示了不同 $Pe$ 值的情况。当 $Pe \ll 1$ 时，扩散占主导，我们观察到从 $\phi_{0} = 0$ 到 $\phi_{L} = 1$ 的线性增加。当 $Pe \gg 1$ 时，对流占主导，流体从左向右流动，会使 $z = 0$ 处的流入流体向右下游推进，导致在 $z = 1$ 处的出口流出附近 $\phi$ 急剧增加。

2.1 中心差分格式（CDS）

为了数值求解这个问题，我们采用有限差分法，在网格点 $0 < z_{1} < z_{2} < \cdots < z_{N} < L$ 上要求微分方程在每个点局部满足：
$-v\frac{d\phi}{dz}\big| {z {j}} + \kappa\frac{d^{2}\phi}{dz^{2}}\big| {z {j}} = 0$

为了简单起见，我们均匀放置网格点，$\Delta z = \frac{L}{N + 1}$，$z_{j} = j(\Delta z)$。对于二阶导数，我们有近似：
$\frac{d^{2}\phi}{dz^{2}}\big| {z {j}} \approx \frac{\phi_{j - 1} - 2\phi_{j} + \phi_{j + 1}}{(\Delta z)^{2}}$

对于一阶导数，在 CDS 中，我们选择更精确的近似：
$\frac{d\phi}{dz}\big| {(CDS)z {j}} = \frac{\phi_{j + 1} - \phi_{j - 1}}{2(\Delta z)} + O[(\Delta z)^{2}]$

这得到网格点 $j$ 的线性方程：
$-v\left(\frac{\phi_{j + 1} - \phi_{j - 1}}{2(\Delta z)}\right) + \kappa\left(\frac{\phi_{j - 1} - 2\phi_{j} + \phi_{j + 1}}{(\Delta z)^{2}}\right) = 0$

两边乘以 $-\frac{2(\Delta z)^{2}}{\kappa}$ 得到：
$\frac{v(\Delta z)}{\kappa}[\phi_{j + 1} - \phi_{j - 1}] - 2\phi_{j - 1} + 4\phi_{j} - 2\phi_{j + 1} = 0$

我们定义局部佩克莱特数：
$Pe_{loc} \equiv \frac{v(\Delta z)}{\kappa} = (Pe)\left(\frac{\Delta z}{L}\right) = \frac{Pe}{N + 1}$

将 CDS 方程写为：
$-(Pe_{loc} + 2)\phi_{j - 1} + 4\phi_{j} + (Pe_{loc} - 2)\phi_{j + 1} = 0$

令 $\alpha \equiv Pe_{loc} - 2$，$\beta \equiv Pe_{loc} + 2$，CDS 线性系统为：
$\begin{bmatrix}
4 & \alpha & 0 & \cdots & 0 \
-\beta & 4 & \alpha & \cdots & 0 \
0 & -\beta & 4 & \cdots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
0 & 0 & 0 & -\beta & 4
\end{bmatrix}
\begin{bmatrix}
\phi_{1} \
\phi_{2} \
\phi_{3} \
\vdots \
\phi_{N}
\end{bmatrix}
=
\begin{bmatrix}
\beta\phi_{0} \
0 \
0 \
\vdots \
-\alpha\phi_{1}
\end{bmatrix}$

对于所有 $Pe_{loc} \geq 0$，$\beta > 0$，但 $\alpha$ 在 $Pe_{loc} = 2$ 时改变符号。当这种情况发生时，我们预计解会有定性变化。图 6.7 展示了对于 50 个点的网格，不同 $Pe_{loc}$ 值下的 CDS 解。当 $Pe_{loc} > 2$ 时，数值解会出现真实解中不存在的振荡。而且，由于该方程模拟的是密度场的对流和扩散，$\phi$ 取负值在物理上是不合理的，但在数值解中当 $Pe_{loc} > 2$ 时会出现这种情况。为了消除这种虚假振荡，我们可以使网格更细，从而降低固定 $Pe$ 下的 $Pe_{loc}$。但通常我们事先不知道要使网格多小才能避免振荡，而振荡对数值算法收敛的影响可能很严重。

2.2 迎风差分格式（UDS）

我们可以通过对对流项中的一阶导数使用不同的近似来消除振荡，即迎风差分。对于 $v > 0$（从左向右流动），定义为：
$\frac{d\phi}{dz}\big| {(UDS)z {j}} = \frac{\phi_{j} - \phi_{j - 1}}{\Delta z} + O[(\Delta z)]$

这是单边差分，对于有限的 $\Delta z$，这个公式不如中心差分近似准确。网格点 $j$ 的 UDS 方程为：
$-v\left(\frac{\phi_{j} - \phi_{j - 1}}{\Delta z}\right) + \kappa\left(\frac{\phi_{j - 1} - 2\phi_{j} + \phi_{j + 1}}{(\Delta z)^{2}}\right) = 0$

两边乘以 $-\frac{(\Delta z)^{2}}{L}$ 并定义 $Pe_{loc}$ 如前，得到：
$-(Pe_{loc} + 1)\phi_{j - 1} + (2 + Pe_{loc})\phi_{j} - \phi_{j + 1} = 0$

令 $\beta \equiv Pe_{loc} + 2$，$\gamma \equiv Pe_{loc} + 1$，UDS 线性系统为：
$\begin{bmatrix}
\beta & -1 & 0 & \cdots & 0 \
-\gamma & \beta & -1 & \cdots & 0 \
0 & -\gamma & \beta & \cdots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
0 & 0 & 0 & -\gamma & \beta
\end{bmatrix}
\begin{bmatrix}
\phi_{1} \
\phi_{2} \
\phi_{3} \
\vdots \
\phi_{N}
\end{bmatrix}
=
\begin{bmatrix}
\gamma\phi_{0} \
0 \
0 \
\vdots \
\phi_{L}
\end{bmatrix}$

注意，对于所有 $Pe_{loc} \geq 0$，$\beta > 0$，$\gamma > 0$，并且矩阵元素的符号与 $Pe_{loc} < 2$ 时 CDS 系统中的符号相同。如图 6.8 所示，UDS 解对于所有 $Pe_{loc}$ 值都表现良好。

虽然 UDS 方程对对流导数的近似不如 CDS 方程准确，但它们更安全且表现更好。如果使用 CDS 方法时网格点太少，会出现大的振荡，在有非线性源项时可能会很糟糕，但使用 UDS 方法时，$\Delta z$ 太大的影响只是精度损失。

2.3 迎风差分法为何有效

为了理解迎风差分法为何有效，我们通过以下恒等式将迎风和中心差分的一阶导数近似联系起来：
$\frac{\phi_{j} - \phi_{j - 1}}{\Delta z} = \frac{\phi_{j + 1} - \phi_{j - 1}}{2(\Delta z)} + \left(\frac{\phi_{j} - \phi_{j - 1}}{\Delta z} - \frac{\phi_{j + 1} - \phi_{j - 1}}{2(\Delta z)}\right)$

$\frac{d\phi}{dz}\big| {(UDS)z {j}} = \frac{d\phi}{dz}\big| {(CDS)z {j}} + \left(\frac{2\phi_{j} - 2\phi_{j - 1}}{2(\Delta z)} - \frac{\phi_{j + 1} - \phi_{j - 1}}{2(\Delta z)}\right)$

$\frac{d\phi}{dz}\big| {(UDS)z {j}} = \frac{d\phi}{dz}\big| {(CDS)z {j}} + \left(\frac{\Delta z}{2}\right)\left(\frac{-\phi_{j - 1} + 2\phi_{j} - \phi_{j + 1}}{(\Delta z)^{2}}\right)$

认识到方括号内的因子是 $-\frac{d^{2}\phi}{dz^{2}}\big| {z {j}}$ 的近似，我们得到关系：
$\frac{d\phi}{dz}\big| {(UDS)z {j}} = \frac{d\phi}{dz}\big| {(CDS)z {j}} - \left(\frac{\Delta z}{2}\right)\frac{d^{2}\phi}{dz^{2}}\big| {z {j}}$

因此，方程 $v\frac{d\phi}{dz}\big| {(UDS)z {j}} - \kappa\frac{d^{2}\phi}{dz^{2}}\big| {z {j}} = 0$ 的 UDS 离散化等价于带有“额外”扩散的 CDS 离散化：
$v\left(\frac{d\phi}{dz}\big| {(CDS)z {j}} - \left(\frac{\Delta z}{2}\right)\frac{d^{2}\phi}{dz^{2}}\big| {z {j}}\right) - \kappa\frac{d^{2}\phi}{dz^{2}}\big| {z {j}} = 0$

$v\frac{d\phi}{dz}\big| {(CDS)z {j}} - \left(\kappa + \frac{v(\Delta z)}{2}\right)\frac{d^{2}\phi}{dz^{2}}\big| {z {j}} = 0$

迎风差分等价于添加额外的数值扩散或人工扩散，以获得有效扩散常数：
$\kappa_{eff} = \kappa + \frac{v(\Delta z)}{2}$

使得有效局部佩克莱特数始终小于 2：
$Pe_{loc, eff} = \frac{v(\Delta z)}{\kappa_{eff}} = \frac{v(\Delta z)}{\kappa + \frac{v(\Delta z)}{2}} = \frac{2}{1 + \frac{2}{Pe_{loc}}}$

我们也可以通过增加有效扩散常数来避免振荡。在多个空间维度中，通常只在流线方向添加数值扩散，以获得有效的各向异性扩散张量：
$\kappa = \kappa I + \kappa_{num}\frac{vv^{T}}{|v|^{2}}$

3. 最优控制的 HJB 方程的数值解

现在我们来分析最优控制中 HJB 方程的数值解，以最小化成本泛函：
$F\left(u(t); x[0]\right) = \int_{t_{0}}^{t_{H}}\sigma(s, x(s), u(s))ds + \pi(x(t_{H}))$

对于由 $\dot{x} = f(t, x, u)$ 控制的系统。用反向时间 $\tau = t_{H} - t$，我们求解 HJB 方程：
$\frac{\partial\phi}{\partial\tau} = \min_{u(\tau, x)}[\sigma(t_{H} - \tau, x, u) + \nabla\phi \cdot f(t_{H} - \tau, x, u)]$

初始条件为 $\phi(0, x) = \pi(x)$。成本泛函 $F[u(t); x[0]]$ 的最优值是 $\phi(\tau = t_{H} - t_{0}, x[0])$，系统的“最佳”反馈控制律是 $u_{con}(x) = u(\tau = t_{H} - t_{0}, x)$。

我们现在认识到 (6.94) 是一个反应/对流方程，其中 $\sigma$ 充当源项，$f$ 充当 $-v$。因此，在计算的每个点，我们应该检查 $f$ 的局部符号并选择合适的迎风单边差分。在第 5 章的例子中就是这样做的，但没有解释。在实践中，通常还会向 HJB 方程添加人工扩散项 $\epsilon\nabla^{2}\phi$ 以获得“粘性解”。扩散使动态规划更稳健，因为它可以平滑动力学、成本泛函或控制输入中的任何不连续性。

在求解 (6.94) 时，应将空间域扩展到非常大的正和负 $x$ 值，以限制边界条件对感兴趣区域解的影响。要了解如何处理边界，请参考第 5 章的例子。

4. 偏微分方程的特征和类型

从我们对一维对流/扩散方程的讨论中可以看出，方程及其数值解的定性性质会随着参数的变化而改变。这些变化与微分方程在空间和时间中传播场信息的方式有关。信息流动的重要性反映在二阶偏微分方程的常见命名约定中。

考虑二阶微分方程的一般形式：
$a\frac{\partial^{2}\phi}{\partial t^{2}} + b\frac{\partial^{2}\phi}{\partial t\partial z} + c\frac{\partial^{2}\phi}{\partial z^{2}} + f\left(t, z, \phi, \frac{\partial\phi}{\partial t}, \frac{\partial\phi}{\partial z}\right) = 0$

例如，时变一维对流/扩散方程：
$\frac{\partial\phi}{\partial t} - \kappa\frac{\partial^{2}\phi}{\partial z^{2}} + \left(v\frac{\partial\phi}{\partial z}\right) = 0$

在 $Pe \ll 1$ 的极限下，具有以下形式：
$a = 1$
$b = 0$
$c = -\kappa$

此时 $b^{2} - 4ac < 0$，该偏微分方程被称为抛物型。相反，在 $Pe \gg 1$ 的极限下，我们有一个由对流主导的偏微分方程：
$\frac{\partial\phi}{\partial t} + \left(v\frac{\partial\phi}{\partial z}\right) = 0$

对时间求一次导数得到：
$\frac{\partial^{2}\phi}{\partial t^{2}} + v\frac{\partial^{2}\phi}{\partial t\partial z} = 0$

即 $a = 1$，$b = v$，$c = 0$

现在 $b^{2} - 4ac > 0$，方程被称为双曲型。因此，通过改变 $Pe$，我们改变了偏微分方程的类型，这也显著改变了场的传播方式。

稳态扩散方程：
$-\kappa\frac{\partial^{2}\phi}{\partial z^{2}} = 0$

具有二阶偏微分方程的一般形式，其中：
$a = 0$
$b = 0$
$c = -\kappa$

这里 $b^{2} - 4ac = 0$，方程被称为椭圆型。

考虑某个点 $P$，其位置为 $z_{p}$，时间为 $t_{p}$。在时空上，哪些点的场值会影响点 $P$ 的场值，而点 $P$ 的场值又会影响哪些点的场值呢？以一维对流方程为例：
$\frac{\partial\phi}{\partial t} + \left(v\frac{\partial\phi}{\partial z}\right) = 0$

对于 $v > 0$，它描述了 $\phi$ 在 $z$ 增加方向上的纯对流传输。在未来时间 $t > t_{p}$，点 $P$ 只影响下游 $z_{p} < z < z_{p} + v(t - t_{p})$ 的点 $(t, z)$。同样，只有上游 $z_{p} - v(t_{p} - t) < z < z_{p}$ 的过去点 $(t, z)$ 会影响点 $P$。时空图如图 6.9 所示。分隔影响区域和非影响区域的两条线是点 $P$ 的特征线。

我们将特征线与 (6.95) 中的系数 $a$，$b$，$c$ 联系起来。考虑点 $P$ 和其特征线上的点 $Q$，如果 $P$ 和 $Q$ 无限接近：
$z_{q} = z_{p} + dz$
$t_{q} = t_{p} + dt$

我们使用截断的泰勒展开来关联 $P$ 和 $Q$ 处的场值：
$\phi_{q} - \phi_{p} = dt\frac{\partial\phi}{\partial t}\big| {(P)} + dz\frac{\partial\phi}{\partial z}\big| {(P)}$

对一阶导数应用类似的展开：
$\frac{\partial\phi}{\partial t}\big| {(Q)} - \frac{\partial\phi}{\partial t}\big| {(P)} = dt\frac{\partial^{2}\phi}{\partial t^{2}}\big| {(P)} + dz\frac{\partial^{2}\phi}{\partial z\partial t}\big| {(P)}$
$\frac{\partial\phi}{\partial z}\big| {(Q)} - \frac{\partial\phi}{\partial z}\big| {(P)} = dt\frac{\partial^{2}\phi}{\partial t\partial z}\big| {(P)} + dz\frac{\partial^{2}\phi}{\partial z^{2}}\big| {(P)}$

将这两个一阶导数的展开与微分方程在 $P$ 点满足的条件结合起来：
$a\frac{\partial^{2}\phi}{\partial t^{2}}\big| {(P)} + b\frac{\partial^{2}\phi}{\partial t\partial z}\big| {(P)} + c\frac{\partial^{2}\phi}{\partial z^{2}}\big|_{(P)} + f(P) = 0$

得到线性系统：
$\begin{bmatrix}
dt & dz & 0 \
0 & dt & dz \
a & b & c
\end{bmatrix}
\begin{bmatrix}
\frac{\partial^{2}\phi}{\partial t^{2}}\big| {(P)} \
\frac{\partial^{2}\phi}{\partial t\partial z}\big| {(P)} \
\frac{\partial^{2}\phi}{\partial z^{2}}\big| {(P)}
\end{bmatrix}
=
\begin{bmatrix}
\frac{\partial\phi}{\partial t}\big| {(Q)} - \frac{\partial\phi}{\partial t}\big| {(P)} \
\frac{\partial\phi}{\partial z}\big| {(Q)} - \frac{\partial\phi}{\partial z}\big|_{(P)} \
-f(P)
\end{bmatrix}$

如果 $Q$ 在影响区域内部或外部，一阶导数的小变化对应于二阶导数的小变化。但如果 $Q$ 在特征线上，即使一阶导数的无穷小变化也可能对二阶导数产生有限影响，这要求 (6.107) 的矩阵奇异：
$\det\begin{bmatrix}
dt & dz & 0 \
0 & dt & dz \
a & b & c
\end{bmatrix} = 0 = c(dt)^{2} + b(dz)(dt) + a(dz)^{2}$

除以 $(dt)^{2}$，我们得到特征线的斜率：
$\frac{dz}{dt} = \frac{b \pm \sqrt{b^{2} - 4ac}}{2a}$

$b^{2} - 4ac$ 的符号控制着实特征线的数量。对于一维时变对流方程，我们得到：
$\frac{dz}{dt} = \frac{v \pm \sqrt{v^{2} - 4(1)(0)}}{2(1)} = \frac{v \pm v}{2} = v, 0$

因此，我们区分了二阶偏微分方程的三种特征类型：
- 双曲型方程（$b^{2} - 4ac > 0$） ：存在两条实特征线，因此在时空上存在由每个点 $P$ 影响或不影响的不同区域。例如一维对流方程和波动方程 $\frac{\partial^{2}\phi}{\partial t^{2}} - \frac{\partial^{2}\phi}{\partial z^{2}} = 0$。在双曲型方程中，场的信息以行波形式传输，数值误差的影响通常表现为振荡。
- 椭圆型方程（$b^{2} - 4ac < 0$） ：不存在实特征线，无法区分影响区域和非影响区域。在这类问题中，每个点的场值会影响所有其他点的场值。例如二维稳态扩散方程 $\frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\partial^{2}\phi}{\partial y^{2}} = 0$。在椭圆型问题中，场的信息在所有方向上扩散传输，因此数值误差会被“平滑”。
- 抛物型方程（$b^{2} - 4ac = 0$） ：只存在一条实特征线。例如时变一维扩散方程 $\frac{\partial\phi}{\partial t} = \kappa\frac{\partial^{2}\phi}{\partial z^{2}}$。其单条特征线的斜率不确定（因为 $a = 0$），对应于时空图中穿过点 $P$ 的垂直线。因此，每个点 $P$ 会影响所有未来时间的所有空间点，并且会受到过去所有时间所有空间位置的场值影响。随着时间增加，过去的数值误差会被平滑。

5. 具有分散作用的管式化学反应器建模

许多感兴趣的问题涉及多个场，每个场都有自己的控制方程。考虑一个长度为 $L$ 的管式化学反应器（图 6.10），其中发生反应 $A + B \rightarrow C$ 和 $B + C \rightarrow D$。初始时仅包含 $A$ 和 $B$ 的流体以平均轴向速度 $v_{z}$ 流过反应器。除了对流，由于扩散与速度场不均匀性的耦合，还存在类似扩散的混合，称为分散。我们为每个物种使用共同的分散系数 $D$ 作为有效轴向扩散率。分散通常比对流次要，佩克莱特数 $Pe = \frac{v_{z}L}{D}$ 较高。

在稳态下，$A$，$B$，$C$ 和 $D$ 的浓度场由一组耦合的偏微分方程控制：
$\frac{\partial c_{A}}{\partial t} = -v_{z}\frac{\partial c_{A}}{\partial z} + D\frac{\partial^{2}c_{A}}{\partial z^{2}} - k_{1}c_{A}c_{B} = 0$
$\frac{\partial c_{B}}{\partial t} = -v_{z}\frac{\partial c_{B}}{\partial z} + D\frac{\partial^{2}c_{B}}{\partial z^{2}} - k_{1}c_{A}c_{B} - k_{2}c_{B}c_{C} = 0$
$\frac{\partial c_{C}}{\partial t} = -v_{z}\frac{\partial c_{C}}{\partial z} + D\frac{\partial^{2}c_{C}}{\partial z^{2}} + k_{1}c_{A}c_{B} - k_{2}c_{B}c_{C} = 0$
$\frac{\partial c_{D}}{\partial t} = -v_{z}\frac{\partial c_{D}}{\partial z} + D\frac{\partial^{2}c_{D}}{\partial z^{2}} + k_{2}c_{B}c_{C} = 0$

在反应器入口和出口，我们使用 Danckwerts 边界条件。在入口，物种 $j = A$，$B$，$C$，$D$ 进入反应器的通量为 $v_{z}c_{j0}$，但在反应器内部存在分散时，通量为 $v_{z}c_{j} - D\frac{dc_{j}}{dz}$。平衡这两个通量在 $z = 0$ 处得到入口边界条件：
$v_{z}[c_{j}(0) - c_{j0}] - D\frac{dc_{j}}{dz}\big|_{0} = 0$

当 $D \rightarrow 0$ 时，$c_{j}(0) = c_{j0}$，但当 $D \rightarrow \infty$ 时，该边界条件正确地强制 $\frac{dc_{j}}{dz}\big| {0} = 0$。如果反应在流体离开反应器后停止，浓度分布变得均匀，我们使用出口边界条件：
$\frac{dc {j}}{dz}\big|_{L} = 0$

我们使用迎风有限差分法求解这组耦合的偏微分方程。我们放置一个由 $N$ 个均匀间隔的点组成的网格：
$0 < z_{1} < z_{2} < \cdots < z_{N} < L$，$z_{k} = k(\Delta z)$，$\Delta z = \frac{L}{N + 1}$

将每个偏微分方程写为：
$0 = -v_{z}\frac{dc_{j}}{dz} + D\frac{d^{2}c_{j}}{dz^{2}} + s_{j}[{c_{m}(z)}]$

其中 ${c_{m}(z)}$ 表示局部浓度集，每个场的源项为：
- $s_{A}[{c_{m}(z)}] = -k_{1}c_{A}c_{B}$
- $s_{B}[{c_{m}(z)}] = -k_{1}c_{A}c_{B} - k_{2}c_{B}c_{C}$
- $s_{C}[{c_{m}(z)}] = k_{1}c_{A}c_{B} - k_{2}c_{B}c_{C}$
- $s_{D}[{c_{m}(z)}] = k_{2}c_{B}c_{C}$

应用迎风有限差分法，对于每个网格点 $z_{k}$ 和每个物种 $j = A$，$B$，$C$，$D$，我们得到一个非线性代数方程：
$0 = -v_{z}\left(\frac{c_{j}(z_{k}) - c_{j}(z_{k - 1})}{\Delta z}\right) + D\left(\frac{c_{j}(z_{k - 1}) - 2c_{j}(z_{k}) + c_{j}(z_{k + 1})}{(\Delta z)^{2}}\right) + s_{j}[{c_{m}(z_{k})}]$

整理项得到：
$0 = \alpha_{lo}c_{j}(z_{k - 1}) + \alpha_{mid}c_{j}(z_{k}) + \alpha_{hi}c_{j}(z_{k + 1}) + s_{j}[{c_{m}(z_{k})}]$

其中系数为：
- $\alpha_{lo} = \frac{v_{z}}{\Delta z} + \frac{D}{(\Delta z)^{2}}$
- $\alpha_{mid} = -\frac{v_{z}}{\Delta z} - \frac{2D}{(\Delta z)^{2}}$
- $\alpha_{hi} = \frac{D}{(\Delta z)^{2}}$

我们通过对 (6.117) 和 (6.118) 进行离散化来消除不存在的未知数 $c_{j}(z_{0})$ 和 $c_{j}(z_{N + 1})$ 以强制执行边界条件：
$0 = v_{z}[c_{j}(z_{0}) - c_{j0}] - D\left(\frac{c_{j}(z_{1}) - c_{j}(z_{0})}{\Delta z}\right)$
$c_{j}(z_{N + 1}) = c_{j}(z_{N})$

从入口边界条件，我们得到：
$c_{j}(z_{0}) = \frac{Pe_{loc}c_{j0} + c_{j}(z_{1})}{1 + Pe_{loc}}$，其中 $Pe_{loc} = \frac{v_{z}z_{1}}{D}$

因此，对于 $k = 1$，我们有 (6.122) 的修改版本：
$0 = \left(\alpha_{mid} + \frac{\alpha_{lo}}{1 + Pe_{loc}}\right)c_{j}(z_{1}) + \frac{\alpha_{lo}Pe_{loc}c_{j0}}{1 + Pe_{loc}} + \alpha_{hi}c_{j}(z_{2}) + s_{j}[{c_{m}(z_{1})}]$

类似地，对于 $k = N$，我们有 (6.122) 的修改版本：
$0 = \alpha_{lo}c_{j}(z_{N - 1}) + [\alpha_{mid} + \alpha_{hi}]c_{j}(z_{N}) + s_{j}[{c_{m}(z_{N})}]$

我们使用一种标记系统将所有未知数堆叠成一个单一的状态向量，该系统为每个未知数分配一个唯一的整数。方程 (6.122) 形成一组具有稀疏雅可比矩阵的非线性代数方程，为了避免第 1 章中讨论的填充问题，我们应使带宽尽可能小（即，将非零值聚集在主对角线周围）。我们看到 (6.122) 关联了 $c_{j}(z_{k})$，同一场 $c_{j}$ 在相邻点 $c_{j}(z_{k - 1})$ 和 $c_{j}(z_{k + 1})$ 的值，以及所有其他场在同一点 $c_{m \neq j}(z_{k})$ 的值。因此，从两种堆叠状态向量的方式：
- 方案 I：$x^{T} = [c_{A}(z_{1}) c_{A}(z_{2}) \cdots c_{A}(z_{N}) c_{B}(z_{1}) \cdots c_{B}(z_{N}) \cdots c_{D}(z_{N})]$
- 方案 II：$x^{T} = [c_{A}(z_{1}) c_{B}(z_{1}) c_{C}(z_{1}) c_{D}(z_{1}) c_{A}(z_{2}) c_{B}(z_{2}) \cdots c_{C}(z_{N}) c_{D}(z_{N})]$

我们选择方案 II，因为它产生的雅可比矩阵带宽更小（图 6.11）。

我们将 $c_{j}(z_{k})$ 标记为 $n = N_{f}(k - 1) + j$，其中 $N_{f} = 4$ 是场的数量。

在编写涉及多个场的边界值问题的代码时，最好尽可能使用自然名称和索引，例如 $c_{j}(z_{k})$，并依靠例程来堆叠和展开状态向量，以及将矩阵和向量元素放置在适当的位置。这种方法在 tubular reactor 2rxn SS.m 中使用。示例结果如图所示。

综上所述，本文详细介绍了球形催化剂颗粒的相关计算、对流/扩散方程的有限差分法、最优控制的 HJB 方程的数值解、偏微分方程的特征和类型，以及具有分散作用的管式化学反应器建模等内容。通过对这些内容的学习，我们可以更好地理解和处理各种复杂的物理和化学过程中的数学模型。

下面用 mermaid 流程图展示对流/扩散方程有限差分法的求解流程：

graph TD;
    A[开始] --> B[选择方程类型:对流/扩散方程];
    B --> C[确定边界条件];
    C --> D[选择差分格式:CDS或UDS];
    D -- CDS --> E[计算CDS方程系数];
    D -- UDS --> F[计算UDS方程系数];
    E --> G[构建CDS线性系统];
    F --> H[构建UDS线性系统];
    G --> I[求解CDS线性系统];
    H --> J[求解UDS线性系统];
    I --> K[检查解的稳定性:是否有振荡];
    J --> K;
    K -- 有振荡 --> L[细化网格];
    K -- 无振荡 --> M[输出结果];
    L --> D;
    M --> N[结束];

再用表格总结二阶偏微分方程的类型：
| 方程类型 | $b^{2}-4ac$ 的符号 | 特征线数量 | 示例方程 | 信息传播特点 | 数值误差表现 |
| ---- | ---- | ---- | ---- | ---- | ---- |
| 双曲型 | $>0$ | 2 条实特征线 | 一维对流方程、波动方程 | 以行波形式传输 | 通常表现为振荡 |
| 椭圆型 | $<0$ | 无实特征线 | 二维稳态扩散方程 | 各方向扩散传输 | 误差被“平滑” |
| 抛物型 | $=0$ | 1 条实特征线 | 一维时变扩散方程 | 影响未来所有时间空间点 | 过去误差随时间平滑 |

球形催化剂颗粒与偏微分方程数值解法

6. 数值求解过程总结与对比

在前面的内容中，我们介绍了多种数值求解方法，这里对它们进行总结与对比，以便更好地理解和选择合适的方法。

方法	适用方程	优点	缺点	注意事项
中心差分格式（CDS）	对流/扩散方程	一阶导数近似精度较高	当局部佩克莱特数 $Pe_{loc} > 2$ 时，数值解会出现振荡，可能导致物理上不合理的结果；难以事先确定合适的网格大小以避免振荡	需关注 $Pe_{loc}$ 值，必要时细化网格
迎风差分格式（UDS）	对流/扩散方程	解始终表现良好，不会出现振荡；即使网格点较少，也只是精度损失，不会出现灾难性结果	一阶导数近似精度不如 CDS	可适当放宽对网格密度的要求
迎风有限差分法	管式化学反应器的耦合偏微分方程	能有效处理对流项，结合边界条件可求解多场问题；通过合理堆叠未知数可减小雅可比矩阵带宽	涉及多个场和非线性源项，计算相对复杂	注意未知数的标记和堆叠方式，使用自然名称和索引编写代码

7. 实际应用中的考虑因素

在实际应用这些数值方法时，还需要考虑以下因素：

7.1 网格密度

网格密度对数值解的精度和稳定性有重要影响。如在 CDS 方法中，当 $Pe_{loc}$ 较大时，需要更细的网格来避免振荡。但网格过细会增加计算量，因此需要在精度和计算效率之间进行平衡。一般可以先采用较粗的网格进行初步计算，观察解的稳定性和精度，再根据需要逐步细化网格。

7.2 边界条件处理

边界条件的选择和处理直接影响数值解的准确性。不同的问题可能需要不同的边界条件，如在管式化学反应器中使用 Danckwerts 边界条件。在处理边界条件时，要确保其物理意义合理，并且能够正确地离散化到数值计算中。

7.3 非线性源项

当方程中存在非线性源项时，数值求解会变得更加复杂。非线性源项可能导致数值解的不稳定性，如 CDS 方法在非线性源项存在时，振荡问题可能会更加严重。在这种情况下，UDS 方法的优势更加明显，因为它相对更稳健。同时，可能需要采用迭代方法来求解非线性代数方程组。

7.4 计算资源

不同的数值方法对计算资源的需求不同。较复杂的问题和较细的网格会增加计算量，需要更多的内存和计算时间。在实际应用中，要根据可用的计算资源选择合适的方法和网格密度。

8. 未来研究方向

随着科学技术的不断发展，数值求解偏微分方程的方法也在不断进步。以下是一些可能的未来研究方向：

8.1 自适应网格方法

目前我们使用的是均匀网格，但在实际问题中，不同区域的物理特性可能差异很大。自适应网格方法可以根据解的变化情况自动调整网格密度，在解变化剧烈的区域使用更细的网格，在解变化平缓的区域使用较粗的网格，从而在保证精度的同时提高计算效率。

8.2 高精度数值格式

虽然 CDS 和 UDS 是常用的差分格式，但它们的精度还有提升空间。研究更高精度的数值格式可以更准确地模拟物理过程，减少数值误差。例如，一些高阶差分格式或谱方法可以提供更高的精度，但计算复杂度也相应增加。

8.3 多物理场耦合问题

实际问题往往涉及多个物理场的相互作用，如热传导、流体流动、化学反应等。如何有效地处理多物理场耦合问题是一个重要的研究方向。需要开发新的数值方法和算法，以确保不同物理场之间的耦合能够准确模拟。

8.4 并行计算

随着计算机硬件的发展，并行计算成为提高计算效率的重要手段。将数值求解方法并行化，可以充分利用多核处理器和分布式计算资源，大大缩短计算时间。因此，研究适合并行计算的数值算法是未来的一个重要方向。

9. 代码实现示例

以下是一个简单的 MATLAB 代码示例，用于实现对流/扩散方程的中心差分格式（CDS）和迎风差分格式（UDS）的求解：

% 定义参数
L = 1; % 空间长度
N = 50; % 网格点数
dz = L / (N + 1); % 网格间距
v = 1; % 对流速度
kappa = 0.1; % 扩散系数
Pe = v * L / kappa; % 佩克莱特数
Peloc = Pe / (N + 1); % 局部佩克莱特数

% 定义边界条件
phi0 = 0;
phiL = 1;

% 中心差分格式（CDS）
A_CDS = zeros(N, N);
b_CDS = zeros(N, 1);

A_CDS(1, 1) = 4;
A_CDS(1, 2) = Peloc - 2;
b_CDS(1) = (Peloc + 2) * phi0;

for j = 2:N-1
    A_CDS(j, j-1) = -(Peloc + 2);
    A_CDS(j, j) = 4;
    A_CDS(j, j+1) = Peloc - 2;
    b_CDS(j) = 0;
end

A_CDS(N, N-1) = -(Peloc + 2);
A_CDS(N, N) = 4;
b_CDS(N) = -(Peloc - 2) * phiL;

phi_CDS = A_CDS \ b_CDS;

% 迎风差分格式（UDS）
A_UDS = zeros(N, N);
b_UDS = zeros(N, 1);

A_UDS(1, 1) = Peloc + 2;
A_UDS(1, 2) = -1;
b_UDS(1) = (Peloc + 1) * phi0;

for j = 2:N-1
    A_UDS(j, j-1) = -(Peloc + 1);
    A_UDS(j, j) = Peloc + 2;
    A_UDS(j, j+1) = -1;
    b_UDS(j) = 0;
end

A_UDS(N, N-1) = -(Peloc + 1);
A_UDS(N, N) = Peloc + 2;
b_UDS(N) = phiL;

phi_UDS = A_UDS \ b_UDS;

% 绘制结果
z = (1:N) * dz;
figure;
plot(z, phi_CDS, 'b-o', 'DisplayName', 'CDS');
hold on;
plot(z, phi_UDS, 'r-s', 'DisplayName', 'UDS');
xlabel('z');
ylabel('\phi');
title('Convection/Diffusion Equation Solution');
legend;
grid on;

这个代码示例展示了如何使用 CDS 和 UDS 方法求解对流/扩散方程，并绘制了两种方法的解。通过运行这个代码，你可以直观地观察到两种方法的差异。

10. 总结与展望

本文详细介绍了球形催化剂颗粒的相关计算、对流/扩散方程的有限差分法、最优控制的 HJB 方程的数值解、偏微分方程的特征和类型，以及具有分散作用的管式化学反应器建模等内容。我们了解了不同数值方法的优缺点和适用范围，以及在实际应用中需要考虑的因素。

未来，随着数值计算技术的不断发展，我们可以期待更高效、更精确的数值方法的出现，以解决更复杂的物理和化学问题。同时，结合实验数据和数值模拟结果，我们可以更好地理解和优化各种物理和化学过程，为实际应用提供更有力的支持。

下面用 mermaid 流程图展示实际应用中选择数值方法的流程：

graph TD;
    A[开始] --> B[确定问题类型:对流/扩散、多场耦合等];
    B --> C[评估对流和扩散强度:计算Pe值];
    C -- Pe较小 --> D[可优先考虑CDS];
    C -- Pe较大 --> E[考虑UDS或其他稳健方法];
    D --> F[检查是否有非线性源项];
    E --> F;
    F -- 有非线性源项 --> G[UDS更合适];
    F -- 无非线性源项 --> H[根据精度和计算资源选择];
    G --> I[确定网格密度和边界条件];
    H --> I;
    I --> J[编写代码实现数值求解];
    J --> K[检查解的合理性和稳定性];
    K -- 不合理或不稳定 --> L[调整参数或方法];
    K -- 合理稳定 --> M[输出结果并分析];
    L --> I;
    M --> N[结束];

希望本文能够帮助读者更好地理解和应用这些数值方法，解决实际问题。如果你有任何疑问或建议，欢迎留言讨论。