☰
软物理信息神经网络求解平板间对流传热方程的实践与踩坑
2026/10/18 5:51:14 网站建设 项目流程

说实话,第一次看到“平板间二维稳态对流传热方程的软物理信息神经网络实现”这个题目时,我的第一反应是:这不是一个已经有经典级数解的老问题吗,为什么还要上神经网络?但真的把软物理信息神经网络(Soft PINN)用在传热问题上折腾过几轮之后,我才发现这类看似“传统”的传热问题,恰恰是检验软约束策略的绝佳试验场。入口处的热边界层很薄,出口是开边界,上下壁面温度跳变剧烈,任何一个环节处理不干净,损失函数都会给你一顿毒打。这篇博文记录的就是我如何用软物理信息神经网络求解平板间二维稳态对流传热方程的一次完整实践,包括方程推导、无量纲化、Python代码实现、让人劝退的踩坑过程,以及最终如何让网络学出一个合理温度场的全部经验。

1. 为什么经典问题还需要“软”求解器?

1.1 传统方法哪里让人不舒服

如果只讨论“计算平板间稳态温度场”这件事,传统有限差分或者有限体积法在规则几何下其实非常强悍,结构化网格一发,边界层加密一做,结果又快又准。那PINN的价值到底在哪里?我个人的体会是:它的价值不在正向计算的精度,而在于把正向问题、反向问题和多源数据融合放进同一个框架里。

举个实际场景:如果我不满足于给定边界条件去求温度场,而是反过来——通过部分壁面温度测量数据,反推入口温度分布或者壁面热流,传统求解器通常需要伴随方法、梯度求解,或者反复迭代调用正向求解器,工程成本不低。而PINN只需要在损失函数里面加一项数据残差,训练完成之后,网络不仅拟合了测量数据,还能自动满足控制方程。这种“同一套代码,换物理问题就换损失项”的方式,对工程人员来说有实实在在的吸引力。

1.2 硬约束与软约束的本质差异

PINN处理边界条件有两种路线。硬约束的思路是:直接把网络输出构造成一个解析表达式,让它在边界上自动满足条件。比如算一个温度场时,可以令预测值等于边界函数加上一个到边界的距离函数与网络输出的乘积。这样网络输出再怎么乱跑,边界条件都严格成立,网络只需要专心最小化PDE残差就行。

听起来很美,但实际操作时会遇到两个麻烦:一是构造试探解高度依赖几何外形,平板这种简单几何还好说,稍微复杂一点的曲线边界,距离函数怎么写就够折腾半天;二是我们这里有一个出口的Neumann边界条件,也就是∂θ/∂X = 0,这种法向导数边界条件很难用简单的乘积形式光滑地嵌进网络结构里。一旦以后边界条件从定温改成热流、从入口定温改成抛物线分布,整个试探解又得推倒重来。

软约束则简单粗暴得多:网络直接输出裸的温度场,边界条件不强制满足,而是作为正则化项放进损失函数里,通过权重系数压着网络去逼近。它的本质是最小化加权残差,也就是把偏微分方程约束和边界约束揉成一个多目标优化问题。这么做的好处是代码极其通用,改边界条件基本只是改几行采样和损失项的事;坏处是边界条件只能渐近满足,权重、采样、学习率稍微配合不好,就会出现边界对、内部错,或者内部尚可、边界离谱的尴尬局面。

1.3 为什么平板间传热适合当软约束的基准问题

选这个题目做实践,不是因为平板问题本身需要PINN,而是因为它作为软约束策略的基准测试非常严苛。

第一,平板几何足够简单,参考解可以独立求解,我们能够严格量化软PINN预测的误差到底有多大,而不是像某些复杂工程问题那样,算出来也不知道准不准。第二,这个物理问题里有两个典型难点:入口角点存在温度不连续,出口是开放边界。这两个地方恰好都是软约束容易暴露问题的地方,把这些问题在一个可控的benchmark上调通,换到更复杂的几何和边界条件时,心里才有底。第三,传热方程的PDE残差涉及二阶导数,对自动微分的稳定性和网络激活函数的光滑性要求都很高,能检验一套PINN代码的真实功底。

