pdeSolver 工业 PDE 计算教程

长为 1 的均匀棒材，端点维持零温度基准，内部持续均匀加热。归一化方程为 $u_t=\alpha u_{xx}+2\alpha$ ，初值为

$u(x,0)=0$ ，端点为 $u(0,t)=u(1,t)=0$ 。

独立解析基准。在稳态 $u_{xx}=-2$ ，故 $u_{*}(x)=x(1-x)$ ，中心温度恰为 0.25。这是解析值，不是 pdeSolver 已运行输出。

设置为 pdeCoeff1D(diff=alpha, source=2*alpha)，boundary=0.0，initial=0.0。固定终止时间和空间节点计算 $L_{\infty }$ 误差，再分别加密 $N_x=51/101/201$ 与 $N_t=100/200/400$ ，同时检查非负性和对称性。

完整代码见第 3.3 节“快速跑通”。本例只是将源项、边界和解析参考写成一个最小的一维热传导问题。

运行结果：N_x=101、N_t=400 时，maxError=1.107879e-09，symmetryError=1.998401e-15。空间和时间收敛结果见第 7 章。

通道内存在正向流速和内部源项，方程为 $c_t=D{c}_{xx}-v{c}_x-kc+s(x)$ 。在 pdeCoeff1D 中设置 conv=-v，与物理速度正方向相反。

总流量（单位截面）为 $J=vc-D{c}_x$ ，所以 $c_x=0$ 只表示零扩散梯度，不表示零总通量。本例用入口 $c(0,t)=0$ 、出口 $c_x(1,t)=0$ 近似开放通道。若边界封闭，应直接建模 $vc-D{c}_x=0$ ，再按 API 的外法向符号转换为 Robin 条件。

检查入口、出口、非负性、浓度峰和质心位移，并与 v=0 的结果对照。高 Peclet 数时优先试 upwind，但要通过网格加密评估其数值扩散。

def pollutantSource(x, t) {
    normalizedDistance = (x - 0.2) / 0.03
    return exp(-normalizedDistance * normalizedDistance)
}

xGrid = (0..100) / 100.0
domain = pdeDomain1D(xGrid)

outletZeroGradient = dict(STRING, ANY)
outletZeroGradient["type"] = "neumann"
outletZeroGradient["value"] = 0.0
boundary = [0.0, outletZeroGradient]  // xMin: 固定入口浓度; xMax: 零梯度近似出流

coeff = pdeCoeff1D(
    diff=0.003,
    conv=-0.3,
    react=-0.2,
    source=pollutantSource
)

method = dict(STRING, ANY)
method["convectionScheme"] = "upwind"

setting = dict(STRING, ANY)
setting["resultMode"] = "detail"

result = pdeSolver(
    coeff,
    0.0,
    domain,
    (0..200) * (0.5 / 200.0),
    boundary,
    method,
    setting
)

运行结果：最终场浓度范围约为 [0.0, 0.154389]。该数值对应当前 upwind、网格和时间步；改变 Peclet 数或网格后应重新检查峰值和非负性。

考虑封闭一维介质中的局部高斯脉冲，方程为 $u_t=D{u}_{xx}$ 。端点没有对流，零法向导数同时表示零扩散通量和零总通量。

数学性质。积分质量 $M(t)={\int }_0^{1}u(x,t)dx$ 常数；峰值不增加，空间分布逐渐平缓，质量归一化方差趋向均匀分布的极限。

用与网格匹配的积分权重计算 $|M(t)-M(0)|/M(0)$ 。非均匀网格不能直接用节点值求和代替积分。

def gaussianPulse(x) {
    normalizedDistance = (x - 0.5) / 0.04
    return exp(-0.5 * normalizedDistance * normalizedDistance)
}

xGrid = (0..200) / 200.0
tGrid = (0..100) * (0.2 / 100.0)

zeroFlux = dict(STRING, ANY)
zeroFlux["type"] = "neumann"
zeroFlux["value"] = 0.0

setting = dict(STRING, ANY)
setting["resultMode"] = "detail"
setting["output"] = "all"

result = pdeSolver(
    pdeCoeff1D(diff=0.002),
    gaussianPulse,
    pdeDomain1D(xGrid),
    tGrid,
    [zeroFlux, zeroFlux],
    NULL,
    setting
)

