基于降階模型和數據驅動的動態結構數字孿生方法

王青山, 嚴 波, 陳 巖, 鄧 茂, 蔡源斌

(重慶大學 航空航天學院, 重慶 400044)

0 引 言

2003年Grieves首先提出了數字孿生概念[1],近年來其在各個領域受到廣泛關注并成為研究熱點[2-3]．數字孿生充分利用物理模型、傳感器數據、運行歷史數據等,在虛擬空間中完成映射,從而反映相對應的物理結構的全生命周期過程．數字孿生體跟隨物理結構一同演化,可以實現物理結構狀態實時智能監測和預測等．結構數字孿生方法的研究剛剛起步,其構建方法是重要的基礎,具有非常廣闊的應用前景．

近年來,已有關于結構數字孿生的研究報道．Haag和Anderl[4]給出了結構數字孿生的概念,并建立了一個彎曲梁試驗平臺以驗證數字孿生概念,但沒有給出建立數字孿生的具體過程．Moi等[5]用有限元模擬了一個小型轉向節臂起重機的實時變形,基于應變傳感器數據,利用逆方法建立了用于狀態監測的數字孿生體．Wang等[6]結合數值仿真、傳感器數據、多保真代理模型和可視化技術給出了建立數字孿生的方法,并利用桁架結構物理模型對方法進行了驗證．Lai等[7]建立了起重機吊臂的三維有限元高保真模型和空間梁有限元低保真模型,進而建立了吊臂的多保真模型,構建了其數字孿生體以實現結構安全性監測．上述這些工作都沒有考慮結構的損傷變化和所受載荷的動態慣性效應．Kapteyn等[8-9]提出了一種基于組件的降階模型和機器學習算法構建數據驅動數字孿生的方法,給出了一固定翼無人機機翼數字孿生體的建立過程,考慮了機翼的損傷．Milanoski等[10]針對一個航空復合材料加筋板,結合有限元模擬和布拉格光柵傳感器應變測量數據,建立了該加筋板的數字孿生體,其可預測結構所受靜載荷和損傷位置．但文獻[8-10]的工作沒有考慮結構動力載荷的影響．

Wang等[11]針對一旋轉軸,給出了其簡化系統的空間梁參數有限元模型,提出了一種基于參數靈敏度分析的故障模型更新方案,建立了用于故障診斷的數字孿生模型,但對于復雜結構,用于敏感性分析的參數難以確定,難以實現模型的實時更新．Ganguli和Adhikari[12]以單自由度振動系統為對象,提出了同時描述結構服役長時效應和結構動力響應短時效應兩個時間尺度的數字孿生方法．通過測量結構的無阻尼固有頻率、阻尼系數或有阻尼固有頻率的變化,反映結構質量和剛度的變化．在此基礎上,Chakraborty和Adhikari[13]采用多專家模型學習以獲取描述質量和剛度變化的系統參數的演化,但這種方法難以反映結構局部損傷的位置和大小的變化．

數字孿生涉及物理結構不同狀態下大規模的數值模擬,高效的算法是關鍵．模型降階是一有效途徑．模型降階方法主要有基于物理的模型簡化法、投影法和代理模型法等．其中投影法得到了非常廣泛的應用,包括本征正交分解法[14-15]、縮減基方法[16-19]、平衡截斷[20-21]、模態截斷[22]、Krylov子空間方法[23-25]等．其中,Krylov子空間方法具有算法穩定、容易實現、計算量小且適用于大型系統的優點．

本文針對受動載荷作用的結構,考慮結構服役過程中的損傷及其演化,研究了結構數字孿生體的構建方法．設計制作了考慮不同損傷位置和損傷程度的空間框架結構對方法進行驗證,該方法為動態結構數字孿生體的構建奠定了基礎．

1 數字孿生構建及模型降階方法

1.1 數字孿生構建方法

