在工程和科學(xué)領(lǐng)域，AI與科學(xué)方法的結(jié)合正在解決更多經(jīng)典科學(xué)問題，并在固體、流體、傳熱、材料等越來越多的領(lǐng)域得到了驗證。為了更好支撐科研人員開展AI與基礎(chǔ)學(xué)科的交叉融合研究，百度飛槳不斷完善框架能力，提供支持AI+科學(xué)計算的復(fù)數(shù)算子、高階自動微分機制以及高性能編譯等等。本期我們聚焦結(jié)構(gòu)領(lǐng)域中的“高階”物理機理，向大家說明飛槳如何利用“高階自動微分算子”處理結(jié)構(gòu)領(lǐng)域中的經(jīng)典問題。在案例中，我們將介紹如何使用飛槳框架對4階以上偏微分方程進行無監(jiān)督求解，同時結(jié)合DeepXDE或賽槳PaddleScience等科學(xué)計算工具組件獲得更好的開發(fā)體驗。

01 結(jié)構(gòu)領(lǐng)域背景與痛點

在結(jié)構(gòu)領(lǐng)域中，飛機機翼桁架、汽車底盤、輪船甲板等機械結(jié)構(gòu)件的受力變形、破壞以及疲勞損傷等都是最典型的“力學(xué)”工程難題，解決好這些經(jīng)典工程挑戰(zhàn)無論對基礎(chǔ)科研或工程應(yīng)用都具有十分重要的意義。

傳統(tǒng)的結(jié)構(gòu)分析方法包括有限元法、有限差分法、邊界元法等，通常需要大量的計算資源和人力投入，這會限制模型的規(guī)模和精度。近年來，基于物理信息神經(jīng)網(wǎng)絡(luò)（PINN）的AI方法逐漸被應(yīng)用于結(jié)構(gòu)領(lǐng)域的物理方程求解，這類方法普遍被認(rèn)為兼具提升計算速度與降低人力投入的優(yōu)勢。

不同于經(jīng)典NLP、CV等領(lǐng)域問題，其中多數(shù)任務(wù)都是基于梯度下降算法進行優(yōu)化，只需要一階導(dǎo)數(shù)（即梯度）來更新模型參數(shù)。物理問題通常需要使用高階微分方程（ODEs/PDEs/fPDEs/IDEs）及相應(yīng)的初邊值條件來描述物理機理，這要求融合物理機理的AI方法可處理高階微分。目前深度學(xué)習(xí)框架中高階微分算子往往通過手寫、組合等方式實現(xiàn)，尤其是4階以上微分算子的定義難度及實現(xiàn)工作量都非常大。

02 結(jié)構(gòu)領(lǐng)域問題求解原理

PART1 結(jié)構(gòu)領(lǐng)域問題描述
結(jié)構(gòu)領(lǐng)域問題通?？捎闷胶馕⒎址匠?、幾何方程和物理方程來描述，其通用形式一般很難求解。實際應(yīng)用中多使用簡化方程，如廣泛存在的平面問題，從而降低求解難度。

針對二維線彈性問題，我們可以找到某些標(biāo)量來代表應(yīng)力函數(shù)，從而可獲得同時滿足平衡微分方程和相容方程的標(biāo)量微分方程。在忽略體積力的情況下，二維線彈性問題可由以下偏微分方程描述：

(1)式中，φ是極坐標(biāo)系(r,θ)下的Airy應(yīng)力函數(shù)。在直角坐標(biāo)系中，該偏微分方程可描述如下：? ? ? ? ? ? ? ? ? ?? ? ? ? ??

式中，(x,y)=r(sinθ,cosθ)。以上方程也被稱為雙調(diào)和方程，它通過一個統(tǒng)一的表達式來描述平面線彈性問題。

理論上，基于方程(1)或(2)以及充分的邊界條件，即可實現(xiàn)問題求解。然而，由于方程包含四階雙調(diào)和算子?4，傳統(tǒng)數(shù)值方法需要使用更高階的單元和算法才能進行求解，而其計算時間和成本是極高的，大型復(fù)雜結(jié)構(gòu)問題中則普遍面臨著“維度爆炸”、精度低、求解難等問題。對比傳統(tǒng)求解方法，PINN方法具備自適應(yīng)采樣、非線性激活函數(shù)、并行計算等特點，可以將PDEs嵌入到神經(jīng)網(wǎng)絡(luò)中，可有效地解決高維空間和復(fù)雜幾何形狀的問題。

一旦求解獲得Airy應(yīng)力函數(shù)，則可計算應(yīng)力分量如下：

?采用胡克定律，則可計算應(yīng)變分量如下：

