付書山， 孫廣開， 何彥霖， 初大平

（1.北京信息科技大學光電測試技術(shù)及儀器教育部重點實驗室，北京 100192；2.北京信息科技大學光纖傳感與系統(tǒng)北京實驗室，北京 100016）

變形監(jiān)測（Shape sensing）是飛行器結(jié)構(gòu)健康監(jiān)測的重要內(nèi)容，指利用現(xiàn)有傳感技術(shù)實時獲取航空器在飛行狀態(tài)下的各種信息，借助先進的信息處理方法和力學建模方法，計算出飛行器結(jié)構(gòu)的形狀、大小和位置在空間域或者時間域的變化。通過對比結(jié)構(gòu)參數(shù)在狀態(tài)前和狀態(tài)后變化量，可以分析結(jié)構(gòu)的損傷程度以及損傷位置，評估當前健康狀態(tài)、預(yù)測飛行器的有效壽命，以降低系統(tǒng)故障風險、減少結(jié)構(gòu)安全維護成本[1-4]。

飛機蒙皮是包圍在飛機骨架結(jié)構(gòu)之外，形成飛機氣動力外形的維形構(gòu)件。其中機翼蒙皮與機翼骨架結(jié)構(gòu)作為產(chǎn)生和傳遞氣動載荷的主要部件，本身承受有較大承載力。機翼蒙皮的形狀直接決定了飛行器的氣動特性，是影響飛機安全性能的關(guān)鍵部件，加上機翼蒙皮直接與外界接觸，受力情況較為復(fù)雜，所以對機翼蒙皮變形測量技術(shù)是飛機結(jié)構(gòu)健康監(jiān)測的技術(shù)支撐，也是飛機安全監(jiān)控的重要組成部分。隨著機翼尺寸逐漸增大，結(jié)構(gòu)越來越復(fù)雜，對飛機的安全性和可維護性提出了更高的要求。隨著對飛行器結(jié)構(gòu)變形監(jiān)測的深入研究，基于光纖傳感的結(jié)構(gòu)變形監(jiān)測方法成為該領(lǐng)域研究的核心內(nèi)容[5-7]。光纖布拉格光柵（Fiber bragg grating，F(xiàn)BG）傳感器是以光信號為傳感方式的新型傳感器，與傳統(tǒng)傳感器相比，除了具有體積小、質(zhì)量輕、靈敏度高、抗電磁干擾、耐腐蝕能力強等優(yōu)點外，還可通過波分復(fù)用的方式串聯(lián)組成傳感網(wǎng)絡(luò)[8-9]，其解調(diào)系統(tǒng)采集頻率高、性能穩(wěn)定，特別適用于航空航天動態(tài)監(jiān)測領(lǐng)域。

目前，國內(nèi)外一些高校和研究機構(gòu)利用光纖傳感技術(shù)和曲率信息對結(jié)構(gòu)變形監(jiān)測展開了研究。Clements 等[10]則提出了利用FBG傳感器實現(xiàn)精確三維定位和形狀測量系統(tǒng)，通過布置在一個由多芯光纖構(gòu)成的細小智能電纜的FBG傳感器測量電纜上的各個測點的彎曲曲率和撓率，推導(dǎo)獲得電纜上的各個點的相對位置和方位，從而重構(gòu)整條智能電纜的形態(tài)；Jutte等[11]提出利用FBG傳感器陣列，在機翼地面載荷試驗中進行全尺度的機翼彎曲和扭曲形態(tài)實時測量研究；Glaser等[12]利用FBG測量的簡易梁結(jié)構(gòu)曲率數(shù)據(jù)，進行了結(jié)構(gòu)形態(tài)實時重構(gòu)的試驗分析和研究；Sun等[13]采用光纖Bragg光柵采集了模擬機翼表面的48處應(yīng)變數(shù)據(jù)，使用曲面擬合的方法對機翼的變形場進行了重構(gòu)，并使用視覺測量的方法對重構(gòu)結(jié)果進行評價。Yi等[14]提出一種基于曲率信息的太陽帆板的變形重構(gòu)算法，通過在帆板結(jié)構(gòu)中布置光纖光柵傳感器，結(jié)合可視化技術(shù)，實現(xiàn)了結(jié)構(gòu)的空間曲面擬合和動態(tài)變形重構(gòu)。張合生等[15]采用了正交FBG網(wǎng)絡(luò)獲取二維正交曲率信息，提出針對空間曲面的變形重構(gòu)算法，并與利用單向曲率信息的重構(gòu)算法進行了比較。結(jié)果表明其在彎曲變形和扭轉(zhuǎn)變形狀態(tài)下均有更高的重構(gòu)精度。