物理實體結構在服役過程中可能出現各種狀態．可以通過預先分析確定這些狀態,包括承受的外載荷、結構原始無損傷狀態、不同位置的損傷及其演化狀態等．針對這些狀態可以建立一個基于物理的模型庫L={L1,L2,…,LN},其中N表示模型數量,庫中每一個模型為對應于物理實體結構一種可能狀態的數值模型．即使有些物理狀態事先可能無法預料,也可以通過后續不斷向模型庫中添加,進行豐富和更新．在物理結構上特定位置布置傳感器測量物理量Xt,如各個測量點的應變．利用模型庫中的每一個模型計算其在不同載荷作用下的響應,獲得訓練數據,并結合機器學習算法訓練獲得一個數據驅動的模型選擇器T．在物理實體結構服役階段,該模型選擇器根據物理實體上傳感器的實測數據Xt推測模型庫中的哪一個模型能最好地解釋測量數據,從而將該模型用作當前時刻t的最新的數值模型dt．該過程可表達為

dt=T(Xt)∈L．

(1)

基于當前時刻的數值模型,快速給出計算結果,如結構位移、應力、應變．整個時程中均重復相同的過程,即根據各個時刻傳感器測量數據不斷推斷當前物理實體結構的狀態,從模型庫中選擇最匹配的數值模型計算給出結果．通過前述的整個過程,可以構建基于物理實體的數字孿生體,其能夠捕捉物理結構上的變化并跟隨其一同演化,快速甚至實時給出計算結果．

模型庫中的模型包含描述物理實體結構在服役期間可能出現的狀態,這些模型可以用微分方程描述,如有限元模型．一方面,精確地模擬一個完整系統需要很大的計算模型和極大的計算量;另一方面,數字孿生體需要跟隨物理實體結構實時響應和演化,幫助做出診斷和決策;最后,要利用這些模型建立一個包含豐富數據的數據集以訓練機器學習模型,獲得模型選擇器,這些都需要進行大量的模擬計算．因此有必要采用模型降階方法,對所有模型進行降階．

1.2 Krylov子空間模型降階

為了對模型庫中的數值模型進行降階,本文采用Krylov子空間模型降階方法．結構的動力學有限元方程可表示為

(2)

其中M,C,K∈n×n分別為結構有限元模型的質量、阻尼和剛度矩陣,F為等效結點載荷列向量,和分別為結點的位移、速度和加速度．對載荷列向量做變換,將式(2)改寫成如下形式：

(3)

其中B∈n×m是常數矩陣,其描述結構所受載荷的分布,f∈m是結構的輸入．

Krylov子空間降階方法的基本思想是通過投影找到一個低階模型,使其一定數量的“矩”和原模型匹配以近似原高階模型．因此,需要找到一個可以直接應用于二階微分方程(3)的投影．考慮如下投影：

a=Var,V∈n×r,a∈n,ar∈r,

(4)

其中r?n．將式(3)乘以矩陣W∈n×r的轉置,并將式(4)代入,得到一個階數為r的降階模型：

(5)

將上式寫成

(6)

其中

Mr=WTMV,Cr=WTCV,Kr=WTKV,Br=WTB．

(7)

根據文獻[25]的推導可知,如果式(5)中的矩陣V是輸入二階Krylov子空間Kr(-K-1C,-K-1M, -K-1B)的基,且其秩為r,矩陣W的選擇使得矩陣WTKV非奇異,那么降階模型和原模型的前r階矩匹配．因此,通過構造輸入二階Krylov子空間,可以計算出W和V．為了簡便,通常取W=V．最終得到一個階數為r的降階模型,其前r階矩和全階模型匹配．

文獻[25]對Arnoldi算法[26]進行了擴展,給出了計算轉換矩陣V的數值算法,其步驟如下：

2) 循環：i=2,3,…,r

(a) 計算向量

(b) 正交化,循環：j=1,2,…,i-1

(c) 標準化

(d) 返回步驟(a)進行下一次循環直到r．

針對1.1小節模型庫中的所有數值模型,均利用Krylov子空間方法建立其對應的降階模型,用這些降階模型替換庫中的原始模型,最終構成一個降階模型庫L．

2 數據驅動數字孿生

2.1 訓練數據集生成