运行结果：output="all" 返回 201×101=20301 个网格—时间值。均匀网格梯形积分得到初始质量 0.100265、末时刻质量 0.100265，相对漂移 3.737095e-15。

一根长度为 $L=0.20\mathrm{m}$ 的细长实心钢棒，横截面足够小，可采用一维轴向温度场 $T(x,t)$ 近似。棒材取常数热物性，体内存在均匀体积热源，两端与环境发生对流换热。忽略辐射、横向温度梯度、材料相变及热物性随温度变化；这只是可复现的教学工程模型，不可未经校准用作具体设备设计结果。

这里使用的是教学参数，不对应某种特定钢材或设备测量值。真实棒材若存在横向散热，应把侧壁散热折算为体积源/汇项，并使用实际横截面周长和面积。

在常数物性和无轴向流动条件下：

${\rho }_m{c}_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial x^2}+Q,x\in (0,L).$

其中 ${\rho }_m{c}_p$ 的单位为 J/(m³·K)， $k{T}_{xx}$ 和 Q 的单位均为 W/m³。直接以物理量作为 PDE 系数：mass=rho*cp，diff=k，source=Q，无需预先归一化。温度可以使用 °C，因为该模型只涉及温差和温度导数。

两端向环境散热： $-k{\partial }_{n}T=h(T-T_{\infty })$ ，即

$\underbrace{\frac{h}{k}}\alpha T+\underbrace{1}\beta {\partial }_{n}T=\underset{\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}}{\underbrace{\frac{h}{k}{T}_{\infty }}}.$

注意左端外法向为 -x，右端外法向为 +x；两端使用相同 Robin 字典的原因是 pdeSolver 的导数以外法向定义，而不是两端统一使用 $+{\partial }_x$ 。初值 $T(x,0)=T_{\infty }$ 在 t=0 满足两端边界。

// 参数均采用 SI 单位：m, s, kg, J, W, K（温差）；显示温度为摄氏度。
length = 0.20
rhoMaterial = 7850.0
heatCapacity = 470.0
conductivity = 45.0
volumetricHeatSource = 12000.0
heatTransferCoefficient = 25.0
ambientTemp = 20.0
simulationSeconds = 86400.0

xGrid = (0..120)*(length/120.0)
tGrid = (0..480)*(simulationSeconds/480.0)

// 原始量形式：rho*c_p*T_t = k*T_xx + Q
coeff = pdeCoeff1D(
    diff=conductivity,
    source=volumetricHeatSource,
    mass=rhoMaterial*heatCapacity
)

// Robin alpha*T + beta*dT/dn = value，外法向 n 随端点变化。
convectiveEnd = dict(STRING, ANY)
convectiveEnd["type"] = "robin"
convectiveEnd["alpha"] = heatTransferCoefficient/conductivity
convectiveEnd["beta"] = 1.0
convectiveEnd["value"] = (heatTransferCoefficient/conductivity)*ambientTemp
boundary = [convectiveEnd, convectiveEnd]

method = dict(STRING, ANY)
method["scheme"] = "implicit"
setting = dict(STRING, ANY)
setting["resultMode"] = "detail"
setting["interpAt"] = length/2.0
setting["diagnostics"] = true

result = pdeSolver(
    coeff, ambientTemp, pdeDomain1D(xGrid), tGrid,
    boundary, method, setting
)

// 解析稳态：仅作为长时极限基准，不应冒充 t_f 时刻的精确解。
steadyTemp = ambientTemp + volumetricHeatSource*length/(2.0*heatTransferCoefficient) + volumetricHeatSource*xGrid*(length-xGrid)/(2.0*conductivity)
steadyCenter = ambientTemp + volumetricHeatSource*length/(2.0*heatTransferCoefficient) + volumetricHeatSource*length*length/(8.0*conductivity)

summary = dict(STRING, ANY)
summary["tempCenterAtFinalTime"] = result["interpValue"]
summary["steadyCenterReference"] = steadyCenter
summary["maxDistanceToSteady"] = max(abs(result["value"]-steadyTemp))
summary

