鄭勇剛, 吳哲同, 張涵博, 劉振海, 葉宏飛, 張洪武

(大連理工大學 工程力學系 工業(yè)裝備結構分析優(yōu)化與CAE軟件全國重點實驗室 大連 116024)

1 引 言

梁、板和殼等結構力學模型由于其形式簡潔和具有準確表征物體主要變形模式的能力,其在工程中獲得了廣泛應用。在有限元法FEM[1]和高性能計算機尚未普及之前,這些模型的建立為結構的力學響應分析帶來了極大便利。雖然這些模型的控制方程目前可以從彈性力學基本方程退化得來[2],然而其建立及發(fā)展成熟花費了科學家數(shù)十年乃至上百年的時間。同時,隨著科學的不斷發(fā)展,需要解決的力學問題變得愈加復雜,這使得傳統(tǒng)建模方式面臨巨大的挑戰(zhàn)。近年來,數(shù)據(jù)科學作為繼實驗、理論和模擬之后的第四個科學范式,其與計算力學的結合得到了廣泛的關注。結合大數(shù)據(jù)[3]、人工智能[4]和機器學習[5,6]等技術的數(shù)據(jù)驅動計算力學新范式發(fā)展迅猛[7,8],各種機器學習方法廣泛應用于建立材料本構模型[9,10]和高效求解力學響應[11,12]等科學問題以及智能仿真與制造[13]、異常管理監(jiān)測[14,15]和數(shù)字孿生[16]等工程問題。

對于數(shù)據(jù)驅動力學建模問題,其主要難點在于構建快速穩(wěn)健、自動化的控制方程識別方法[7]。目前獲得廣泛應用的控制方程識別方法主要有兩類,即符號回歸方法和稀疏學習方法。符號回歸方法采用二叉樹形式的拓撲結構來描述輸出與輸入變量之間的函數(shù)關系,并對結構進行迭代優(yōu)化,得到具有較小誤差的函數(shù)表達式,從而實現(xiàn)對輸出變量的非線性擬合[17]。近幾年,符號回歸與貝葉斯學習[18]和圖神經(jīng)網(wǎng)絡[19]等技術相結合,在具有強非線性的擬合問題中也展現(xiàn)出強大的能力。稀疏學習方法則假設在整個狀態(tài)搜索空間中,控制方程實際包含的只有其中少量的項,因此約束控制方程的稀疏性不僅可以平衡模型的精度和復雜性,也可以提升控制方程的可解釋性和泛化能力。目前應用最廣泛的稀疏學習方法有LASSO[20]和SINDy[21]。除了上述兩大類識別方法外,深度學習[6]也廣泛用于擬合變量之間隱藏的非線性關系[22,23],如基于物理信息的神經(jīng)網(wǎng)絡PINN[24]通過自動微分技術[25,26]將已知具體形式的控制微分方程的殘差加入到目標函數(shù)中,可從已知數(shù)據(jù)中識別出控制方程的系數(shù)。此外,量綱分析和最小作用量原理也結合用于方程的識別[27],采用此方法已從三維靜力學模擬數(shù)據(jù)中準確識別出了細長梁的Euler-Bernoulli撓曲線方程[28]。

盡管上述方法成功應用于一些力學建模問題,但它們都需要預先提供關于控制方程的結構化先驗知識,如控制方程的個數(shù)及需要進行回歸的左端項表達式,在本質上都屬于有監(jiān)督的數(shù)據(jù)驅動方法。在面對一個嶄新且復雜的力學系統(tǒng)時,上述方法所需的先驗知識往往都是未知的,這使得其應用范圍和可靠性大大降低?；谝环N無監(jiān)督的控制方程識別方法,序列奇異值過濾Seq-SVF[29]算法,本文提出了一種基于彈性力學第一性原理的數(shù)據(jù)驅動力學建模方法。使用此方法從細長結構的高精度有限元計算數(shù)據(jù)中自動識別出了Timoshenko梁形式的撓曲線控制微分方程;區(qū)分了三種不同加載條件下剪切影響系數(shù)關于結構尺寸和力學參數(shù)的函數(shù)表達式,并對經(jīng)典理論公式進行了修正。

