【科学计算】使用飞桨高阶自动微分功能探索AI+结构领域科研

在工程和科学领域，AI与科学方法的结合正在解决更多经典科学问题，并在固体、流体、传热、材料等越来越多的领域得到了验证。为了更好支撑科研人员开展AI与基础学科的交叉融合研究，百度飞桨不断完善框架能力，提供支持AI+科学计算的复数算子、高阶自动微分机制以及高性能编译等等。

本期我们聚焦结构领域中的“高阶”物理机理，向大家说明飞桨如何利用“高阶自动微分算子”处理结构领域中的经典问题。在案例中，我们将介绍如何使用飞桨框架对4阶以上偏微分方程进行无监督求解，同时结合DeepXDE或赛桨PaddleScience等科学计算工具组件获得更好的开发体验。

01 结构领域背景与痛点

在结构领域中，飞机机翼桁架、汽车底盘、轮船甲板等机械结构件的受力变形、破坏以及疲劳损伤等都是最典型的“力学”工程难题，解决好这些经典工程挑战无论对基础科研或工程应用都具有十分重要的意义。

传统的结构分析方法包括有限元法、有限差分法、边界元法等，通常需要大量的计算资源和人力投入，这会限制模型的规模和精度。近年来，基于物理信息神经网络（PINN）的AI方法逐渐被应用于结构领域的物理方程求解，这类方法普遍被认为兼具提升计算速度与降低人力投入的优势。

不同于经典NLP、CV等领域问题，其中多数任务都是基于梯度下降算法进行优化，只需要一阶导数（即梯度）来更新模型参数。物理问题通常需要使用高阶微分方程（ODEs/PDEs/fPDEs/IDEs）及相应的初边值条件来描述物理机理，这要求融合物理机理的AI方法可处理高阶微分。目前深度学习框架中高阶微分算子往往通过手写、组合等方式实现，尤其是4阶以上微分算子的定义难度及实现工作量都非常大。

02 结构领域问题求解原理

结构领域问题描述

结构领域问题通常可用平衡微分方程、几何方程和物理方程来描述，其通用形式一般很难求解。实际应用中多使用简化方程，如广泛存在的平面问题，从而降低求解难度。针对二维线弹性问题，我们可以找到某些标量来代表应力函数，从而可获得同时满足平衡微分方程和相容方程的标量微分方程。在忽略体积力的情况下，二维线弹性问题可由以下偏微分方程描述：

(1)

式中，φ是极坐标系(r,θ)下的Airy应力函数。在直角坐标系中，该偏微分方程可描述如下： 
                           

(2)

式中，(x,y)=r(sinθ,cosθ)。以上方程也被称为双调和方程，它通过一个统一的表达式来描述平面线弹性问题。

理论上，基于方程(1)或(2)以及充分的边界条件，即可实现问题求解。然而，由于方程包含四阶双调和算子∇4，传统数值方法需要使用更高阶的单元和算法才能进行求解，而其计算时间和成本是极高的，大型复杂结构问题中则普遍面临着“维度爆炸”、精度低、求解难等问题。对比传统求解方法，PINN方法具备自适应采样、非线性激活函数、并行计算等特点，可以将PDEs嵌入到神经网络中，可有效地解决高维空间和复杂几何形状的问题。

一旦求解获得Airy应力函数，则可计算应力分量如下：

（3）

采用胡克定律，则可计算应变分量如下：

(4)

式中，对于平面应变和平面应力问题，k分别取3-v和(3-v)(1+v)；μ可取为剪切模量G，v为泊松比。通过以上高阶偏微分方程，我们可以对结构领域的典型问题进行定义。不同于传统的数值差分的计算方法，基于给定的偏微分方程，利用PINN方法可在无监督的情况下直接对结构领域问题进行求解，具体如下。

PINN方法求解结构问题

围绕结构领域中高阶微分方程所描述的工程和科学问题，飞桨框架可以支持PaddleScience和DeepXDE科学计算工具组件构建相应的PINN求解模型，并针对前面描述的二维线弹性结构问题，构建了求解应力函数的PINN模型，如图1所示。

图1 结构领域PINN求解模型