2. 平板通道的对流传热方程:从物理到无量纲形式

2.1 物理假设与控制方程推导

我考虑的场景是两片无限宽平行平板构成的水平通道,下板温度高,上板温度低,流体沿x方向流动,在流动过程中被上板冷却。做如下简化假设:

  • 常物性、不可压缩、层流,流动已经充分发展,速度剖面固定为抛物线分布;
  • 忽略黏性耗散、体积力、内热源;
  • 流动与传热均为二维稳态,z方向无限长,不考虑展向变化。

在这个前提下,二维稳态能量方程可以写为:

ρcp (u ∂T/∂x + v ∂T/∂y) = k (∂²T/∂x² + ∂²T/∂y²)

因为充分发展的平板通道流中法向速度v = 0,再考虑Peclet数较大时轴向分子导热相对轴向对流可以忽略,方程进一步简化为:

u(y) ∂T/∂x = α ∂²T/∂y²

其中α = k/(ρcp)是热扩散系数。这里速度剖面取无量纲平均速度归一化后的形式:

u(y) = 6 Um y(H − y) / H²

其中H是板间距,Um是截面平均速度。选择这个剖面的原因是它对y从0到H积分正好是Um,方便后续无量纲化。

2.2 无量纲化过程与关键参数

为了减少参数对网络训练的干扰,我把方程改写成无量纲形式。定义:

θ = (T − Tc) / (Th − Tc)

X = x / L

Y = y / H

U = u / Um

其中Th是下板温度,Tc是上板温度,L是通道加热段长度。代入简化方程:

(Th − Tc)/L · U ∂θ/∂X = α (Th − Tc)/H² · ∂²θ/∂Y²

整理得到无量纲控制方程:

U ∂θ/∂X = C ∂²θ/∂Y²

其中无量纲系数:

C = αL / (Um H²)

这个C非常关键,它本质上是轴向对流与法向导热的竞争关系。C越大代表法向扩散越强,温度场更容易被抹平;C越小代表对流占主导,入口段热边界层就越薄,网络需要解析的梯度就越陡峭。

用一组工程参数来感受一下量级:空气的α大约2×10⁻⁵ m²/s,H取0.02 m,Um取0.1 m/s,L取1 m,那么C = 2×10⁻⁵ × 1 / (0.1 × 0.0004) = 0.5。如果流速提高一个量级,C就变成0.05,热边界层明显变薄。我在后面的实验中把C控制在0.05到0.1之间,这样既能看到清晰的入口段温度发展,又不会让边界层薄到完全无法训练。

2.3 边界条件与参考解的构造

无量纲后的计算域是X ∈ [0,1],Y ∈ [0,1],边界条件为:

  • 入口 X = 0:θ = 1
  • 下板 Y = 0:θ = 1
  • 上板 Y = 1:θ = 0
  • 出口 X = 1:∂θ/∂X = 0

这里有一个物理上的不连续点:入口与上板的交界处,也就是坐标(X,Y) = (0,1),一边要求θ = 1,一边要求θ = 0。在真实物理中,这个温度跳变会被法向导热瞬间抹平,形成一段极陡的温度过渡区。软PINN不会也无法精确复现这种间断,它只能给一个光滑近似,这个光滑近似的宽度由网络分辨率和采样密度决定。

为了评估网络效果,我没有直接用传统CFD软件,而是用一个高分辨率的迭代有限差分程序,在500×200的网格上求了数值参考解。计算时把入口上角点略微做了一点过渡处理,避免数值振荡。这个参考解被保存成numpy数组,训练完成后用来计算相对误差和逐剖面比对。在软PINN的验证流程里,独立参考解是必须的,否则光看损失函数下降根本判断不了学出来的温度场是物理正确还是数值幻觉。

