基于iFEM 和虛實結合的加筋板應變場構建研究

陳國材，汪雪良，2，楊華偉，2，蔣鎮濤，2，張 濤，2，張 正

（1.中國船舶科學研究中心，江蘇無錫 214082；2.深海技術科學太湖實驗室，江蘇無錫 214082；3.武漢理工大學，武漢 430070）

0 引 言

船舶在長期服役過程中，遭受風、浪、流多種復雜載荷的作用，極易發生結構疲勞損傷，從而導致船體安全性能下降，通過船體監測系統能夠感知船舶結構自身安全，實現對船舶結構健康狀態的把握[1]。應變作為結構最常采集的物理量之一，同時也是評價船體結構健康狀態的重要依據，通常以實際的離散應變測點為基礎，開展對結構的全局應變場構建。目前，國內外比較成熟的應變場重構方法有模態法[2-4]、Ko 等人提出的位移理論[5-6]、人工智能[7]和iFEM 法[8-9]等。模態法的核心思想是將結構變形視為各階模態的線性組合，其重構精度受模態分析精度影響較大；Ko 的位移理論是基于撓度曲線推導的，所以只適用于重構單方向結構變形[10]；人工智能的方法主要基于人工神經網絡模型，變形重構模型與被測物的材料屬性、結構參數以及環境載荷的分布形式關聯性弱，具有較強的通用性，但其重構位移的精確性嚴重依賴于訓練網絡時所使用載荷集的豐富性[11]；iFEM 通過構建實測應變值與理論值之間的最小二乘誤差函數，結合形函數矩陣與幾何矩陣，組裝成結構全局應變場，iFEM 能擺脫傳統重構方法的局限性，忽略結構的材料屬性和受載荷情況，具有較高的工程應用前景。自2003年Tessler[12]正式提出iFEM 的概念和理論以來，國內外學者針對板、梁、殼等結構開展了深入研究，并在機翼、船體等結構健康監測，復合材料變形場重構等領域進一步驗證了iFEM的魯棒性和可行性[13-14]。

重構方法的精度通常正相關于應變傳感器布置數量，但是實際工程應用中無法通過布置大量傳感器來提升應變場的重構精度，主要有以下原因：一是特殊位置不滿足工程施工條件，如結構角隅處等；二是傳感器數量的增加，將造成系統成本的上升和可靠性的下降。

面對重構精度提升與降低傳感器數量的矛盾，本文提出以虛擬應變傳感器和物理傳感器虛實結合的應變采集方法。所謂虛擬應變傳感器，是指利用數值分析方法給出結構目標位置在不同受載條件下的應變響應，以替代物理應變傳感器的技術手段。在iFEM 技術框架下，通過對船舶結構的典型組成單元加筋板布置虛擬應變傳感器進行應變場重構，可突破少量物理測點難以實現應變場重構的困境，為全船結構應變場重構提供有效方法。

1 船舶加筋板應變場重構理論框架

在iFEM 技術框架下的應變場構建流程可概括為“兩支一主”，“兩支”為根據逆殼單元構建單元理論應變和根據物理模型采集的實測應變，“一主”則是構建關于單元節點位移的最小二乘誤差函數，對誤差函數求極值后，求得單元節點位移，再結合單元幾何矩陣組裝總體矩陣得到應變場，簡化求解流程，如圖1所示。

圖1 iFEM求解流程Fig.1 Solution process of iFEM

為了適應不同的結構形式和物理模型，需要構建不同的逆向單元，由于加筋板模型是應變呈線性變化的板殼結構，并考慮剪切變形對結構的影響，故采用iQS4 作為逆殼單元和Mindlin 板理論進行應變場重構[15]。如圖2所示建立一個板單元的局部坐標系（x,y,z），坐標原點位于板單元的中性面形心處。u、v和w分別是該局部坐標系下x、y和z方向的位移，θx、θy和θz分別為繞x、y和z軸方向的轉角。單元厚度為2h，點1、2、3、4分別為單元的4個節點。