该求解模型引入二维线弹性结构问题的控制方程和边值条件，这些物理信息约束不仅降低了神经网络对数据的依赖，也降低了对网络复杂度的要求，使得采用一个简单的全连接网络即可实现二维线弹性结构问题的求解。

同时为了解决该领域中高维偏微分方程的定义与求解，飞桨一方面在框架中增加高阶微分算子，另一方面构建基础算子体系，可支持不限阶数的自动微分。目前，飞桨已经支持绝大多数算子的无限阶自动微分，以及部分算子的有限阶自动微分，主要算子简述如下：

无限阶自动微分
包括标量与Tensor之间的加、减、乘、除、幂运算，elementwise系列运算、矩阵乘matmul、激活函数tanh，以及assign、concat、cumsum、expand_v2、reverse、squeeze、unsqueeze、scale、tile、transpose、sign、sum、mean、flip、cast、slice等；

有限阶自动微分
sin、cos、sigmoid等算子支持三阶计算。基于飞桨最新提供的无限高阶自动微分（AD）算子，我们也可在飞桨全量支撑的DeepXDE科学计算工具中，通过导入paddle.fluid中的core模块，利用deepxde.gradients.jacobian和deepxde.gradients.hessian可分别计算一阶和二阶微分，以及通过它们的任意组合来实现高阶微分。对此，我们基于飞桨实现了PINN方法求解结构领域中1维梁变形、2维平面受载两类典型问题求解，并得到与理论解一致的结果。

03 典型案例介绍

下面根据前面给出的二维线弹性问题PINN求解模型，介绍基于飞桨联合科学计算工具DeepXDE构建的结构领域典型案例，包括欧拉梁挠度计算、矩形平板集中受载基准案例、圆形平板分布受载基准案例（具体可见DeepXDE结构领域科学计算案例）。

一维欧拉梁问题

欧拉梁（Euler beam）是被用来描述长条形物体弯曲运动的物理模型。在桥梁、飞机机翼、吊桥等众多情况下，都可以用欧拉梁模型来描述这些物体的弯曲、扭转和振动等运动。欧拉梁模型假设物体是细长且刚性的，可以在一个平面上自由弯曲，其几何通常可被简化为一维。DeepXDE的官方目录examples / pinn_forward下包含该类案例的一个具体实例Euler_beam.py，其对应的控制方程如下：

(5)

式中，w为挠度，也即选取的应力函数φ。其边界条件可描述如下：左端点（x=0）固定，因而其扰度和转角为零，即：

（6）

右端点（x=1）承载，因而其弯矩和剪切力为零，即：

（7）

根据以上给出的挠度四阶方程和边界条件，基于DeepXDE可以简单地构建相应的PINN求解模型，主要步骤的代码如图 2所示。

图 2 一维欧拉梁问题的PINN求解模型代码

采用以上代码可以求解欧拉梁的挠度，结果如图 3所示。可以看出，梁的左端因为固定而挠度为0，从左到右逐渐增大。

图 3 欧拉梁问题的PINN求解结果

矩形平板分布受载基准案例

图 4 分布载荷作用下的矩形平板基准案例Vahab M, Haghighat E, Khaleghi M, et al. A physics-informed neural network approach to solution and identification of biharmonic equations of elasticity, 2022.

考虑更复杂的二维矩形平板受分布载荷的基准案例，如图 4所示。平板的长、宽和厚分别为a=2m、b=3m和t=0.01m，厚度为四周被简单支撑，表面则被施加一个正弦分布载荷 ππ ，相应的控制方程可描述如下：

(8)

式中，w为平板挠度，D为抗弯刚度，可计算如下：

(9)

式中，E为弹性杨氏模量，v为泊松比。根据平板挠度w，可计算扭矩和剪切力如下：

(10)

该问题的边界条件可描述如下：在x=0和x=а两条边上，有挠度w和力矩My为0，即：

(11)

在x=0和x=b两条边上，有挠度w和力矩Mx为0，即：

(12)

根据以上微分方程和边界条件，基于DeepXDE构建相应的PINN求解模型，核心步骤的代码如图 5所示。