值得一提的是,對于靜力學問題,結構所受載荷和其位移呈線性關系,且一一對應．但對于動力學問題,當前結構位移與載荷歷程有關．對于動力學問題,根據d’Alembert原理,動力學方程可寫成如下形式：

(8)

對于線彈性系統,方程右端等效載荷和位移呈線性關系,且一一對應,可等效為靜力學問題．因此,在實際計算生成訓練數據時,只要載荷范圍足夠大即可涵蓋包含慣性力產生的等效載荷．實際計算時采用降階模型進行計算．

2.2 基于機器學習的模型選擇器

數字孿生體跟隨物理結構一同演化時,需要根據物理結構上的傳感器測量數據實時推斷物理結構當前的狀態,不斷更新計算模型．因此,為了準確推斷物理結構當前狀態,從模型庫中獲取最匹配的模型,本文采用隨機森林機器學習算法去訓練一個模型選擇器T．

隨機森林機器學習算法[27]是一種基于決策樹的集成學習算法．該算法是以決策樹為基分類器的一個集成學習模型{h(X,Θk),k=1,2,…},其參數集{Θk}是獨立同分布的隨機向量．如圖1所示,首先利用Bootstrap抽樣從原始訓練集中抽取k個樣本,且每個樣本的容量都與原始訓練集一樣;其次對k個樣本分別建立k個決策樹模型得到k種分類結果,最終的分類結果由單個決策樹的輸出結果投票決定．

圖1 隨機森林算法示意圖Fig. 1 The flowchart of the random forest algorithm

2.1小節得到的訓練數據集D將作為隨機森林機器學習的輸入,通過分類訓練獲得模型選擇器T,即建立起從傳感器測點處的模擬數據到計算模型的映射關系．

2.3 數據驅動數字孿生及其演化

T：Xt→Lj．

(9)

通過前述方法可創建基于物理結構的數字孿生體,如圖2所示．下標ti表示當前時刻,Sti表示該時刻的物理結構狀態,Xti表示傳感器測量數據,Lti為當前物理結構的計算模型,而Oti表示該時刻的計算結果,可以包括位移、應變和應力等．數字孿生體跟隨物理結構一同演化,自動更新,實時完成計算．通過傳感器的合理設置,還可以實時感知損傷的出現和演化,幫助決策者對結構的安全性做出實時預測和評估,確定是否需要采取相應措施．

圖2 物理結構與數字孿生體Fig. 2 The physical structure and the digital twin

事實上,根據實際情況和需要,數字孿生體還可以向物理實體結構發送反饋信息,通過驅動器實現主動控制,得到更廣泛的應用．

假設傳感器的采樣頻率為f,可以根據結構的響應頻率范圍設定驅動數字孿生模型的時間間隔Δt,對于前述柔性結構,Δt一般遠大于傳感器數據采樣周期1/f．利用傳感器數據驅動模型選擇器后,立即利用選擇的降階模型計算結構響應,只要數據驅動孿生響應時間小于Δt,即可實現數據孿生體跟隨物理實體結構實時演化．

如果由傳感器數據驅動的當前時刻t+Δt的降階模型為Li,該時刻的響應可以利用前一時刻t的響應,用Newmark方法時間積分得到．根據Newmark法,動力學有限元方程(2)的求解格式為[28]

(10)

對于當前的降階模型Li,可利用如下降階方程求解：

(11)

則降階模型的速度和加速度為

(12)

最后利用映射關系(4)可計算全場量

(13)

進一步可以計算結構的應變和應力等．由于利用了降階模型,式(11)—(13)的計算效率非常高,可實現實時反應．

值得一提的是,上述方法取決于結構的響應頻率和數字孿生體的反應時間．對于受洋流或風荷載作用下的海底纜索、輸電導線、桿塔、風力發電機等響應頻率較低的結構,如風力發電機在風荷載作用下的響應周期達到10 s量級[29],結合模型降階該方法可以實現數字孿生體的實時更新．另一方面,對于模型降階后自由度很小的結構,該方法也適用于響應頻率較高的情況．

3 例子：框架結構數字孿生體構建