圖2 四節點逆殼單元iQS4Fig.2 Four-node inverse shell unit iQS4

圖3 通過iQS4內的應變花測量的離散表面應變Fig.3 Discrete surface strain measured by strain flower in iQS4 element

式中：u0、v0和w0分別表示節點在x,y,z方向上的位移；θx0和θy0表示中線面上一點繞x軸和y軸方向旋轉的角度，其中z方向的位移w沿板厚方向基本不變，z∈[-h,+h]。

根據求解流程及有限元理論，四節點單元的形函數如式（2），單元的每個節點可表示為式（3）：

則四節點單元位移矩陣可表示為

根據線彈性理論，加筋板結構某單元的理論應變可大致分為面內的拉壓應變e(ue)、彎曲應變k(ue)、剪切應變s(ue)的線性組合，分別如式（5）～（7）所示。

另外，可通過形函數的偏導數矩陣建立起理論應變與單元節點位移矩陣的關系，即加筋板結構的表面應變εb可以用四節點逆殼單元形函數的偏導數矩陣Bm、Bk與單元節點位移矩陣ue表示為

同樣，橫向剪切應變εs也可以表示為四節點逆殼單元形函數的偏導數矩陣Bs與單元節點位移矩陣ue的乘積，

形函數偏導矩陣Bm、Bk、Bs的詳細計算過程及原理見文獻[8]。

構建實測應變與理論應變之間的最小二乘誤差函數，然后求誤差函數對單元節點位移的偏導，當誤差函數為0時，取得理論上該單元各節點的位移值，并認為此時的單元各節點位移分布對應此時實測應變狀態下的實際位移分布。

式中，e(ue)、k(ue)、s(ue)是由中性面位移場表示的單元理論拉壓、彎曲、剪切應變；ec、kc、sc是通過仿真試驗或者應變片測量得到表面應變信息推導得到的中性面的拉壓、彎曲、剪切應變；we、wk、ws為無量綱加權系數，與各截面應變相關，控制著理論中性面應變與實際中性面應變的強弱關系。如果通過應變片或者光纖光柵測量的表面應變信息，可以得到所有單元的中性面應變ec、kc和sc，那么加權系數可取we=wk=ws=1，對于缺失實際中性面應變的情況，如上述中的剪切應變，可以通過調整對應的系數來控制結果的精度以及合理性[10，16]。經過對加筋板某板構造及材質的綜合分析，取ws=λ=10-3，其中0<λ≤1，λ為罰參數[16]。

式中，Ae是板單元內的面積，n為單元內傳感器的數量。

對誤差函數φe(ue)求關于位移向量ue的偏導，并令偏導等于零，求得誤差函數φe(ue)的極小值。

整理上式可得

與有限單元法的單元剛度方程類似，式中：ke為單元剛度矩陣，fe為單元載荷矩陣。

矩陣ke由矩陣Bm、Bk、Bs及其對應的加權系數we、wk、ws確定，由于we=wk=1，ws=λ=10-3，具體形式如下：

矩陣fe與矩陣Bm、Bk、Bs及其對應的加權系數we、wk、ws以及應變傳感器數量n、單元板厚h有關，由于we=wk=1，ws=λ=10-3，具體形式如下：

組裝單元矩陣構建整體矩陣，通過坐標變換矩陣，可將單元在局部坐標系下的剛度矩陣方程轉化為在全局坐標系下的剛度矩陣方程，進而將離散結構的剛度矩陣方程進行整合，得到總體結構的總體剛度矩陣方程，此過程類似于有限元思想[14]。

具體求法與形式如下：

