基于物理信息神经网络的传热过程物理场代理模型的构建

陆至彬，瞿景辉，刘桦，何畅，3，张冰剑，3，陈清林，3

（1 中山大学材料科学与工程学院，广东广州510275； 2 中山大学化学工程与技术学院，广东珠海519082； 3 广东省石化过程节能工程技术研究中心，广东广州510275）

引 言

复杂流体在自然界和工业过程中普遍存在，精确地模拟流体基本传热行为在化工、热力和航空航天等诸多学科领域的研究中不可或缺。通常，研究人员需要测量大量的实验数据（例如温度、速度、压力等）来构建模型[1]，以便进行参数分析、性能表征和模型优化等工作。传统的机器学习方法（深度/卷积/递归神经网络）[2-4]可以最大化利用上述数据蕴含的信息[5]，构建高保真的传热模型以准确描述“输入—输出”间关联关系，进而可用于状态监视、产品质量控制和风险预测。但是，对于大部分具有微观、封闭特征的复杂换热器件（如：微尺度换热器和高温辐射换热器），其内部物理场的关键信息数据的测量成本高昂甚至难以测量，研究人员往往需要在数据稀疏的小数据条件下构建模型或者做出决策，这导致绝大多数传统机器学习方法缺乏鲁棒性，无法提供收敛的保证[6-7]。另外，需要注意的是换热器件内部流体系统通常受高度非线性的Navier-Stokes 方程组控制，物理信息间受到基本物理规律的约束，盲目地以数据驱动方法学习物理场信息，可能会得到难以解释的甚至违背物理规律的结果。

近年来，随着深度残差神经网络（ResNet）[8-9]以及物理信息神经网络（PINN）[7]的相继提出，深度学习算法在求偏微分方程和认识物理规律方面取得了重大进展，突破了传统机器学习受限于数据驱动和黑箱建模的局面。特别是，PINN将神经网络的预测结果约束在物理规律之内，使机器学习方法摆脱了对实验或模型数据的根本性依赖，在提高精度同时也大大提高了模型的可理解性。例如：Raissi等[7,10-11]将PINN 用于求解确定性的一维Burger 偏微分方程，以及求解具有中等数量训练数据的二维/三维偏微分方程约束反问题，可反演推导Navier-Stokes 方程中的不确定参数。 在此基础上Karniadakis等[12-14]成功将PINN应用于反应扩散方程和欧拉方程的求解。Wang 等[15]利用PINN 来预测基于雷诺时均方程模型的流体雷诺应力的差异，并取得了良好的预测效果。类似地，PINN还被应用于学习达西流中的本构关系[16]、材料的原子模拟[17]、多物理场亚表面传输问题[18]，以及不确定性量化问题[19-21]。

PINN耦合了换热器件内部的物理规律，因而仅仅需要实际测得的边界条件即可准确预测整个物理场的信息。但仍需注意到，在训练PINN 时使用不同的边界条件设置方法，对训练结果有着显著的影响。PINN 的损失函数包含域损失项和边界损失项，通常损失函数通过“软边界”的方式设置边界损失项，由惩罚系数来控制边界损失项和域损失项的权重，以加速优化问题的收敛。但是，惩罚系数的选择往往依靠经验来调整，不当的惩罚系数容易导致求解结果出现畸形解。针对上述问题，Sun等[22]引入距离函数，使用“硬边界”的方式设置边界条件，将边界损失和域损失统一在同一项中，避免了惩罚系数的使用，并提升了预测精度。

PINN提供了求解偏微分方程的新途径，同时其耦合物理定律的特性呈现出其直接作为代理模型的潜力[22-25]。本工作基于PINN 理论，以二维空间坐标和控制方程的参数构成三维输入，预测了物理场的云图分布，并与商业软件模拟结果进行对比验证。在这一过程中，首先通过简单的具有内热源的二维稳态导热案例分析对比软边界PINN 与硬边界PINN 的性能，接着将硬边界PINN 应用到板式换热器的流动传热案例当中，讨论其在构建复杂物理场代理模型方面的性能。

1 基于物理信息的神经网络