以一個空間框架結構為例,制作物理實體模型,設計制作不同損傷程度的桿件,采用前述方法構建其數字孿生體,通過實驗驗證數字孿生體構建方法的正確性．

3.1 框架結構物理實體及測試系統

框架結構如圖3所示．組成結構的每一根圓截面桿,半徑為8 mm,構成結構邊框的每根桿長度為400 mm,斜撐為542 mm．材料為鋁合金6061,其密度、彈性模量、Poisson比分別為2.75 g/cm3,74 GPa,0.33．框架與4根剛性很大的立柱連接,立柱固定在基礎上．為了模擬結構的損傷,除了完好桿件外,設計制作了部分帶損傷的桿件．在桿件上特定位置改變圓桿的直徑來反映損傷,直徑越小表示損傷越大．無損傷桿件的直徑為16 mm,表示損傷程度1和損傷程度2的桿件在損傷位置的直徑分別減小到14 mm和12 mm．

圖3 空間框架結構示意圖Fig. 3 Schematic diagram of the frame structure

假設框架結構在外載荷作用下,在特定的位置容易出現損傷．如圖3所示,在損傷位置1和2附近分別設置兩個應變傳感器1、2和3、4,其測量數據作為模型選擇器的驅動數據．為了驗證該框架結構數字孿生體的準確性,在兩個損傷桿件的損傷位置處各布置一個應變傳感器5和6,如圖3所示．

圖4為該框架結構及其測試系統．為了方便桿件的替換,桿件之間采用連接板連接,連接板材料為Q235鋼．可以將不同損傷程度的桿替換在框架中．載荷通過加載機構施加,加載機構與框架之間安裝一個力傳感器,其型號為JLBU-1-50KG,量程和輸出靈敏度分別為50 kg和1.0 mV/V．在傳感器位置沿桿軸向黏貼電阻應變片,其通過DH8302多通道動態應變儀連接筆記本電腦測量應變,數據采樣頻率為50 Hz．

圖4 框架結構物理實體及測試系統Fig. 4 The physical frame structure and the test system

3.2 數字孿生體構建

按前述方法構建框架結構的數字孿生體．假設施加在結構上的載荷位置不變,結構總共有5種狀態,構成有5個模型的模型庫．如圖5所示,對應于桿件無損傷的模型L1;損傷位置1發生損傷程度1(損傷處桿直徑為14 mm)的模型L2;損傷位置1發生損傷程度2(損傷處桿直徑為12 mm)的模型L3;損傷位置2發生損傷程度1和2對應的模型分別為L4和L5．

圖5 框架結構模型庫Fig. 5 The model library of the framework structure

針對這5個模型,分別建立有限元模型．用空間梁單元B31離散該框架結構,每一個模型均劃分2 760個單元,自由度個數為16 560．分別計算得到各有限元模型的質量、阻尼和剛度矩陣M,C,K．其次,利用Krylov子空間模型降階方法對各有限元模型降階,得到對應的降階模型Li,降階模型的自由度r為100．數值計算表明,該框架結構采用原模型進行一次時間積分的時間為0.071 s,采用降階模型為0.003 6 s,效率顯著提高．值得一提的是,可以采用三維實體單元對結構模型進行離散,獲得更高精度的有限元模型,并對其進行降階．

利用該降階模型庫L,按照2.1小節的方法將計算模型上特征點處的應變模擬值作為訓練數據,得到包含37 510個樣本的數據集D,其中劃分的數據集的數量k取20．選擇其中70%的樣本用作隨機森林機器學習的訓練樣本,訓練得到模型選擇器T．余下30%作為測試集,利用訓練得到的模型選擇器對測試數據進行分類預測,模型預測分類與實際分類的混淆矩陣如圖6所示,可見模型選擇器分類的錯誤率僅為1.32%,具有較高分類精度．

圖6 模型選擇器測試數據混淆矩陣Fig. 6 The confusion matrix of the testing sample

3.3 實驗驗證

