我正在尝试求解二维对流扩散方程，在无量纲变量中可以将我的问题表示为：

$\frac{\partial c}{\partial t} = \lambda^2\frac{\partial^2c}{\partial x^2} + \frac{\partial^2c}{\partial y^2} - 6\, Pe\, y (1 - y)\frac{\partial c}{\partial x}$

这里，$\lambda$是一个几何参数（各种纵横比）和$Pe$是佩克莱特数。这$y\,(1 - y)$由于入口（左壁）处的抛物线速度分布而出现。边界条件为：

• 左壁：狄利克雷型
• 顶壁和右壁：没有助焊剂，所以$y$和$x$梯度分别为零
• 底壁：Robin 型边界条件。基本上是遵循朗缪尔动力学的表面反应。

我已经使用空间的中心差异对方程进行了离散化，并且正在使用 Crank-Nicholson 方案来隐式地推进时间。我的代码崩溃了（c 的值非常大，~$10^{134}$） 为了$Pe > 100$. 它工作正常时$Pe$低于那个。

为什么会这样？我该怎么办？