1.1 基础框架

物理过程均可以用多个偏微分方程进行描述，求解得到的原函数，无论其形式是解析解还是数值解，均包含了大量的物理场信息。对原函数T(0)(x,y,θ)，其二维偏微分方程g(T)如下：

其中，T(n)(x,y,θ)表示T(0)(x,y,θ)的n阶导数项，n={0,1,2}，C为常数，x、y为坐标，θ为参数。

如图1 所示，求解T(0)(x,y,θ)可按以下三个步骤进行：（1）首先将输入{x,y,θ}转换成神经元的形式，生成初始的神经网络t0(x,y,θ)；（2）借助自动微分算法[26]求t0(x, y, θ)的梯度，构成偏微分方程g(t0)，并以g(t0)与g(T)之间的偏差模量作为损失函数；（3）最后利用梯度下降算法[27]，调整神经网络的权重w 和b，使g(t0)逼近g(T)，同时也使t(x, y, θ)逼近T(0)(x, y, θ)。多次迭代训练后得到的t(x, y, θ)即是T(0)(x, y, θ)的PINN代理模型。

图1 PINN结构示意图Fig.1 PINN structure diagram

其中，上角标l 代表神经网络的层序号，下角标i，j代表神经元索引，例如：和分别为第l-1 层中第i 个神经元到第l 层中第j 个神经元的权重和偏置；a=σ(z)为激活函数[28-29]，可选择双曲正切函数或sigmoid函数：

其中，z 和a 分别为神经元的输入值和输出值。最后一层第L层为输出层，输出拟合函数t：

由此，基于式(2)给定任意的输入值即可快速计算出t(x,y,θ)。

1.2 软边界PINN

神经网络各层之间的权重、偏置可以通过最小化均方误差损失来学习，定义损失函数为

对于给定的参数，{xb,yb,T(xb,yb)}为用于训练的边界点数据，{xf,yf}为区域内部用于训练的配置点数据，λ 是惩罚系数。损失函数Loss 第一项对应于边界数据；而第二项为充当惩罚作用的正则化项，保证了神经网络在有限的配置点集合上强制执行偏微分方程式(1)所施加的物理定律约束。例如，在不可压缩流体动力学问题中，任何违反质量守恒原理的流动解都会被丢弃，从而将解空间约束在物理框架下[7,30]。

1.3 硬边界PINN

软边界PINN 对边界的处理采用的是有约束优化，其局限性在于损失函数的优化性能取决于每个项的相对大小。因此，如何确定λ 值使得损失函数每一项在优化过程中处于相同的数量级是很有挑战性的工作。Sun 等[22,31]提出一种灵活处理边界条件的方法，即构造一个满足边界条件特解的函数，具体的表达式如下：

其中，G(x,y,θ)为边界函数，是边界数据的光滑扩展函数，D(x,y,θ)为距离函数，是与内部点到边界最小距离有关的光滑函数。

对于Dirichlet 边界条件，边界上各点的边界值为T(xb, yb, θ)，坐标为(xb, yb)。此时，由于边界上的D(xb, yb, θ)=0，因此有h(xb, yb, θ)=G(xb, yb, θ) = T(xb, yb,θ)。而对于不在边界的内部点(xf,yf)，其D(xf,yf,θ)取值不为0，则h(x,y,θ)受到边界函数和距离函数共同控制。通过将Dirichlet边界条件耦合在所构造的扩展函数h(x, y, θ)当中，损失函数被进一步简化为只包括偏导数的形式：

其中，{xt,yt}为用于训练的所有点数据（包括边界）。在上述情况下，距离函数与DNN 输出的乘积D(x,y,θ)aL(x,y,θ;w,b)在数学上相当于此偏微分方程的通解，而边界函数则包含了偏微分方程的特解。这种“通解+特解”的数学结构，赋予求解结果自动满足特解的先行条件，加速PINN 向合理梯度方向收敛。

值得注意的是，上述推导是基于Dirichlet 边界条件，对于梯度形式的Neumann 边界（h(1)(x, y, θ)），则将边界梯度直接引入损失函数当中，以硬边界PINN损失函数为例：