3. 把方程塞进网络的损失函数:结构、微分与权重

3.1 网络结构怎么选才适合二阶PDE

网络结构方面,我采用了一个标准的全连接MLP,输入是二维坐标(X,Y),输出是一维无量纲温度θ。结构上我选了四层隐藏层,每层64个神经元,这属于一个折中方案:太浅拟合不了入口段热边界层的陡峭变化,太深则训练难度和计算成本都上去了。

比层数更重要的是激活函数。这里强烈不建议用ReLU,原因是PDE残差里包含∂²θ/∂Y²,ReLU的二阶导几乎处处为0,网络根本无法通过损失反向传播给二阶导数项提供有效梯度。Tanh因为二阶导连续且非零,成为处理这类PDE残差的默认选择。

如果你发现Tanh收敛太慢,可以试试正弦激活函数,它在某些边界层问题里有奇效,但学习率要调小。我的经验是:先老老实实用Tanh,把损失、权重、采样这几件事理顺了,再去做激活函数的横向对比才有意义,否则你根本分不清是结构问题还是调参问题。

3.2 损失函数:PDE残差、边界残差与自动微分

软PINN的损失函数是一个加权求和。在没有实测数据的情况下,我用的形式是:

L = λr · Lr + λbc · (L_inlet + L_low + L_up) + λout · Lout

其中各项的含义如下:

  • Lr是内部采样点上的PDE残差均方值,也就是(U ∂θ/∂X − C ∂²θ/∂Y²)²的均值。
  • L_inlet、L_low、L_up分别是入口、下板、上板的Dirichlet边界残差均方值。
  • Lout是出口Neumann边界残差,也就是(∂θ/∂X)²的均值。

计算这些残差依赖自动微分。PyTorch里要用两次torch.autograd.grad来求二阶导数,第一次对网络输出求一阶导,第二次对一阶导结果再次求导。这里有一个代码上的坑:二阶导必须在第一次求导时设置create_graph=True,否则计算图会被释放,第二次求导直接报错。我第一次写的时候漏了这个参数,报错信息一度看得我摸不着头脑。

3.3 权重系数:边界优先还是残差优先

权重系数是软约束的老大难问题。最开始我采用了λbc = 1、λr = 1的“公平”配比,结果训练出来的温度场在边界附近一律变形。原因很简单:边界残差的误差会通过光滑函数传播到整个计算域,如果边界一开始就没压住,PDE残差再小也没有意义。

后来我调整了策略:初期让λbc = 20左右,λr = 1,先把边界条件学个大概;等边界损失降到和残差损失同量级之后,再把λr逐步提高到5到10。这背后的逻辑相当于让网络先记住“在哪里该是多少度”,然后再去抠“内部的方程到底对不对”。边界都学不对的时候,强行加大PDE权重只会让网络在全局范围内制造一个温和而持续的偏差来换取残差下降。

这里还要提醒一个尺度陷阱:当C特别小,比如0.001时,PDE残差里导热项的量级远小于对流项,网络几乎感知不到二阶导的影响。对策是把原方程重整化,残差写成U ∂θ/∂X / C − ∂²θ/∂Y²,把导热项前面的系数人为归一。我在实验中把C控制在0.05左右,残差各项量级不至于差太多,但方程重整化这个思路仍然值得记在脑子里。

4. Python实现:一个可以直接运行的稳态对流传热PINN

4.1 网络与采样代码框架

下面这段代码是我最终跑通的精简版本。它不是一个教学用的伪代码,而是一个可以直接运行、用来验证核心思路的最小实现。首先定义网络和采样函数:

import torch import torch.nn as nn device = torch.device("cuda" if torch.cuda.is_available() else "cpu") class MLP(nn.Module): def __init__(self, layers=[2, 64, 64, 64, 1]): super().__init__() seq = [] for i in range(len(layers) - 2): seq.append(nn.Linear(layers[i], layers[i + 1])) seq.append(nn.Tanh()) seq.append(nn.Linear(layers[-2], layers[-1])) self.net = nn.Sequential(*seq) def forward(self, xy): return self.net(xy) def sample_points(n_internal, n_boundary): # 内部点:x均匀,y向壁面方向加密 x_in = torch.rand(n_internal, 1, device=device) z = torch.rand(n_internal, 1, device=device) y_in = 0.5 + 0.5 * torch.tanh(2.0 * (z - 0.5)) # 入口 x = 0 x_ent = torch.zeros(n_boundary, 1, device=device) y_ent = torch.rand(n_boundary, 1, device=device) # 下板 y = 0 x_low = torch.rand(n_boundary, 1, device=device) y_low = torch.zeros_like(x_low) # 上板 y = 1 x_up = torch.rand(n_boundary, 1, device=device) y_up = torch.ones_like(x_up) # 出口 x = 1 x_out = torch.ones(n_boundary, 1, device=device) y_out = torch.rand(n_boundary, 1, device=device) return (x_in, y_in), (x_ent, y_ent), (x_low, y_low), (x_up, y_up), (x_out, y_out)

内部点采样时,我用了tanh变换把更多采样点压向上下壁面。这个操作很多人会忽略,但入口段的热边界层恰恰聚集在壁面附近,均匀采样会导致壁面附近的PDE残差约束不足,网络就会在那里自由振荡。

4.2 残差计算与损失函数

接下来是PDE残差和边界损失。注意看autograd.grad的用法,尤其是create_graph的设置:

def pde_loss(model, xy, C): x, y = xy x = x.clone().requires_grad_(True) y = y.clone().requires_grad_(True) xy_cat = torch.cat([x, y], dim=1) theta = model(xy_cat) grad1 = torch.autograd.grad( theta, xy_cat, grad_outputs=torch.ones_like(theta), create_graph=True )[0] theta_x = grad1[:, 0:1] theta_y = grad1[:, 1:2] grad2x = torch.autograd.grad( theta_x, xy_cat, grad_outputs=torch.ones_like(theta_x), create_graph=True )[0][:, 0:1] grad2y = torch.autograd.grad( theta_y, xy_cat, grad_outputs=torch.ones_like(theta_y), create_graph=True )[0][:, 1:2] U = 6.0 * y * (1.0 - y) r = U * theta_x - C * grad2y return (r ** 2).mean()

这里有三个容易出错的细节。

第一,x和y必须先设置为叶子张量再cat。很多初学者直接对torch.cat的结果调用requires_grad_,结果RuntimeError说非叶子张量不能请求梯度。正确做法是让参与cat的两个叶子张量各自requires_grad_,然后再拼接。

第二,需要求二阶导时,第一次grad调用必须带上create_graph=True,否则计算图销毁,第二次求导就没得玩了。

第三,grad_outputs不能少。很多人写成torch.autograd.grad(theta, xy_cat, create_graph=True),这里默认输出的梯度是标量对输入的梯度,但theta是向量,直接这样会报错,必须提供和theta形状一致的ones作为grad_outputs。

边界损失函数如下:

def boundary_losses(model, points, C): _, (x_ent, y_ent), (x_low, y_low), (x_up, y_up), (x_out, y_out) = points # 入口 θ = 1 theta_ent = model(torch.cat([x_ent, y_ent], dim=1)) loss_ent = ((theta_ent - 1.0) ** 2).mean() # 下板 θ = 1 theta_low = model(torch.cat([x_low, y_low], dim=1)) loss_low = ((theta_low - 1.0) ** 2).mean() # 上板 θ = 0 theta_up = model(torch.cat([x_up, y_up], dim=1)) loss_up = (theta_up ** 2).mean() # 出口 ∂θ/∂X = 0 xo = x_out.clone().requires_grad_(True) yo = y_out.clone().requires_grad_(True) xy_out = torch.cat([xo, yo], dim=1) theta_out = model(xy_out) dtheta_out = torch.autograd.grad( theta_out, xy_out, grad_outputs=torch.ones_like(theta_out), create_graph=True )[0][:, 0:1] loss_out = (dtheta_out ** 2).mean() return loss_ent + loss_low + loss_up, loss_out