式中，K是總體剛度矩陣，與結構離散的單元位置和數量有關，nel是結構離散的單元數量；U是總體結構位移向量；F是總體載荷矩陣，與結構離散的單元中面應變有關；Te為坐標轉換矩陣。而k矩陣是一個對稱矩陣，與應變測量值無關，與應變測量點所在逆向單元的單元節點位置和應變測量點位置有關，再結合單元邊界條件，系數矩陣k將簡化為一個正定矩陣，求逆后可得到單元節點的全局位移U。

在忽略板單元的剪切變形的情況下，利用單元的形函數偏導數矩陣與單元節點位移矩陣，加筋板結構表面應變可表示為

如計入板單元的結構剪切應變，可根據剪切應變與拉壓應變的關系，合理取罰參數λ，得到

同樣可以經過對離散應變單元的組裝后得到全局應變ε。

2 模型分析與試驗

2.1 加筋板幾何參數和材料屬性

加筋板作為船體結構中的主要構件，在保證船體結構強度和剛度的情況下可大幅減輕鋼材用量。本文中加筋板模型采用低溫船用鋼設計，板面以下為“九縱四橫”結構，即9 根縱骨，4 根橫梁，板面以上布置肋板（如圖4所示）。詳細尺寸及材料屬性見表1。

表1 加筋板幾何參數和材料屬性Tab.1 Geometric parameters and material properties of stiffened plates

圖4 加筋板模型Fig.4 Stiffened plate model

2.2 數值仿真

加筋板以Abaqus 有限元軟件進行建模，板、肋板、固定邊采用殼單元建模，縱骨橫梁采用梁單元建模，共有23 518個單元、71 002個節點，單元均為四邊形單元，即為逆殼單元iQS4（如圖5所示）。

圖5 加筋板有限元模型Fig.5 Finite element model of stiffened plate

為了模擬船舶甲板發生中垂時船舶甲板的載荷響應，試驗時采用一端固定，一端加載，固定端邊界條件為U1=U2=U3=UR1=UR2=UR3=0，加載端邊界條件為U2=U3=UR1=UR2=UR3=0，如圖6所示。

圖6 加筋板加載示意圖Fig.6 Schematic diagram of stiffened plate loading

2.3 模型試驗研究

加筋板模型充分考慮了船舶上層甲板建筑物及工程實際，在加筋板邊緣布置傳感器，以模擬船舶左右舷甲板列板。試驗共布置6 個單向應變片、3 個電阻傳感器和3 個光纖傳感器，本文內容僅用到6 個單向應變片，E1、A1、E3、A2 距最近的固定邊15 cm，E2、E4位于加筋板長度方向的對稱線上，所有測點均距加筋板邊緣2 cm，加筋板試驗模型如圖7所示。

圖7 加筋板模型試驗Fig.7 Stiffened plate model test

試驗過程中，不可避免地存在應變片粘貼位置不夠精確、粘貼工藝不完全相同等狀況。以應變片粘貼最終位置為準，在有限元模型中選取相對應位置的單元，使數值仿真與試驗保持一致，如圖8 所示。實測點與有限元單元號對應表如表2所示。

表2 實測點與單元號對應表Tab.2 Measuring points and corresponding unit numbers

圖8 加筋板測點布置示意圖Fig.8 Schematic diagram of measuring point arrangement of stiffened plate

試驗過程中，先以100 kN為階梯加載至1400 kN后，再以200 kN為階梯加載至2000 kN，如此加載重復兩次，應變數據基本吻合，工況如表3所示。

表3 加筋板加載工況Tab.3 Loading condition of stiffened plate

由于物理測點有限，基于有限測點難以開展高精度應變場重構工作，針對以上困境開展了Xgboost和虛實結合兩種數據補充方法的研究。

3 實測數據分析與處理