上述學者提出的變形重構(gòu)方法都較好地滿足了特定領(lǐng)域變形監(jiān)測需求，但還存在一些不足。基于曲率信息的重構(gòu)方法適用于大變形的場景中，在處理微變形的情況時，易出現(xiàn)重構(gòu)曲率不連續(xù)導(dǎo)致重構(gòu)的曲面不夠光滑和連續(xù)，重構(gòu)精度不理想。針對此問題本研究提出一種基于逆有限元法和應(yīng)變信息的變形重構(gòu)算法，該算法采用三節(jié)點三角形單元進行結(jié)構(gòu)離散，并結(jié)合Kirchhoff板殼彎曲理論進行力學建模，理論推導(dǎo)了逆有限元算法的可行性，最后使用翼形蒙皮結(jié)構(gòu)和方形蒙皮結(jié)構(gòu)驗證了變形重構(gòu)算法的可行性和通用性。

1 算法原理

有限元法（Finite element method，F(xiàn)EM）是求解數(shù)理方程的數(shù)值計算方法，是將一個連續(xù)體的求解區(qū)域離散為有限個形狀簡單的單元（子區(qū)域），各單元通過有限的節(jié)點進行連接，載荷通過節(jié)點在單元之間傳遞。根據(jù)節(jié)點的平衡條件建立節(jié)點的載荷場方程，然后將所有節(jié)點組合進行綜合求解，以獲得復(fù)雜工程的近似數(shù)值解。逆有限元法（inverse finite element method，iFEM）與有限元法思想相同，其本質(zhì)為利用分片插值逼近全域問題[16-17]，但在某些環(huán)節(jié)的構(gòu)造與計算上存在區(qū)別。首先根據(jù)待求解對象的幾何特點，將其視為板、桿、梁、殼一種或多種結(jié)構(gòu)的組合體，然后選擇有限個一維、二維或三維單元對其進行離散。單元內(nèi)任一點的位移d（x）可以由沿坐標軸方向的分位移表示，即

式中，x≡（x,y,z）表示單元上任一點的坐標；u、v、w分別為沿x軸、y軸、z軸方向的位移分量；N（x）為單元形函數(shù)；qe為單元全部節(jié)點位移組成的位移向量。

針對結(jié)構(gòu)的幾何特點和受力特點，選擇合適的受力模型，如Kirchhoff薄板理論、Mindlin中厚板理論等來描述結(jié)構(gòu)在變形狀態(tài)下應(yīng)變和位移的關(guān)系。結(jié)構(gòu)體的應(yīng)變場可以由n個獨立的應(yīng)變分量來描述，每個應(yīng)變分量表示一個方向的應(yīng)變，即ε≡{ε1,ε2,ε3, ...,εn},將應(yīng)變場用單元位移向量qe和單元的插值函數(shù)來表示為

式中，x≡（x,y,z），表示應(yīng)變測點的坐標；[]為求導(dǎo)算子，表示對單元形函數(shù)N（x） 求導(dǎo)；S（x）為應(yīng)變-位移矩陣。

通過理論應(yīng)變值ε與實測應(yīng)變值εε構(gòu)建的誤差函數(shù)Φ，其中εε≡{}，使用最小二乘法求解誤差函數(shù)的極值。誤差函數(shù)表達式為

式中，k為整體結(jié)構(gòu)設(shè)置的應(yīng)變測點數(shù)量；假設(shè)結(jié)構(gòu)體離散為M個單元，每個單元設(shè)置m個應(yīng)變測點，則總的誤差Φ可以看作所有單元誤差Φe之和：

對于每個單元，當Φe取值最小時，可以認為此時的模型最能反映出結(jié)構(gòu)的變形情況。將式（2）帶入式（4）可知，單元誤差函數(shù)Φe是關(guān)于單元位移向量qe的函數(shù)，根據(jù)最小二乘法的求解原理，當Φe取得最小值時，此時有