4.3 训练循环与参数配置

训练阶段我使用的是Adam优化器,初始学习率1e-3,每6000轮衰减一次,乘以0.1。每轮训练都重新随机采样,相当于给网络制造了源源不断的新数据,对防止过拟合参考解或者陷入局部解都有帮助。训练轮数设置在15000左右,内部点4000个,每条边界线300个点。

def train(model, C=0.05, n_internal=4000, n_boundary=300, lr=1e-3, epochs=15000, lambda_bc=20.0, lambda_r=1.0): opt = torch.optim.Adam(model.parameters(), lr=lr) sched = torch.optim.lr_scheduler.StepLR(opt, step_size=6000, gamma=0.1) for it in range(epochs): opt.zero_grad() pts = sample_points(n_internal, n_boundary) loss_r = pde_loss(model, pts[0], C) loss_bc, loss_out = boundary_losses(model, pts, C) loss = lambda_r * loss_r + lambda_bc * loss_bc + 1.0 * loss_out loss.backward() opt.step() sched.step() if it % 2000 == 0: print(f"iter {it:6d} | Lr {loss_r.item():.3e} | " f"Lbc {loss_bc.item():.3e} | Lout {loss_out.item():.3e}") model = MLP().to(device) train(model, C=0.08)

训练完成后,在均匀网格上做一次前向推理,就能画出温度云图。评估时需要注意:训练时采样点是随机的,评估时则要固定在规则网格上,我一般用200×100的网格,然后和参考解做逐点相对误差:

xx, yy = torch.meshgrid( torch.linspace(0, 1, 200), torch.linspace(0, 1, 100), indexing="xy" ) xy_eval = torch.cat([xx.reshape(-1, 1), yy.reshape(-1, 1)], dim=1).to(device) with torch.no_grad(): theta_pred = model(xy_eval).reshape(200, 100).cpu().numpy()

如果说整个训练过程只能留下一句话,那就是:同时监控各项损失,而不是只看总损失。总损失下降不代表物理正确,各项损失分开打印才能看出是哪一部分在拖后腿。

5. 踩坑记录:权重失衡、入口尖角与出口假约束

5.1 残差降了,温度场却振荡

我第一次训练时,总损失一路降到非常漂亮的数量级,但画出来的温度云图在入口下游一段距离出现了波浪形振荡。这个结果差点让我怀疑整个PINN方法,后来才发现问题出在采样密度上。

均匀采样时,近壁区域的点太少,而入口段热边界层的温度梯度非常大。网络在梯度大的地方如果缺少PDE残差约束,就会表现得不稳定,可以理解为它找到了一个全局折中,让边界和大部分区域都大致正确,却允许壁面附近的小范围区域任意振荡。解决方案就是我在采样函数里写的tanh加密变换。这个细节看似不起眼,却是整个实验从“能训练”到“能看”的分水岭。

5.2 入口角点的不连续处理

入口和上板的交界处,也就是(X,Y)=(0,1),存在温度从1跳到0的不连续。如果你不做任何处理,直接把采样点打到这个角点上,网络会被迫用一个极陡的过渡去拟合矛盾条件,进而污染整个入口区域。

我的做法有三重保险:一是在入口边界采样时,让角点附近稍微多留一点空隙,不把采样点精确压在角点上;二是在角点附近的内部区域加密采样,给网络充分的自由度去构造光滑过渡;三是只要工况允许,就把入口温度条件改成一个带过渡的解析形式,比如θ(0,Y)=1−Y^m,用m控制过渡宽度,这样物理上更平滑,网络也更好学。