长时间后 $T_t=0$ ，因此 $k{T}^{″}+Q=0$ ，结合两端对流 Robin 边界得到：

$T_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{y}}(x)=T_{\infty }+\frac{QL}{2h}+\frac{Qx(L-x)}{2k}$

示例参数下，稳态两端温度为 $68{}^{\circ }\mathrm{C}$ ，中心温度为 $69.333333{}^{\circ }\mathrm{C}$ 。

独立 Python 有限体积（隐式 Euler）瞬态参考，t=86400 s：

这些是独立 FVM 数值结果，不是 DolphinDB 运行结果，也不能代替时间趋于无穷时的解析稳态。实际 pdeSolver 精度验收应与同时间、同物性、同边界的参考瞬态比较，而非要求 24 小时结果严格等于稳态。

这两个数字是公式计算的解析稳态基准，不是 $t=86400s$ 时 DolphinDB 的求解结果。瞬态结果应随时间趋近稳态，初始时刻严格为 $20{}^{\circ }\mathrm{C}$ 。

• 单位：PDE 各项为 W/m³，Robin 左右两项为 K/m，热通量还要乘以 k。

• 网格：至少比较 $N_x=61,121,241$ 和对应时间步的中心温度、端点温度及稳态距离。

• 参考：validation/reference_regression.py 使用独立 FVM 隐式推进；真实设备还需校准热物性、侧壁换热和不确定性。

验证结果：解析关系、独立 FVM 参考、DolphinDB 运行结果以及空间/时间收敛数据均已记录。更换部署版本后，应在相同时间点重新做 FVM 对照，并复核生产边界条件。

在 (0.3,0.7) 附近持续注入窄高斯热源，方程为 $u_t=D_x{u}_{xx}+D_y{u}_{yy}+s(x,y)$ 。本例取 $D_x=D_y=0.002$ ，四边定温为 0。

方法。构造 41×41 均匀网格和 41×41 双方向 sinh 聚焦网格，再以更细网格参考比较热点中心温度。局部加密提高关注区域精度，但不一定降低全域最大误差。

比较时保持两组网格的总节点数接近，记录聚焦区域的最小步长，并与细网格结果比较。时间步也要单独加密，避免把时间误差误认为网格收益。聚焦网格是否更准，以运行结果为准。

def hotspotSource(x, y, t) {
    scaledX = (x - 0.3) / 0.03
    scaledY = (y - 0.7) / 0.03
    return exp(-scaledX * scaledX - scaledY * scaledY)
}

finalTime = 0.05
tGrid = (0..100) * (finalTime / 100.0)

uniformGrid = (0..40) / 40.0
focusedX = pdeGrid1D(0.0, 1.0, 41, "sinh", 0.3, 3.0)
focusedY = pdeGrid1D(0.0, 1.0, 41, "sinh", 0.7, 3.0)
referenceGrid = (0..120) / 120.0

coeff = pdeCoeff2D(
    xx=0.002,
    yy=0.002,
    source=hotspotSource
)

setting = dict(STRING, ANY)
setting["resultMode"] = "detail"
setting["interpAt"] = 0.3 0.7

method = dict(STRING, ANY)
method["parallel"] = false

uniform = pdeSolver(
    coeff, 0.0,
    pdeDomain2D(uniformGrid, uniformGrid),
    tGrid, 0.0, method, setting
)

focused = pdeSolver(
    coeff, 0.0,
    pdeDomain2D(focusedX, focusedY),
    tGrid, 0.0, method, setting
)

reference = pdeSolver(
    coeff, 0.0,
    pdeDomain2D(referenceGrid, referenceGrid),
    tGrid, 0.0, method, setting
)

运行结果：均匀网格、聚焦网格和参考网格的插值值分别为 0.043145、0.041562 和 0.041574。当前结果显示聚焦网格更接近参考值，正式网格仍应按同样方法复核。

考虑局部污染源和二维流场，方程为 $c_t=D_x{c}_{xx}+D_y{c}_{yy}-v_x{c}_x-v_y{c}_y+s(x,y)$ 。源点为 (0.2,0.5)，比较 $(v_x,v_y)=(1,0)$ 和 (1, 0.4)。