式中，對于平面應(yīng)變和平面應(yīng)力問題，k分別取3-v和(3-v)(1+v)；μ可取為剪切模量G，v為泊松比。通過以上高階偏微分方程，我們可以對結(jié)構(gòu)領(lǐng)域的典型問題進行定義。不同于傳統(tǒng)的數(shù)值差分的計算方法，基于給定的偏微分方程，利用PINN方法可在無監(jiān)督的情況下直接對結(jié)構(gòu)領(lǐng)域問題進行求解，具體如下。

PART2 PINN方法求解結(jié)構(gòu)問題
圍繞結(jié)構(gòu)領(lǐng)域中高階微分方程所描述的工程和科學(xué)問題，飛槳框架可以支持PaddleScience和DeepXDE科學(xué)計算工具組件構(gòu)建相應(yīng)的PINN求解模型，并針對前面描述的二維線彈性結(jié)構(gòu)問題，構(gòu)建了求解應(yīng)力函數(shù)的PINN模型，如圖1所示。

該求解模型引入二維線彈性結(jié)構(gòu)問題的控制方程和邊值條件，這些物理信息約束不僅降低了神經(jīng)網(wǎng)絡(luò)對數(shù)據(jù)的依賴，也降低了對網(wǎng)絡(luò)復(fù)雜度的要求，使得采用一個簡單的全連接網(wǎng)絡(luò)即可實現(xiàn)二維線彈性結(jié)構(gòu)問題的求解。

同時為了解決該領(lǐng)域中高維偏微分方程的定義與求解，飛槳一方面在框架中增加高階微分算子，另一方面構(gòu)建基礎(chǔ)算子體系，可支持不限階數(shù)的自動微分。目前，飛槳已經(jīng)支持絕大多數(shù)算子的無限階自動微分，以及部分算子的有限階自動微分，主要算子簡述如下：

無限階自動微分：包括標(biāo)量與Tensor之間的加、減、乘、除、冪運算，elementwise系列運算、矩陣乘matmul、激活函數(shù)tanh，以及assign、concat、cumsum、expand_v2、reverse、squeeze、unsqueeze、scale、tile、transpose、sign、sum、mean、flip、cast、slice等；

有限階自動微分：sin、cos、sigmoid等算子支持三階計算?；陲w槳最新提供的無限高階自動微分（AD）算子，我們也可在飛槳全量支撐的DeepXDE科學(xué)計算工具中，通過導(dǎo)入paddle.fluid中的core模塊，利用deepxde.gradients.jacobian和deepxde.gradients.hessian可分別計算一階和二階微分，以及通過它們的任意組合來實現(xiàn)高階微分。對此，我們基于飛槳實現(xiàn)了PINN方法求解結(jié)構(gòu)領(lǐng)域中1維梁變形、2維平面受載兩類典型問題求解，并得到與理論解一致的結(jié)果。

03 典型案例介紹

下面根據(jù)前面給出的二維線彈性問題PINN求解模型，介紹基于飛槳聯(lián)合科學(xué)計算工具DeepXDE構(gòu)建的結(jié)構(gòu)領(lǐng)域典型案例，包括歐拉梁撓度計算、矩形平板集中受載基準(zhǔn)案例、圓形平板分布受載基準(zhǔn)案例（具體可見DeepXDE結(jié)構(gòu)領(lǐng)域科學(xué)計算案例）。

PART1 一維歐拉梁問題

歐拉梁（Euler beam）是被用來描述長條形物體彎曲運動的物理模型。在橋梁、飛機機翼、吊橋等眾多情況下，都可以用歐拉梁模型來描述這些物體的彎曲、扭轉(zhuǎn)和振動等運動。歐拉梁模型假設(shè)物體是細長且剛性的，可以在一個平面上自由彎曲，其幾何通?？杀缓喕癁橐痪S。DeepXDE的官方目錄examples / pinn_forward下包含該類案例的一個具體實例Euler_beam.py，其對應(yīng)的控制方程如下：

式中，w為撓度，也即選取的應(yīng)力函數(shù)φ。其邊界條件可描述如下：

• 左端點（x=0）固定，因而其擾度和轉(zhuǎn)角為零，即：

• 右端點（x=1）承載，因而其彎矩和剪切力為零，即：

根據(jù)以上給出的撓度四階方程和邊界條件，基于DeepXDE可以簡單地構(gòu)建相應(yīng)的PINN求解模型，主要步驟的代碼如圖 2所示。

采用以上代碼可以求解歐拉梁的撓度，結(jié)果如圖 3所示。可以看出，梁的左端因為固定而撓度為0，從左到右逐漸增大。

PART2 矩形平板分布受載基準(zhǔn)案例

考慮更復(fù)雜的二維矩形平板受分布載荷的基準(zhǔn)案例，如圖 4所示。平板的長、寬和厚分別為a=2m、b=3m和t=0.01m，厚度為四周被簡單支撐，表面則被施加一個正弦分布載荷?ππ?，相應(yīng)的控制方程可描述如下：