求解方程會得到Keqe=Fe形式的方程，其中Ke為單元剛度方程；Fe為單元載荷向量。通過相應(yīng)的坐標變換，對各個單元剛度矩陣、單元載荷向量進行組裝，再添加結(jié)構(gòu)體的邊界限制條件，最終得到結(jié)構(gòu)體的總體平衡方程為

式中，K為結(jié)構(gòu)總體剛度矩陣，只與單元形函數(shù)和測點坐標有關(guān)；F為總體載荷向量，只與測點應(yīng)變值有關(guān)；q為待求節(jié)點位移向量，通過計算q=K-1F即可求出節(jié)點位移。至此，通過結(jié)構(gòu)表面的應(yīng)變值反解出結(jié)構(gòu)位移值的全部過程，通過實時讀取應(yīng)變值并計算q=K-1F即可得到結(jié)構(gòu)體各個節(jié)點的實時位移，對于非節(jié)點處的位移，可以通過式（1）中節(jié)點位移和形函數(shù)插值得到。機翼蒙皮變形重構(gòu)逆有限元算法實現(xiàn)流程如圖1所示。

圖1 機翼蒙皮變形重構(gòu)逆有限元算法實現(xiàn)流程Fig.1 Inverse finite element algorithm for wing skin deformation reconstruction

2 試驗方案

2.1 機翼蒙皮變形特點

分析機翼蒙皮的幾何尺寸和受力特點，由于其厚度遠小于長寬尺寸，如圖2所示，可以將其看作薄板結(jié)構(gòu)。機翼蒙皮承受空氣動力作用后將作用力傳遞到相連的機身機翼骨架上，機翼蒙皮的位移為機翼的位移，作用于機翼蒙皮上的載荷分為兩種類型： （1）平行于蒙皮表面的縱向載荷，屬于平面應(yīng)力問題； （2）垂直于蒙皮表面的橫向載荷，屬于薄板的彎曲問題。由于機翼的實際受力主要為垂直于機翼表面豎直方向的氣動載荷，機翼蒙皮的受力變形問題可以采用Kirchhoff彎曲理論來分析，該理論滿足3個假設(shè)。

圖2 機翼蒙皮結(jié)構(gòu)模型Fig.2 Wing skin structure model

（1）直法線假設(shè)。薄板中面法線變形后仍保持為法線且長度不變，有

（2）板厚不變假設(shè)，忽略板中面的法線應(yīng)力分量，且不計其引起的應(yīng)變。板內(nèi)各點的撓度w與z坐標無關(guān)，只是x，y的函數(shù)，即

（3）中面不變形假設(shè)，即薄板中面內(nèi)的各點沒有平行于中面的位移，即

按照Kirchhoff計算假設(shè)，機翼蒙皮結(jié)構(gòu)的變形彎曲問題等價于蒙皮中面的變形彎曲問題，中面上的一點實際上代表著一個長度為板厚2t法線段，t表示中面到薄板表面的距離，即中面的變形撓度就是蒙皮整體的變形撓度。對于中面上的任意一點位移都有3個位移分量w、θx、θy,分別表示撓度、法線繞x的轉(zhuǎn)角和法線繞y的轉(zhuǎn)角。當中面的位移為小撓度位移時，其轉(zhuǎn)角位移分量θx、θy與撓度w存在如下關(guān)系：

機翼蒙皮表面的應(yīng)變場ε包含3個應(yīng)變分量{εx,εy,εxy}，分別表示x方向的應(yīng)變、y方向應(yīng)變、xy面內(nèi)的切向應(yīng)變。根據(jù)幾何方程，應(yīng)變場ε可表示為

由式（11）可知，中面（z=0）內(nèi)既沒有面內(nèi)應(yīng)變，也沒有面內(nèi)應(yīng)力，只需測量蒙皮表面（z=±t）的應(yīng)變數(shù)據(jù)，根據(jù)應(yīng)變位移的幾何關(guān)系，可求解機翼蒙皮的變形量。

2.2 結(jié)構(gòu)離散方案

逆有限元法求解蒙皮彎曲變形問題，首先需要對機翼蒙皮結(jié)構(gòu)進行結(jié)構(gòu)離散。常用的離散單元有三角形和矩形單元，由于三角形單元有較好的邊界適應(yīng)性和靈活性，工程實用價值較大。本研究選取9自由度三角形單元來對蒙皮結(jié)構(gòu)進行離散，如圖3所示，三角形單元節(jié)點位移向量qe可表示為