2 基于彈性力學第一性原理的數(shù)據(jù)驅動力學建模方法

2.1 Seq-SVF方法介紹

Seq-SVF針對于解決如下形式的控制方程識別問題。

(1)

X·Λ=0

(2)

式中X為通過直接測量及計算派生項得到的數(shù)據(jù)矩陣,Λ為控制方程的系數(shù)陣。由于控制方程個數(shù)及需要進行回歸的左端項都是未知的,因此LASSO和SINDy等經(jīng)典的系統(tǒng)識別方法不再適用。Seq-SVF基于奇異值分解方法SVD能夠提取矩陣本征線性結構的特點,設計了一套高度自動化的迭代流程,用于準確識別方程個數(shù)及從候選函數(shù)庫中提取方程中實際存在的項。此方法分為預處理、奇異值過濾、識別左端項和篩選相關項四個步驟。

2.1.1 預處理

由式(1,2)可知,每個控制方程的系數(shù)對應的列向量都是方程X·x=0的一個非平凡解,且系數(shù)列向量之間是線性無關的。因此,數(shù)據(jù)矩陣的零空間的維數(shù)等于控制方程的個數(shù)P?；谏鲜鲈?則可使用SVD識別控制方程個數(shù)。SVD將一個矩陣分解為三個矩陣的乘積,即

X=UΣVT

(3)

式中U和V均為正交矩陣,Σ為奇異值矩陣,其存放了原始矩陣的所有奇異值,即

(4)

當矩陣X奇異時,只有前rank(X)個奇異值是非零的,即σrank(X)+1=σrank(X)+2=…=σN=0。由于數(shù)據(jù)矩陣中存在噪聲,使得這些嚴格為0的奇異值變?yōu)榉橇?其數(shù)量級與其他奇異值數(shù)量級的相對大小取決于噪聲的水平。因此,在信號原有噪聲相對較小并且對數(shù)據(jù)進行去噪處理后,可以通過SVD確定相對較小的奇異值個數(shù),即為控制方程的個數(shù)P。Seq-SVF選用高斯濾波方法對數(shù)據(jù)集進行去噪處理。此外,為了避免某項因為數(shù)量級太小而在后續(xù)識別過程中忽略,對數(shù)據(jù)矩陣的每一列使用L2正則化處理,

(i=1,…,M;j=1,…,N)

(5)

2.1.2 奇異值過濾

由于在候選函數(shù)庫中的函數(shù)數(shù)量可能遠大于控制方程中實際存在的項數(shù),因此需要構造一個自動化的迭代流程用于將無關項過濾掉。在獲得控制方程個數(shù)P的基礎上,Seq-SVF依次刪除數(shù)據(jù)集矩陣中的各列,進一步實施SVD并識別缺失此列的矩陣的控制方程個數(shù)來實現(xiàn)奇異值的過濾,即如果仍然識別出P個方程,說明此列對應的項不包含在方程組里面;如果識別出P-1個方程,說明這項在方程組里面。將有關項對應的列重組成數(shù)據(jù)矩陣,確保了所識別方程的稀疏性,能夠提升結果的可解釋性和泛化能力。

2.1.3 識別左端項

在過濾掉與控制方程無關的項后,數(shù)據(jù)矩陣的零空間中任何P個線性無關的向量理論上都可以作為方程組的系數(shù)。然而如果不加約束地隨意選取,所識別出的每一個方程都很有可能只是最簡潔形式方程的復雜線性組合,這為理解其中蘊含的物理意義和進一步的應用增加了難度。Seq-SVF的策略是選擇最為線性無關的P列對應的項作為方程的左端項,其原理是由于待定的左端項可以由其余的項線性表示,因此保持左端項之間的獨立性可以確保其與其余項具有強線性相關性。在線性代數(shù)領域中,上述問題稱為列子集選擇問題CSSP[30]或內插分解ID[31],Seq-SVF采用了求解此問題的代表性方法sRRQR[32],只需輸入數(shù)據(jù)矩陣以及待提取列的個數(shù)P,即可識別出對應左端項在原數(shù)據(jù)矩陣中的位置。將識別出的左端項重新排列在數(shù)據(jù)矩陣的左側,記為Y1-YP;其余項為右端項,記為YP+1-YN′。