图 5 矩形平板承受分布载荷问题的PINN求解模型代码

采用以上代码可以直接求解矩形平板的挠度，结果如图6所示。可以看出，平板四周的挠度为0，越靠近平板中心的挠度越大，最大值达0.004m左右。

图6 PINN求解的矩形平板挠度结果

此外，以上微分方程和边界条件定义的问题存在精确解（S. P. Timoshenko, S. Woinowsky-Krieger, Theory of plates and shells, McGraw-hill, 1959），这里将计算结果与理论结果进行比较。图7为本文PINN求解的弯矩（Mx、Mxy、My）与理论解的对比，图8为本文PINN求解的弯矩（Mx、Mxy、My）与理论解的对比。可以看出，PINN求解结果在定性和定量上均与理论结果吻合。

图7 PINN求解的矩形平板弯矩（Mx、Mxy、My）与理论解比较

图8 PINN求解的矩形平板剪切力（Qx、Qy）和挠度（w）与理论解比较

圆形平板集中受载基准案例

图 9 集中载荷作用下的圆形平板基准案例Vahab M, Haghighat E, Khaleghi M, et al. A physics-informed neural network approach to solution and identification of biharmonic equations of elasticity, 2022.

进一步地，考虑二维圆形平板受集中力的基准案例，如图 9所示。平板的半径和厚度分别为ro=1m和t=0.03 m，四周为刚性支撑，表面则被施加一个集中载荷p=100kN。由于圆形平板的几何存在轴对称性，可以将Eq. (8)转化为极坐标，并在圆周方向进行积分，可获得简化的三阶控制方程如下：

(13)

相应地，扭矩和剪切力的计算公式简化如下：

(14)

该问题的边界条件可描述如下：在r = 0处，转角为0，即：

(15)

在r = ro处，挠度和转角均为0，即：

(16)

根据以上微分方程和边界条件，基于DeepXDE构建相应的PINN求解模型，核心步骤的代码如图10所示。

图 10 圆形平板承受集中载荷问题的PINN求解模型代码

采用以上代码可以直接求解圆形平板的挠度，结果如图11所示。可以看出，平板圆周的挠度为0，越靠近平板中心的挠度越大，圆心处最大达0.04m左右。

图 11 PINN求解的圆形平板挠度结果

图12给出了本文PINN求解的圆形平板挠度与理论结果。可以看出，本文采用PINN求解的挠度与理论解吻合较好。

图 12 PINN求解的圆形平板挠度（w）与理论解比较

根据求解的挠度w，可以进一步计算出给圆形平板的弯矩和剪切力，如图13所示。与扰度类似，弯矩和剪切力也在圆心处取得最大值。

图 13 PINN求解的圆形平板弯矩（Mr）与剪切力（Qr）

04 总结

本期结合三个典型的结构领域案例，重点介绍了飞桨高阶自动微分机制在结构领域问题上的求解能力，这是深度学习在求解由高阶微分方程所描述物理问题上的最新尝试，进一步验证了AI for Science可应用于广泛的工程和科研领域。结合当前提供的案例，用户也可通过扩充时空维度和融合更通用的微分方程，实现三维复杂结构问题的求解。相关案例正在准备中，敬请期待～

引用
[1] 飞桨全量支持业内AI科学计算工具——DeepXDE!https://mp.weixin.qq.com/s/yBsuLWozN4AJCFO5grJ8WA
[2] DeepXDE介绍文档
https://deepxde.readthedocs.io/en/latest/
[3] DeepXDE结构领域科学计算案例https://aistudio.baidu.com/aistudio/projectdetail/5678395

拓展阅读
[1] AI+Science系列（三）：赛桨PaddleScience底层核心框架技术创新详解
[2] 飞桨科学计算实训示例
https://aistudio.baidu.com/aistudio/projectoverview/public?to...
[3] 飞桨黑客松第四期任务—科学计算专题https://github.com/PaddlePaddle/Paddle/issues/51281#science

相关地址
[1]飞桨AI for Science共创计划https://www.paddlepaddle.org.cn/science
[2]飞桨PPSIG Science小组https://www.paddlepaddle.org.cn/specialgroupdetail?id=9