式中，{xB, yB}为用于训练的Neumann 边界点数据，C0为Neumann边界处的梯度值。

2 案例一：具有内热源的二维稳态导热方程

2.1 案例模型

本节首先以一个简单的案例探究软、硬两种边界条件设置方式下，PINN预测结果的差异。以具有内热源的二维稳态导热方程为例，其表示式如下

其中，q为热源，κ为热导率。案例选取边长为1的正方形区域作为研究对象，以方形中心为原点建立笛卡儿坐标系。

假设热源保持稳定（q=1），材料热导率在{κ ∈R:0.5 ≤κ ≤1}内变化，相应的参数θ 变化范围为{θ = q/κ ∈R:1 ≤θ ≤2}。以空间坐标和控制方程的参数{x, y, θ}作三维输入，建立PINN 代理模型。训练集中，样本坐标值（x 和y）使用均匀采样生成51×51 的网格，参数θ 在区间内参数化，样本量为10。对于软边界PINN，边界上使用Monte Carlo采样方法生成1000 个样本。此外，设置PINN 含有3 个隐藏层，每个隐藏层含有50 个神经元，激活函数为sigmoid函数，见式(4)。

软边界PINN 的损失函数设置为式(6)的形式，且方形四周的边界温度恒为1。而硬边界PINN 的损失函数设置为式(8)的形式，其中边界函数G(x,y,θ)和距离函数D(x,y,θ)分别如下：

2.2 温度场预测结果

本工作所涉及的采样方案和神经网络训练均基于Python3.7.7 和Tensorflow1.14 环境，在Spyder 中编译完成。有限元模拟均使用COMSOL 5.5。计算机配置为Intel i7-8550 CPU。训练后使用θ=1.5 时的温度场与COMSOL 模拟结果对比，以验证PINN的预测效果。如图2所示，两种PINN 方法构建的代理模型预测结果均与模拟结果非常接近，均准确地捕捉到了温度的变化趋势。结合相对误差场图（软边界PINN 见图3，硬边界PINN 见图4），软边界PINN 和硬边界PINN 最大相对误差分别为0.60%和0.015%，再次表明两种PINN 方法训练得到的代理模型都具有很高的精度，但硬边界PINN 的相对误差要更低。

图2 PINN代理模型和COMSOL模拟对温度场的预测结果对比Fig.2 Comparsion of temperature fields predicted by soft PINN,hard PINN,and COMSOL simulation

图3 软边界PINN温度等值线的预测结果对比Fig.3 Comparsion of temperature contours predicted by soft PINN

图4 硬边界PINN温度等值线的预测结果对比Fig.4 Comparsion of temperature contours predicted by hard PINN

在分别对比软边界和硬边界与COMSOL 模拟结果的等值线图可见，中心区域附近，使用这两种边界设置方法的预测结果均与模拟结果基本重合。但仍需注意到，软边界PINN 在四个角落区域的预测效果较差，而硬边界PINN 在四个角落区域的预测效果依然良好。造成这种现象的原因可能是软边界PINN 用于训练的边界点数据量较少。为此，将边界点的样本量增加到4000后重新训练，预测精度提升13.3%，但是角落区域预测误差较大的问题依然存在。此外，硬边界PINN 的计算时间更小。硬边界PINN 所用的训练时间为7 min，远低于软边界PINN 所用的训练时间39 min（1000 个边界样本点）和55 min（4000 个边界样本点）。因此可以认为，硬边界PINN要优于软边界PINN。

3 案例二：平板间的二维稳态对流传热方程

3.1 案例模型

在换热器件内部，动量传递、热量传递和质量传递往往联系密切，因此研究多物理场对指导实际换热过程具有重要的意义。例如在板式换热器中，平板间的温度传递不仅由热传导引起，还受到流体流动影响，本节以此为例利用PINN 来预测二维稳态对流传热方程。所使用模型几何如图5 所示，换热工质在两块平板间流动，并存在热传导和热对流。假设壁面无滑移，且两平板间的不可压缩稳态流动已经充分发展，则垂直于平板方向的速度可由Navier-Stokes方程推导得出抛物线形的解析解：