圖3 9自由度三角形單元Fig.3 Nine degree of freedom triangle element

2.3 應(yīng)變測點布設(shè)方案

由于蒙皮的上下表面關(guān)于中面對稱，當蒙皮發(fā)生彎曲變形時，其上下表面的應(yīng)變數(shù)值理論上是相等的，通過在蒙皮表面粘貼應(yīng)變花，獲取該測點3個方向的應(yīng)變值εx、εy、γxy。單個應(yīng)變花由3個FBG傳感器構(gòu)成，其放置方向與x軸分別成0°、45°和90°夾角，如圖4所示，其中各方向FBG測量值與應(yīng)變分量的對應(yīng)關(guān)系為

為了使應(yīng)變測點能夠較為均勻地分布在蒙皮表面，對進行結(jié)構(gòu)離散的每個三角形單元均設(shè)置7個應(yīng)變測點，測點位置如圖4所示，測點的坐標值在表1中列出。其中測點1在三角形單元的形心位置，測點2、3、4分別在形心到3個頂點的三等分處，測點5、6、7分別在線段2—3、3—4、2—4的中點。通過選取不同位置、不同數(shù)量的應(yīng)變測點來探究重構(gòu)精度與測點位置和數(shù)量之間的關(guān)系。

圖4 單元應(yīng)變測點布設(shè)位置和FBG布設(shè)方式Fig.4 Location of strain measuring points and FBG layout in unit

表1 應(yīng)變測點位置坐標Table 1 Position coordinates of strain measuring points

2.4 機翼蒙皮變形有限元仿真

為驗證逆有限元算法的變形重構(gòu)性能，選擇基于Ansys平臺對機翼蒙皮結(jié)構(gòu)進行有限元變形仿真。蒙皮材料選擇鋁合金（楊氏模量E=7.1e+4MPa,泊松比v=0.33），幾何尺寸如圖5（a）所示，厚度為2mm。采用1843個三角形單元對其進行離散，機翼左端固定，整體為懸臂板受力模型，并在蒙皮表面施加50N的均布載荷,方向為z軸正方向（圖5（b））。機翼蒙皮變形有限元仿真云圖如圖5（c）～（f）所示，選取蒙皮表面的仿真應(yīng)變值作為試驗測量的應(yīng)變值，帶入iFEM形變重構(gòu)算法，把重構(gòu)出的位移量與有限元仿真得到的位移量進行對比，可以得到算法的重構(gòu)精度。

圖5 機翼蒙皮變形有限元仿真結(jié)果Fig.5 Finite element simulation results of wing skin deformation

2.5 算法重構(gòu)精度評價

使用FEM的分析結(jié)果作為對照真值，對iFEM重構(gòu)算法的性能進行評估。評價指標包括各個節(jié)點的重構(gòu)誤差ei以及各節(jié)點重構(gòu)誤差的均方根值（RMSE），其表達式為

式中，ei為節(jié)點重構(gòu)誤差的表達式；為iFEM算法重構(gòu)的節(jié)點位移；為有限元仿真得到的節(jié)點位移；m為節(jié)點總個數(shù)。

3 結(jié)果與討論

3.1 測點密度α與重構(gòu)誤差的關(guān)系

為探究單元內(nèi)應(yīng)變測點密度與重構(gòu)精度的關(guān)系，分別將機翼蒙皮結(jié)構(gòu)離散為4、6、8個單元，如圖6（a）～（c）所示，并定義測點密度變量α，表示每個單元內(nèi)設(shè)置α個應(yīng)變傳感器，α取值為1～7。當測點密度α=1時，該測點默認在三角形的形心處。當測點密度α＞1時，則在圖4所示的位置中隨機選擇測點。3種離散方案對應(yīng)的測點密度與重構(gòu)誤差的關(guān)系如圖6（d）～（f）所示，對于每種離散方案，機翼蒙皮變形重構(gòu)誤差均隨單元測點密度增大而減小，隨后趨于平穩(wěn)。在單元測點密度α從2增至7的過程，3種離散方案的重構(gòu)誤差分別減少了1.08%、0.79%和0.45%，減小幅度并不明顯。當單元測點密度相同時，重構(gòu)誤差隨單元數(shù)量增加而減小，具體數(shù)據(jù)如表2所示，仿真位移對比如圖7所示。另外，只在形心位置設(shè)置一個測點時（α=1），其重構(gòu)誤差相對更小。