2.1.4 篩選相關項

為了確定各個方程的右端項及其系數(shù),依次對每個控制方程實施2.1.2節(jié)所述的步驟。即對P個左端項做循環(huán),依次將每個左端項與所有的右端項組成數(shù)據(jù)矩陣,此時通過SVD會識別出1個控制方程,之后依次刪除數(shù)據(jù)矩陣中的各右端項并實施SVD:如果仍然識別出了1個控制方程,說明此右端項不包含在當前循環(huán)的左端項對應的控制方程里面;如果識別出了0個控制方程,說明此項應當保留。通過以上步驟,即可獲得所有控制方程的所有項,其形式為

(6)

最后使用最小二乘法計算各項系數(shù),得到最終的控制方程。需指出的是,對于在2.1.2節(jié)中過濾掉的項,其在式(6)中對應的系數(shù)直接設置為零。

綜上所述,Seq-SVF作為一種無監(jiān)督的方程識別方法,能夠在不依賴先驗知識的情況下,自動識別出控制方程的個數(shù)并為每個方程識別一個左端項。借助2.1.2節(jié)和2.1.4節(jié)提出的識別方程稀疏結構方法,提升了結果的泛化能力和可解釋性,相比于LASSO和SINDy結果更為準確并且具有更廣的適用范圍。數(shù)值分析表明,采用Seq-SVF從文獻和數(shù)值模擬數(shù)據(jù)中可成功識別出多種經(jīng)典的線性和非線性系統(tǒng)的控制方程,圖1展示了對Lorenz動力系統(tǒng)的識別結果。

圖1 基于Seq-SVF的Lorenz動力系統(tǒng)控制方程識別結果

2.2 基于彈性力學第一性原理的數(shù)據(jù)驅動力學建模方法流程

基于上述提出的Seq-SVF方法,本文進一步提出了基于彈性力學第一性原理的數(shù)據(jù)驅動力學建模方法,流程如圖2所示,其可用于構建梁、板和殼等結構力學模型。本方法的主要步驟如下,首先對幾何結構進行建模,施加載荷,使用有限元求解力學響應;之后收集感興趣的位置(如梁中軸線和板中截面)處數(shù)據(jù),根據(jù)具體需求設置候選函數(shù)庫的形式;最后使用Seq-SVF自動識別其控制方程。

圖2 基于彈性力學第一性原理的數(shù)據(jù)驅動力學建模方法流程

3 梁彎曲模型的數(shù)據(jù)驅動識別

3.1 梁彎曲模型的建立與簡化

考慮圖3所示的受橫向載荷作用的細長結構,為了從中軸線的橫向位移有限元計算數(shù)據(jù)中建立梁撓曲線微分方程,首先根據(jù)模型的數(shù)學性質對控制方程的形式進行大體的猜測,從而對模型進行簡化,以減少后續(xù)流程的工作量。

圖3 橫向受載梁

控制方程(組)中可能包含坐標x、撓度w(x)和載荷q(x)三個變量,以及結構的尺寸和力學參數(shù),即楊氏模量E、泊松比ν、長l、寬h和高b,因此控制方程(組)可表示為

f(w(x),q(x),x,E,ν,l,b,h)=0

(7)

使用有限元計算得到的撓曲線數(shù)據(jù)具有如下特征,(1) 由于不考慮材料和幾何非線性,因此撓曲線受線性規(guī)律控制,即控制方程是線性的,不含x·w,(?w/?x)2和sin(x)等非線性項; (2) 在不受任何載荷時有w≡0,因此控制方程是齊次的,即方程中不含有常數(shù)項和僅與x有關的函數(shù)項; (3) 由于要保證控制方程對于坐標的平移不變性,即容許剛體位移(滑移不變性),因此方程中不含有w這一項。

綜上所述,控制方程形式上是一個(組)齊次線性常微分方程,

