柱坐标系热传导方程：r=0处奇点与诺依曼边界条件处理咨询

这确实是柱坐标下轴对称热传导数值求解时的经典问题！我来分享几个在实际计算中常用且有效的处理方法，帮你避开r=0处的奇点：

因为你的问题是轴对称的（激光对准中心，与角度无关），r=0处的温度径向梯度必然为0——否则中心位置会出现物理上不可能的温度突变。基于这个对称性，我们可以推导出无奇点的离散格式。

方法1：从微分方程的极限形式推导 #

原方程的径向热传导项可以改写为散度形式：
$$
\frac{\partial^2 T}{\partial r^2} + \frac{1}{r}\frac{\partial T}{\partial r} = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right)
$$
当r→0时，结合$\frac{\partial T}{\partial r}\bigg|{r=0}=0$的对称性条件，对这个散度项取极限：
$$
\lim{r\to0} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right) = 2\frac{\partial^2 T}{\partial r^2}\bigg|{r=0}
$$
因此，r=0处的热传导方程可以替换为：
$$
\frac{\partial T}{\partial t}\bigg|{r=0} = \kappa\left(2\frac{\partial^2 T}{\partial r^2}\bigg|{r=0} + \frac{\partial^2 T}{\partial z^2}\bigg|{r=0}\right) + \frac{1}{\rho c}S(0,z,t)
$$
接下来用显式欧拉法离散：

• 时间导数用前向差分：$\frac{\partial T}{\partial t} \approx \frac{T_{0,z,t+\Delta t} - T_{0,z,t}}{\Delta t}$
• 径向二阶导数用相邻点近似：设r方向步长为$\Delta r$，r=0对应网格点$i=0$，$i=1$对应$r=\Delta r$，则$\frac{\partial^2 T}{\partial r^2}\bigg|{i=0} \approx \frac{T{1,z,t} - 2T_{0,z,t} + T_{-1,z,t}}{(\Delta r)^2}$，但根据对称性$T_{-1,z,t}=T_{1,z,t}$，代入后得到$\frac{2(T_{1,z,t} - T_{0,z,t})}{(\Delta r)^2}$
• z方向二阶导数用常规中心差分即可

把这些代入后，r=0处的显式更新公式就完全没有奇点了。

方法2：直接从离散格式的极限推导 #

对任意内部网格点$i$（$r_i>0$），径向项的离散格式是：
$$
\frac{T_{i+1,z,t} - 2T_{i,z,t} + T_{i-1,z,t}}{(\Delta r)^2} + \frac{1}{r_i}\cdot\frac{T_{i+1,z,t} - T_{i-1,z,t}}{2\Delta r}
$$
当$i=0$时，$r_i=0$，此时第二项是0/0型的不定式，但结合对称性$T_{-1,z,t}=T_{1,z,t}$，分子$T_{i+1}-T_{i-1}=0$，因此第二项的极限为0；第一项代入$T_{-1}=T_1$后变为$\frac{2(T_{1,z,t} - T_{0,z,t})}{(\Delta r)^2}$。最终径向项的离散结果和方法1完全一致，直接用这个结果替换r=0处的径向项即可。

方法3：虚拟网格点（镜像法） #

这是代码实现最友好的方法：在r<0的区域虚拟一个网格点$i=-1$（对应$r=-\Delta r$），根据轴对称性，这个点的温度和$i=1$的点相等，即$T_{-1,z,t}=T_{1,z,t}$。

这样处理后，r=0的点（$i=0$）就可以像普通内部网格点一样使用中心差分格式，计算时直接调用虚拟点的温度值——虚拟点的温度会自动抵消奇点项，不需要额外推导特殊格式。

• 确保显式欧拉法的稳定性：时间步长$\Delta t$必须满足CFL条件，对于二维热传导，大致要求：
  $$
  \Delta t \leq \frac{1}{2\kappa} \left( \frac{1}{(\Delta r)^2} + \frac{1}{(\Delta z)^2} \right)^{-1}
  $$
  具体值需要结合你的材料参数（$\kappa, \rho, c$）和网格步长调整。
• 热源项$S(r,z,t)$在r=0处如果有定义（比如激光聚焦中心），直接代入计算即可，不需要特殊处理。

内容的提问来源于stack exchange，提问作者mjs-wpi