對如圖4所示框架結構加載點施加動態載荷,力傳感器和應變傳感器1～6的采樣信號通過動態應變儀輸入筆記本電腦記錄載荷和應變時程．將加載過程中的任一時刻,應變傳感器1～4測量得到的應變輸入數字孿生體,驅動模型選擇器T迅速推斷當前時刻框架結構的狀態,并從模型庫中選擇最匹配當前框架結構狀態的計算模型,進而利用力傳感器測量得到的當前時刻的載荷快速計算結構的位移、應變和應力,下一時刻重復相同的過程,數字孿生體便可以實時輸出整個加載過程中任意時刻框架結構的變形和應力．將各個時刻位于損傷區域處的傳感器5和6的應變實測值與數字孿生體的輸出結果進行比較,以驗證孿生體的精度．

由于有限元模型沒有考慮結構中的連接板和螺栓, 實驗過程中發現螺栓預緊力影響較明顯．如果要考慮螺栓的影響, 在建立有限元模型時需考慮連接板處的細節等．經反復調節螺栓預緊力實驗得到了較好的結果, 圖7為結構完好和兩個位置分別出現損傷時實測值與數字孿生體預測結果的比較．圖中給出了力傳感器測得的載荷時程、傳感器1～4測得的應變時程以及傳感器5和6的應變實測值與數字孿生的預測值．注意傳感器1～4測得的應變是模型選擇器的驅動數據．從圖中可見, 數字孿生體預測的損傷位置的應變時程與實測值比較接近, 表明建立的數字孿生體能夠反映物理結構的損傷狀態．值得一提的是, 數字孿生體給出的損傷位置的應變曲線, 在特定位置出現了波動, 而對應時刻的載荷卻沒有波動．這可能是模型選擇器T未能100%的正確選擇匹配模型, 從而在這些時刻根據傳感器輸入數據錯誤推斷了物理實體狀態, 導致在載荷未發生波動的情況下應變發生了波動．然而, 從圖中曲線可見, 這種波動并未對數字孿生后續的預測結果產生大的影響．圖8給出了當框架結構損傷位置2出現損傷程度2時, 其數字孿生體給出的典型時刻的應變分布．

圖7 結構實測與數字孿生預測結果比較Fig. 7 Comparison of measured values and predicted results by digital twin

圖8 框架結構數字孿生體給出的典型時刻的應變分布Fig. 8 Strain distributions at typical moments predicted by the digital twin of the frame structure

為了反映該框架結構數字孿生體在結構損傷演化過程中的更新,假設損傷位置2在受載過程中從初始完好狀態變成損傷狀態,首先采用有限元模擬得到傳感器1～4的應變,然后將其作為傳感器的實測數據,模擬數字孿生體隨物理結構的演化過程．圖9給出了結構損傷位置2處損傷程度發生兩次變化時,數值孿生體的更新和輸出結果．結構在開始是無損傷完好狀態,在1.8 s時損傷位置2發生了損傷,之后在4.48 s時該處損傷進一步加大．從圖中可見,損傷位置2處應變傳感器6和其相鄰位置的傳感器4的應變時程曲線在這兩個時刻均發生了明顯變化．結果表明,一方面數字孿生體在這兩個時刻實現了實時更新,另一方面可以通過在易出現損傷的位置附近布置傳感器監測損傷的發生．

圖9 位置2損傷演變過程中數字孿生體的更新Fig. 9 Update of the digital twin during evolution of damage at point 2

4 結 論

本文針對受動載荷作用結構,提出了一種基于降階模型和機器學習的動態結構數字孿生構建方法,并利用框架結構對方法進行了驗證．得到如下結論：

1)采用Krylov子空間模型降階方法對結構動力學有限元方程進行降階,大大降低了模型自由度,提高了模型計算效率,可以用于數字孿生體的建模;

2)本文提出的動態結構數字孿生方法,對于響應頻率較低的柔性結構,以及響應頻率較高且模型降階后自由度較小的情況,均可實現數字孿生體隨物理結構的實時更新;

3)通過傳感器的合理布置,建立的數字孿生體可以預測結構的損傷及其演化．

本文針對動態結構提出了數字孿生構建方法的框架,尚需進一步對不確定性因素和傳感器采樣信號噪聲的影響等進行深入系統的研究．