將相同位置應變值的實測結果與仿真結果進行對比，并作誤差曲線，如圖9 所示。E1 點誤差較大，其它測點趨勢相同，6個測點均隨著載荷增大，誤差逐漸減小。E2點實測值和仿真值最為接近，除起始工況外，其它工況誤差值保持在1.5%～6.7%之間。A1 點和E1 點、E2 點和E4 點、A2 點和E3 點兩兩形成對照，實測值均有一定差異。另外，實測應變初始值基本不為零，原因有以下幾點：一是加筋板在加工過程中存在焊縫不均勻，板內存有內應力等；二是加載過程存在載荷分布不均勻，導致對應點應變值偏差；三是施工工藝不穩定或存在零點漂移，當載荷為零時，5 個測點出現非零值；四是環境因素，試驗室中的溫度變化，其它機械振動等因素也會對應變采集系統產生影響[17]。綜合實測數據分析，試驗值和仿真值數據較為接近，此次試驗仿真值具有較高參考價值。

圖9 加筋板試驗、仿真載荷-軸向應變圖Fig.9 Test and simulation load-axial strain diagram of stiffened plate

4 虛實結合思路及驗證研究

由于此次軸向壓縮最大加載為2000 kN，軸向應變值最大，此工況下最接近結構的破壞狀態，且具有較高的工程研究意義，故對該載荷作用下的軸向應變場進行重構研究，對6個測點位置分析后，分為4個反演點、2個驗證點，如表4所示。

表4 測點分類Tab.4 Classification of measuring points

將板面按照物理區域劃分為反演點區域和驗證點區域，兩者沿加筋板板面對角線對稱，反演點區域包含4個反演點，驗證點區域包含2個驗證點，如圖10所示。

圖10 反演點、實測點分布圖Fig.10 Distribution of inversion points and measured points

基于iFEM 框架，用實測點、虛擬點、虛實結合點三種路徑分別對加筋板進行應變場構建，并用E1點、E2 點進行精度驗證，由于iFEM 法需要測點上下表面的三向應變，其它方向應變由仿真中的虛擬應變傳感器提供補充數據，流程如圖11所示。在實測點反演路徑中，僅使用4個反演點進行應變場構建，另外采用Xgboost 的回歸方法對實測點進行補充，驗證并對比三種方法的精度。在虛擬點反演路徑中，使用有限元仿真值進行反演，在虛實結合路徑中，以4 個實測點補充部分虛擬點再進行反演并構建加筋板應變場。

圖11 虛實結合方法驗證路徑Fig.11 Verification path of virtual-real combination method

4.1 15個重構點

以15個測點為輸入（I），iFEM 重構值（O）與輸入的軸向應變誤差曲線如圖12所示，實測、虛擬、虛實結合3種路徑的平均誤差（δ）分別為1.82%、1.09%和1.8%，驗證點E1和E2的相對誤差見表5。

表5 15個重構點相對誤差Tab.5 Relative errors of 15 reconstruction points

圖12 15個點I-O-δ圖Fig.12 I-O-δ of 15 points

依據實測點、虛擬點、虛實結合點3種路徑構建實測應變場和重構應變場，如圖13所示，圖中（a）、（b）、（c）分別為3種路徑構建的實測應變場，（d）、（e）、（f）分別為3種路徑的重構應變場。

圖13 15個點應變云圖及iFEM應變重構云圖Fig.13 Strain nephogram and iFEM strain reconstruction nephogram of 15 points

4.2 21個重構點

以21個測點為輸入，iFEM 重構值與輸入的軸向應變誤差曲線如圖14所示，實測、虛擬、虛實結合3種路徑的平均誤差（δ）分別為0.47%、2.1%和2.11%，驗證點E1和E2的相對誤差見表6。

表6 21個重構點相對誤差Tab.6 Relative errors of 21 reconstruction points

圖14 21個點I-O-δ圖Fig.14 I-O-δ diagram of 21 points

依據實測點、虛擬點、虛實結合點3種路徑構建實測應變場和重構應變場，如圖15所示，圖中（a）、（b）、（c）分別為3種路徑構建的實測應變場，（d）、（e）、（f）分別為3種路徑的重構應變場。