式中，w為平板撓度，D為抗彎剛度，可計算如下：

式中，E為彈性楊氏模量，v為泊松比。根據(jù)平板撓度w，可計算扭矩和剪切力如下：

該問題的邊界條件可描述如下：

• 在x=0和x=а兩條邊上，有撓度w和力矩My為0，即：

• 在x=0和x=b兩條邊上，有撓度w和力矩Mx為0，即：

根據(jù)以上微分方程和邊界條件，基于DeepXDE構(gòu)建相應(yīng)的PINN求解模型，核心步驟的代碼如圖 5所示。

采用以上代碼可以直接求解矩形平板的撓度，結(jié)果如圖6所示?？梢钥闯?，平板四周的撓度為0，越靠近平板中心的撓度越大，最大值達0.004m左右。

此外，以上微分方程和邊界條件定義的問題存在精確解（S. P. Timoshenko, S. Woinowsky-Krieger, Theory of plates and shells, McGraw-hill, 1959），這里將計算結(jié)果與理論結(jié)果進行比較。圖7為本文PINN求解的彎矩（Mx、Mxy、My）與理論解的對比，圖8為本文PINN求解的彎矩（Mx、Mxy、My）與理論解的對比?？梢钥闯?，PINN求解結(jié)果在定性和定量上均與理論結(jié)果吻合。

PART3 圓形平板集中受載基準(zhǔn)案例

進一步地，考慮二維圓形平板受集中力的基準(zhǔn)案例，如圖 9所示。平板的半徑和厚度分別為ro=1m和t=0.03 m，四周為剛性支撐，表面則被施加一個集中載荷p=100kN。由于圓形平板的幾何存在軸對稱性，可以將Eq. (8)轉(zhuǎn)化為極坐標(biāo)，并在圓周方向進行積分，可獲得簡化的三階控制方程如下：

相應(yīng)地，扭矩和剪切力的計算公式簡化如下：

該問題的邊界條件可描述如下：

• 在r = 0處，轉(zhuǎn)角為0，即：

• 在r = ro處，撓度和轉(zhuǎn)角均為0，即：

根據(jù)以上微分方程和邊界條件，基于DeepXDE構(gòu)建相應(yīng)的PINN求解模型，核心步驟的代碼如圖10所示。

采用以上代碼可以直接求解圓形平板的撓度，結(jié)果如圖11所示。可以看出，平板圓周的撓度為0，越靠近平板中心的撓度越大，圓心處最大達0.04m左右。

圖12給出了本文PINN求解的圓形平板撓度與理論結(jié)果?？梢钥闯?，本文采用PINN求解的撓度與理論解吻合較好。

根據(jù)求解的撓度w，可以進一步計算出給圓形平板的彎矩和剪切力，如圖13所示。與擾度類似，彎矩和剪切力也在圓心處取得最大值。

04 總結(jié)

本期結(jié)合三個典型的結(jié)構(gòu)領(lǐng)域案例，重點介紹了飛槳高階自動微分機制在結(jié)構(gòu)領(lǐng)域問題上的求解能力，這是深度學(xué)習(xí)在求解由高階微分方程所描述物理問題上的最新嘗試，進一步驗證了AI for Science可應(yīng)用于廣泛的工程和科研領(lǐng)域。結(jié)合當(dāng)前提供的案例，用戶也可通過擴充時空維度和融合更通用的微分方程，實現(xiàn)三維復(fù)雜結(jié)構(gòu)問題的求解。相關(guān)案例正在準(zhǔn)備中，敬請期待～

引用

[1]?飛槳全量支持業(yè)內(nèi)AI科學(xué)計算工具DeepXDE!https://mp.weixin.qq.com/s/yBsuLWozN4AJCFO5grJ8WA

[2] DeepXDE介紹文檔
https://deepxde.readthedocs.io/en/latest/

[3]?DeepXDE結(jié)構(gòu)領(lǐng)域科學(xué)計算案例https://aistudio.baidu.com/aistudio/projectdetail/5678395

拓展閱讀

[1]?AI+Science系列（三）：賽槳PaddleScience底層核心框架技術(shù)創(chuàng)新詳解

[2] 飛槳科學(xué)計算實訓(xùn)示例
https://aistudio.baidu.com/aistudio/projectoverview/public?topic=15

[3] 飛槳黑客松第四期任務(wù)—科學(xué)計算專題https://github.com/PaddlePaddle/Paddle/issues/51281#science

相關(guān)地址

[1]飛槳AI for Science共創(chuàng)計劃https://www.paddlepaddle.org.cn/science

[2]飛槳PPSIG Science小組https://www.paddlepaddle.org.cn/specialgroupdetail?id=9