通量为 $\mathbf{J}=\mathbf{v}c-D\nabla c$ 。有穿过边界的对流时， ${\partial }_{n}c=0$ 不等于 $\mathbf{J}·\mathbf{n}=0$ 。本例使用固定零浓度入口和零梯度出口，适用性取决于流向和截断域；真正的无流出壁面应满足 $\mathbf{v}·\mathbf{n}=0$ 或相应 Robin 通量条件。

把等效网格积分质量作为监测量，用浓度质心 $\overline{x},\overline{y}$ 判断输运方向。 $v_y>0$ 时要检查上边界是否确实起到出口作用，并通过增加网格密度和域高度评估截断边界的影响。

def plumeSource(x, y, t) {
    scaledX = (x - 0.2) / 0.035
    scaledY = (y - 0.5) / 0.035
    return exp(-scaledX * scaledX - scaledY * scaledY)
}

xGrid = (0..60) / 60.0
yGrid = (0..40) / 40.0
tGrid = (0..120) * (0.3 / 120.0)
domain = pdeDomain2D(xGrid, yGrid)

outletZeroGradient = dict(STRING, ANY)
outletZeroGradient["type"] = "neumann"
outletZeroGradient["value"] = 0.0
boundary = [0.0, outletZeroGradient, 0.0, outletZeroGradient]  // 入口固定 0，出口零梯度近似

method = dict(STRING, ANY)
method["parallel"] = false
method["convectionScheme"] = "upwind"

reportSetting = dict(STRING, ANY)
reportSetting["resultMode"] = "detail"
reportSetting["diagnostics"] = true

axisAligned = pdeSolver(
    pdeCoeff2D(
        xx=0.001, yy=0.001,
        x=-1.0, y=0.0,
        source=plumeSource
    ),
    0.0, domain, tGrid,
    boundary, method,
    reportSetting
)

oblique = pdeSolver(
    pdeCoeff2D(
        xx=0.001, yy=0.001,
        x=-1.0, y=-0.4,
        source=plumeSource
    ),
    0.0, domain, tGrid,
    boundary, method,
    reportSetting
)

在 3.00.6.2 上，两组代码均执行成功，diagnostics.success=true 且无告警。轴向流和斜向流的最终最大浓度分别为 0.056756 和 0.039527，最小值接近 0。边界适用性仍应结合流向和网格加密判断。

以下数据来自目标版本 3.00.6.2 Server，误差取终止时间的最大绝对误差或指定插值点误差。对应的批量脚本是 validation/convergence_scans.dos；二维制造解扫描使用 validation/mms_convergence.dos。换用生产网格后，按同样方式重新扫描。

一维解析热方程。空间扫描固定 N_t=800、theta 格式：

时间扫描固定 N_x=401、隐式 Euler：

MMS 空间扫描由 validation/mms_convergence.dos 复现。固定 N_t=600 时，N_x=41/81/121/161 的误差为 7.539e-5/1.886e-5/8.382e-6/4.715e-6。validation/convergence_scans.dos 的二维 MCS 空间扫描固定 N_t=200，N=21/41/61/81 的误差为 3.330e-4/8.327e-5/3.701e-5/2.082e-5；二维时间扫描固定 N=81，N_t=50/100/200/400/800 的误差为 2.075e-5/2.080e-5/2.082e-5/2.082e-5/2.082e-5。时间步继续加密收益很小，当前配置主要受空间误差限制。

专项入口为 ZIP 中的 03_2D_Industrial/IND-031/example.dos。模型为

$u_t=0.02(u_{xx}+u_{yy}),(x,y)\in [0,1{]}^2$

初值为 $1+0.2\cos (\pi x)\cos (\pi y)$ ，四边采用零 Neumann 边界，空间网格为 41×41，终止时间为 0.5。解析解用于检查场值误差；二维梯形权重用于检查

$M(t)={\int }_0^1{\int }_0^{1}u(x,y,t)dxdy$

目标生产版本 3.00.6.2 2026.09.18 LINUX x86_64 的运行结果如下：

本专项在当前模型、均匀网格和边界条件下通过；它说明该配置的离散质量误差接近机器精度，但不能外推为所有对流、源项、反应或 Robin 边界问题都严格守恒。非均匀网格应使用对应的二维单元面积或梯形权重重新计算。