(8)

式中i為控制方程的序號(i=1,2,…)。對于微分階次的選取,參考Euler-Bernoulli梁和Timoshenko梁的控制方程,對撓曲線的最高微分階次設置為4,對載荷的最高微分階次設置為2。因此待定的控制方程形式為

(9)

3.2 載荷施加在梁上表面時撓曲線方程識別結果

首先,針對載荷施加在梁上表面的工況進行研究。此種工況可以用于分析梁承放重物的情況,并且許多求解矩形梁的彈性力學位移解的工作都是假設載荷施加于梁上表面的。根據(jù)式(9),候選函數(shù)庫中包含7種函數(shù),分別是撓度對坐標的1階～4階微分和載荷對坐標的0階～2階微分。對于載荷的選取,應當滿足兩個條件以確保識別結果的正確性,一是q,?q/?x和?2q/?x2三者線性無關,否則將會識別出額外的控制方程;二是三者應具有相近的數(shù)量級,以確保其中某項不會因為太小而忽略。因此,將載荷設置為q=1/(sinx+2),施加在梁的上表面,使用有限元法進行計算并收集中軸線上所有節(jié)點的撓度,其中單元類型為平面應力四節(jié)點單元,厚度為1。通過數(shù)值微分計算出7個候選函數(shù)的值并構建數(shù)據(jù)矩陣,將其代入Seq-SVF中識別控制方程。在步驟(1)預處理中,Seq-SVF識別出數(shù)據(jù)矩陣具有一個明顯小的奇異值,說明撓曲線受一個方程控制。之后在步驟(2)中依次刪除各列,發(fā)現(xiàn)刪除?w/?x,?2w/?x2,?3w/?x3和?q/?x四項后,控制方程的個數(shù)依舊是一個,說明這些項不存在于控制方程中;而刪除?4w/?x4,q和?2q/?x2三項后控制方程的個數(shù)變?yōu)榱?說明控制方程只包含這三項。最后,Seq-SVF通過最小二乘法計算了各項的系數(shù)

α(?4w/?x4)-q+β(?2q/?x2)=0

(10)

并且發(fā)現(xiàn)?4w/?x4的系數(shù)α與抗彎剛度EI的大小相等,即

EI(?4w/?x4)=q-β(?2q/?x2)

(11)

這與Timoshenko梁的控制方程形式一致[33,34]。對于?2q/?x2的系數(shù),通過改變模型參數(shù),發(fā)現(xiàn)β與彈性模量E和梁長度l無關,隨泊松比ν和梁高度h的改變而改變。為了進一步驗證結果,將載荷設置為多項式形式,即

q=x3+x2+x

(12)

對撓曲線進行多項式擬合,得到撓曲線的四階微分項的函數(shù)式,即

(13)

將式(12,13)代入式(11),等式左右相等,再次驗證了所識別彎曲方程的準確性。對比兩種載荷下β的值,發(fā)現(xiàn)幾乎一致,列入表1,說明控制方程的系數(shù)不受具體載荷形式的影響。此外,通過表1中β隨ν和h的變化情況,發(fā)現(xiàn)β與ν有線性關系,并且與h的平方成正比,即

β=h2(m+n·ν)

(14)

通過對表1的數(shù)據(jù)進行線性擬合,得到β的函數(shù)表達式為

(15)

表1 β與泊松比ν和梁高度h的定量關系

3.3 載荷施加在梁內部及中軸線上時撓曲線方程識別結果

針對載荷施加在梁內部及中軸線上的工況進行研究。其中,載荷分布在內部可以準確分析梁的自重對于梁彎曲變形的影響;當發(fā)生火災時,受溫度的影響,混凝土與其內部的鋼筋具有不同的強度折減系數(shù),并且有時鋼筋會發(fā)生失效甚至熔化,導致結構的坍塌,此時通過在中軸線上施加等效熱應力載荷可以準確分析梁整體的力學響應[35]。研究方法與3.2節(jié)相同,得到β與泊松比ν和梁高度h的定量關系,列入表2。兩種工況下β仍然與ν有線性關系,并且與h的平方成正比。通過線性擬合得到