5.3 出口Neumann约束的物理矛盾

出口处我要求∂θ/∂X=0,但从物理上讲,这个条件只在出口截面足够远离入口、温度场充分发展时才严格成立。如果通道不够长,出口截面本身还存在明显的轴向温度梯度,这个边界条件就和真实物理产生了冲突,网络会被逼着同时满足两个矛盾的目标。

解决这个问题有两个思路。第一个是从物理建模入手,把计算域加长,比如只观察前60%区域的温度场,出口设在更下游,让出口截面上的温度梯度确实趋于0。第二个是从权重入手,把出口损失权重调小,允许网络在出口附近有一点偏差,优先保证入口段和壁面附近的主流区域正确。第二个思路听起来有点“和稀泥”,但在软约束框架里这就是真实可用的策略。

5.4 学习率与收敛平台期

深度网络在PDE约束下的训练有一个典型现象:前期损失下降很快,中期进入平台期,看似卡死了,但如果你把学习率降一个数量级,又能继续下降。所以我的训练循环里用StepLR做分段衰减,而不是一个学习率用到底。

还要注意一个细节:每轮重新采样会让损失曲线不是单调下降,而是在一个范围内波动。不要被这种波动吓到,只要波动的中心在缓慢下移,训练方向就是对的。如果某一次训练你看到损失曲线完全平滑,反而要怀疑是不是固定了一批采样点导致网络走进了过拟合的死胡同。

6. 从软约束出发还能做什么

6.1 自适应残差加权与重要性采样

固定权重只是软约束的入门玩法。更高效的做法是让采样点和权重跟着残差分布自适应变化:每训练几百轮统计一次PDE残差的绝对值,然后在残差大的区域增加采样点,这相当于让网络把注意力集中在最难学的地方。对入口段这种局部高梯度问题,自适应采样往往能让收敛速度翻一倍。

我后来在代码里加了一版简单的残差加权:每个点的损失不再是均等权重,而是乘上(|r| + 1)的幂次,相当于给残差大的点更高的惩罚权重。这种变体训练出来的温度场,在相同迭代次数下明显更光滑。

6.2 数据与物理的混合:反向问题

如果实验中有热电偶测得若干个点上的温度,软PINN只需要在损失函数里新增一项数据损失:

Ldata = (1/Ndata) Σ |θ_pred(x_i, y_i) − θ_meas_i|²

数据项会和物理方程互相约束,实测数据有噪声时,权重不要开太大,给它和边界条件差不多的量级就行。这种混合方法在CFD里对应的是数据同化,在PINN里只是加一行代码的事,这是软约束框架最迷人的地方。

6.3 从稳态到瞬态,从单场到耦合

在稳态问题上跑通之后,往瞬态扩展其实非常自然:网络输入从二维坐标变成三维坐标,在损失函数里加一项∂θ/∂t的残差就行。真正的难点在于时间方向的因果性,直接在所有时间上同时训练,网络会把时间远处的信息泄漏到早期时刻。文献里常见的做法是按时间窗切片依次训练,先学前0.1秒,再学0.1到0.2秒,逐步外推。

至于从单物理场到流固耦合或者自然对流的完整NS方程,思路是类似的,只是残差函数更复杂,权重更难调,训练时间更长。但底层那一套“把物理方程变成多目标损失”的框架没有变。

如果让我给一个总结性的个人建议,那就是:做这类PINN实践,真正值钱的不是神经网络本身,而是对物理模型尺度的把握和对边界条件一致性的理解。网络只是拟合工具,物理没建对、边界自相矛盾,再漂亮的损失曲线也只是数字游戏。从平板间稳态传热这种小问题入手,把每一个坑都踩明白,换成任何更复杂的PDE和应用场景,你都会多一分底气。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询