圖15 21個點應變云圖及iFEM應變重構云圖Fig.15 Strain nephogram and iFEM strain reconstruction nephogram of 21 points

5 應變場構建精度影響因素評估

應變場構建精度的影響因素諸多，有可定量分析因素，如輸入點的真實度和數量等，有可定性分析因素，如模型加工工藝、試驗加載方式、數據采集技術、環境因素等[18]。通過改變定量物理因素和完善優化相關的定性因素可以更好地評估應變場重構精度。

5.1 Xgboost回歸方法分析

Xgboost（extreme gradient boosting）方法的原理是建立一個新的決策樹模型，該模型正在學習一個新的函數來擬合上次迭代預測的殘差。將所有二分樹模型的回歸值相加，得到最終的回歸結果（如圖16所示）。本研究中，通過設置二分樹的最大深度來限制收斂準則。

圖16 Xgboost回歸法Fig.16 Xgboost regression method

Xgboost預測殘差的目標函數由損失函數和正則化項組成，正則化項為

式中，l(yi,y?i)為損失函數，Ω(fk)為正則化項，yi為樣本，y?i為yi的預測結果，在本研究中，損失函數為線性函數，如式（28）所示：

在Xgboost 回歸方法中，隨機選取16、21、30、40、41、…、48、49、50 個點應變數據，依次放入200 kN數據庫中，以仿真值為基準，插值21個點的應變。圖17為軸向應變在各次插值中的吻合度曲線，在數據庫中，有47 個點時對應的插值結果與仿真值誤差最小，最小誤差為1.92%；當測點超過47 個點時其誤差再次變大，出現過擬合現象。

圖17 軸向應變吻合度曲線Fig.17 Coincidence cure of axial strain data

在此次加筋板模型中，Xgboost回歸方法以仿真與實測值相互補充建立數據庫，同樣基于此方法，在模型試驗時隨機物理測點為47 個時，插值誤差最小，為1.92%，應變場重構結果最佳。數據量豐盈度與插值精度曲線如圖18所示。

5.2 測點數量分析

對此模型9、15、21個測點的情況分別進行了分析研究。通過虛實結合路徑，以9個測點構建的結果與驗證點E1、E2 的相對誤差分別為-2.68%、-4.50%，由于軸向單元劃分區域大，且模型邊界復雜，重構應變場與實際應變場差距大，本文未做詳細分析。以15個測點構建的結果與驗證點E1、E2誤差分別為2.23%、0.56%。以21個測點構建的結果，與E1、E2驗證點誤差則分別為1.62%、-2.13%。實測和虛擬路徑的重構相對誤差如表7所示。

由表可見，虛實結合路徑中以15個點和21個點為輸入時，其結果與驗證點的誤差均小于3%。根據應變場重構結果圖12和圖14，測點為21個點時通過虛實結合路徑構建的應變場與實測值構建的應變場更為吻合。

6 結 論

本文基于iFEM 框架以船舶典型結構加筋板為研究對象，在其邊緣布置傳感器，融入虛實結合思想并結合Xgboost 回歸方法，通過對加筋板模型的數值仿真和加載試驗，解決了同步提升重構精度與降低物理傳感器數量的工程難點，并對試驗過程中可能造成重構應變場精度的因素進行了分析。通過研究得出以下結論：

（1）Xgboost 回歸方法對于補充實測點具有較高的精確性和實用性，在此模型中，通過預測，當物理測點達到30個時插值結果的平均誤差降到2.25%，物理測點達到47個時平均誤差最低，為1.92%。

（2）通過虛實結合路徑快速補充缺失數據的應變場重構結果操作性強、準確度高，以15個點和21個點為輸入時結果與驗證點的誤差均小于3%。

（3）本文針對加筋板結構提供了一套完整的應變場重構方案，對于相似的板架結構具有一定的參考價值，同時虛實結合路徑的驗證也為工程應用提供了有效支撐。