其中，u 表示速度场的x 分量，umax为x 分量的最大速度，R为平板壁面到中轴线的垂直距离。

对于有内热源的稳态层流传热，忽略黏性耗散后，其能量守恒式可简化为

其中，ρ 和Cp分别表示流体的密度和比定压热容，显然该方程是在方程式(10)的基础上引入了第一项对流项。假设平板隔热性能良好，与外界环境绝热，则在平板上下壁面需设置Neumann 边界条件：

本案例采用硬边界PINN 的方法建立代理模型，假设进口温度的Dirichlet边界Tin已知，则可构造以下满足条件的边界函数和距离函数：

其中，L 代表平板长度，x∈[0, L]。基于这一定义，可构造方程式(7)所示的扩展函数：

令空气为板间流体，物性参数为：ρ=1.239 kg/m3，Cp=1,005 J/(kg·K)，κ=0.024 W/(m2·K)。进口温度Tin=30℃，平板间的通道尺寸为L=1 m，R=0.05 m。以{x,y,umax}作为输入变量的三个维度。其中，x，y为两平板间均匀生成的50×100 网格的坐标。最大速度umax在区间{0.01～0.03 m/s}内参数化，取10 个样本。设置PINN 隐藏层4 层，每层含有50 个神经元，使用sigmoid函数[式(4)]作为激活函数。

3.2 温度场预测结果

训练代理模型所用的时间为54 min，训练后损失函数值为2.6×10-4。使用umax=0.02 m/s时的温度场验证预测效果，如图6 所示。由图可知硬边界PINN构建的代理模型所预测的流体流动过程中的温度变化趋势与实际模拟结果是一致的。值得注意的是，在相同的横坐标下，流体在中间的温度要低于两侧，而中间的流速高于两侧，由此说明流动加速了流体的热量交换，流速越大热交换越明显。此外，温度场表现出与速度场类似的对称性，再次说明稳态热量传递与动量传递密切相关。图6的温度场绝对误差云图表明，平板区域内的温度预测误差基本接近于0，仅在四个角落处误差略微变大，但依然在10-3数量级内，最大绝对误差为0.0038℃，再次凸显了硬边界PINN 构建的代理模型精度高的优势。

图5 平板的几何结构示意图Fig.5 Geometric diagram of the flat plate

图6 COMSOL模拟温度场与硬边界PINN预测温度场对比Fig.6 Comparsion of temperature fields predicted by hard PINN and COMSOL simulation

图7 硬边界PINN代理模型（红）与COMSOL模拟（黑）预测温度等值线结果对比Fig.7 Comparsion of temperature contours predicted by hard PINN(red line)and COMSOL simulation(black line)

硬边界PINN 预测的温度等值线与COMSOL 模拟的温度等值线在图7 中给出。显然，等值线在平板间中轴线处的延伸趋势比壁面处的更大，这一变化直观地呈现了流体流动对温度分布的影响。需要指出的是，PINN预测与COMSOL模拟结果的差距沿流体流动方向扩大。在进口处两者的温度等值线基本重合，而在远离进口处两者的差异逐渐增加，并且在中轴线附近的差距最大。这是由于速度场的非线性造成温度场非线性的增加，这会使代理模型中产生更大的偏差，这也是未来PINN 需要面对的主要挑战。

4 结 论

本工作将PINN 模型引入到换热物理场的代理模型的构建中，预测给定参数下的温度场分布，并与商业软件的模拟结果进行验证。通过简单的具有内热源的二维稳态导热案例分析，软边界PINN和硬边界PINN 最大相对误差分别为0.60% 和0.015%，结果表明对边界条件进行耦合的硬边界PINN 在建立代理模型时更胜一筹。而将硬边界PINN应用到板式换热器的流动传热案例当中，所建立的代理模型的预测结果与模拟结果吻合良好，温度场最大绝对误差为0.0038℃，再次证明了硬边界PINN构建高精度代理模型的能力，并为换热器的开发提供了良好的工具。