(16)

表2 載荷施加在梁內部和中軸線上時β與 泊松比ν和梁高度h的定量關系

3.4 與經(jīng)典理論對比

通過將式(11)的系數(shù)與Timoshenko梁模型[30,31]進行對比,并將I=bh3/12,G=E/[2(1+ν)]和A=bh代入,得到剪切影響系數(shù)與載荷二階微分項系數(shù)β的函數(shù)關系為

κ=h2(1+ν)/(6β)

(17)

將式(15,16)代入式(17),得到三種加載位置對應的κ值為

(18)

通過與已有理論對比發(fā)現(xiàn),載荷施加在梁上端時,識別出的控制方程對應Olsson[36]提出的表達式;載荷分布在梁內部時,κ=5/6,與經(jīng)典理論一致;載荷施加在中軸線上時,沒有已有理論與之對應,即通過數(shù)據(jù)驅動方法發(fā)現(xiàn)了一個新的理論公式用于分析此工況。

為了驗證識別結果的準確性,設計了一些算例,將經(jīng)典Euler-Bernoulli梁和Timoshenko梁的解、本文識別方程的數(shù)值求解結果,以及有限元高精度解進行對比。長和高為4×0.4的兩端簡支梁在上端受載荷q=1/(sin3x+2)的結果如圖4所示。對于數(shù)值求解本文識別方程時需設置的邊界條件,分別嘗試了Euler-Bernoulli梁和Timoshenko梁形式的邊界條件。結果表明,采用本文識別出的控制微分方程,并將邊界條件設置為Timoshenko梁的形式時,其數(shù)值求解結果與有限元計算結果最為接近。

圖4 與有限元解和經(jīng)典模型解對比

之后將梁尺寸改為10×1,載荷施加在上表面和中軸線上的結果如圖5所示,再次驗證了識別結果的準確性。

圖5 與有限元解和經(jīng)典模型解對比(尺寸:10×1)

最后,設計了一個三明治型結構,來說明本文識別的梁模型對于分析復雜復合工況的準確性。載荷如圖6(a)所示,在梁的上下表面和中軸線處均受到加載,并且整體受到重力的作用。對于傳統(tǒng)梁模型,是不區(qū)分載荷位置是在上線和中線還是內部的,而本文模型能夠區(qū)分不同加載位置得到的不同撓曲線響應。通過線性疊加原理對撓曲線進行計算

(19)

與有限元解對比如圖6(b)所示,可以看出兩者吻合良好。此外,通過變形后梁的形狀可以發(fā)現(xiàn),梁的截面在變形后依舊保持近似平面,而且截面與中軸線是不垂直的,這在x=2和x=8處的截面體現(xiàn)得尤為明顯。這說明了Timoshenko梁模型保留截面變形后為平面的假設并拋棄軸線與截面垂直的假設,與彈性力學解非常相符。

圖6 復合加載工況下與有限元解的對比

4 結 論

本文提出了一種基于彈性力學第一性原理的數(shù)據(jù)驅動力學建模方法,采用無監(jiān)督的控制方程識別方法Seq-SVF,能夠從彈性力學高精度解中自動識別簡化的力學模型。使用精細有限元計算梁受橫向載荷下的撓曲線數(shù)據(jù),識別出的控制微分方程中包含載荷的二階微分項,其形式與Timoshenko梁的控制方程一致。通過在同一尺寸的模型中施加具有不同形狀的載荷,得出控制方程系數(shù)不隨載荷的形式改變的結論。通過改變模型參數(shù),得到載荷二階微分項的系數(shù)與模型參數(shù)的定量關系,并對三種載荷加載位置下剪切修正因子的函數(shù)表達式進行了區(qū)分,其中載荷施加在上表面和內部時識別出的控制方程分別與兩種已有理論對應,在載荷施加在梁中軸的工況下發(fā)現(xiàn)了一種新的剪切修正因子表達式。最后對識別的控制方程進行數(shù)值求解,得到的結果比經(jīng)典Timoshenko梁方程的數(shù)值解要更為接近彈性力學精確解。