圖7 離散方案3、α=1工況下，iFEM算法重構(gòu)位移與FEM仿真位移對比Fig.7 Comparison of displacement reconstructed by iFEM and simulated displacement of FEM in discrete scheme 3, α=1

表2 離散方案3、α=1工況下，各節(jié)點位移重構(gòu)誤差Table 2 Displacement reconstruction error of each node in discrete scheme 3, α=1

圖6 機翼蒙皮離散方案和測點密度與重構(gòu)誤差的關(guān)系Fig.6 Discrete scheme of wing skin, and relationship between density of measuring point and reconstruction error

3.2 測點密度β與重構(gòu)誤差的關(guān)系

為進一步探究逆有限元算法的應(yīng)變測點密度與重構(gòu)精度之間一般性的關(guān)系，并驗證該算法對蒙皮變形重構(gòu)具有通用性，繼續(xù)選取具有較大尺寸蒙皮結(jié)構(gòu)作為研究對象。大尺寸蒙皮結(jié)構(gòu)是構(gòu)成并維持飛行器氣動外形的重要組件，廣泛應(yīng)用于空間站、衛(wèi)星帆板、深空探測器、大展翼飛機等航空航天器結(jié)構(gòu)組成中[18]，其尺寸通常在數(shù)米以上。本研究選取邊長為6m×6m的方形蒙皮結(jié)構(gòu)，使用Ansys平臺對其進行仿真變形，固支條件為4邊固定，并設(shè)置垂直蒙皮表面大小為20N的均布載荷。蒙皮使用材料和厚度與前文中翼型蒙皮相同，以保持蒙皮受力彎曲特性的一致性。定義測點密度變量β，表示每m2設(shè)置β個應(yīng)變傳感器，如圖8所示。

圖8 6種不同的應(yīng)變測點密度布設(shè)方案Fig.8 Six different layout schemes of strain measuring point density

通過在蒙皮表面設(shè)置不同的測點密度，獲得測點密度與重構(gòu)精度的關(guān)系，如圖9所示，總體重構(gòu)誤差隨測點密度增大而減小。當測點密度β=1.39個/m2時，重構(gòu)誤差＜2%，當測點密度β＞3個/m2后，重構(gòu)誤差基本不再減小并趨于平穩(wěn)。當測點密度β從5.56個/m2減小至2個/m2時，重構(gòu)誤差從最小值0.51%上升至1.2%，誤差未明顯增大，仿真位移如圖10所示。

圖9 蒙皮結(jié)構(gòu)的測點密度與重構(gòu)誤差的關(guān)系Fig.9 Relationship between density of measuring points and reconstruction error of plate structure

圖10 測點密度β=2個/m2工況下，iFEM算法重構(gòu)位移與FEM仿真位移對比Fig.10 Comparison of displacement reconstructed by iFEM and simulated displacement of FEM in measuring point density β=2/m2

4 結(jié)論

（1）采用翼形蒙皮結(jié)構(gòu)和方形蒙皮結(jié)構(gòu)對逆有限元算法進行驗證，試驗結(jié)果表明，逆有限元變形重構(gòu)算法對于蒙皮變形監(jiān)測具有可行性和通用性。

（2）變形重構(gòu)誤差隨離散單元數(shù)量和測點密度增加而減小；測點密度設(shè)置為2個/m2時，整體重構(gòu)均方根誤差為1.2%，最大位移處的重構(gòu)誤差為1.07%，重構(gòu)時間＜20ms；在測點密度一定時，測點位置放置在單元形心處可以獲得更小的重構(gòu)誤差。

（3）基于應(yīng)變信息的逆有限元變形重構(gòu)方法可以實現(xiàn)機翼蒙皮變形動態(tài)監(jiān)測，在飛機蒙皮結(jié)構(gòu)變形監(jiān)測領(lǐng)域具有應(yīng)用前景。采用離散單元分析思想，能夠適用于微小變形和復(fù)雜結(jié)構(gòu)變形監(jiān)測領(lǐng)域，應(yīng)用場景廣泛。