面向船體結構動冰載荷監測的線性形函數識別方法研究

孔 帥， 田于逵， 崔洪宇， 季順迎

(1. 中國船舶科學研究中心 水動力學科研部，江蘇 無錫 214082; 2. 大連理工大學 工業裝備結構分析國家重點實驗室，遼寧 大連 116024)

極地豐富自然資源、便利水運交通和深刻地緣影響是刺激海洋強國研發建造冰區船的內在動力[1]，俄美加日韓及北歐諸國在冰區船研發領域均具有極強的話語權，其中俄羅斯擁有世界上數量最多、噸位最大且推進性能最好的破冰船隊伍，目前仍在建造多艘重型破冰船以鞏固其北極地位。《中國的北極政策》[2]指出“技術裝備是認知、利用和保護北極的基礎”，然而我國破冰船在總體船型、推進系統、冰帶構件優化等方面均存在較大的技術短板[3]。冰載荷作為冰區船設計及安全評估中的關鍵輸入始終是兼具科學及工程意義的研究熱點，其可直接用于船型優化[4]、結構安全評估[5]和疲勞分析[6]等。

實船監測是認知掌握冰載荷的重要技術手段，可有效揭示冰載荷的數值范圍[7]、分布特征[8]及統計特性[9]，其載荷識別算法作為監測系統核心近年來被廣泛關注[10-11]。冰載荷識別中常用的影響系數矩陣法通過向待監測子區域依次施加單位載荷以形成關聯起冰載荷矢量和冰激應變響應矢量的傳遞矩陣[12]，具有形式簡潔、工程應用性強的優點，然而該方法尚未考慮冰載荷動載荷效應且不能解決加載位置偏離預期加載位置而導致求解異常的問題。另外，基于時域反卷積算法的動冰載荷識別模型中遞推連鎖計算格式會導致其求解矩陣維數過大和求解系數過多，影響求解穩定性和快速性[13]，進而導致結構安全狀態誤判和報警時機貽誤。

動載荷識別的形函數法的核心思想為有限元分析理論中場變量的離散化處理，利用一組形函數和權重系數擬合動態荷載，把動態載荷時程的求解轉換為個數有限權重系數的求解，縮減求解矩陣維數并提升載荷識別性能。該方法可有效識別移動載荷[14]、均布載荷[15]及集中載荷，分析識別算例可知其具有載荷識別類型范圍廣和適用性良好的優點。

為對冰區船冰載荷監測系統提供有效的載荷識別方法，本文將首先基于形函數方法建立冰載荷識別的正問題，其次討論對反分析求解中遇到的不適定性問題并采用正則化方法控制求解，最后分別從數值算例和試驗驗證角度對提出的冰載荷識別模型進行驗證性分析。

1 冰載荷識別的形函數理論

在連續時域內，船體結構由冰載荷引起的冰激響應可采用杜哈梅爾積分的形式表示[16-17]，即

(1)

式中:g(t)為對應響應與沖擊載荷關系的Green核函數；p(τ)為冰載荷；ε(t)為冰激應變。

時域反卷積識別模型在反演冰載荷過程中，其識別系統對應反問題的適定性主要由式(1)中的小奇異值和測試信號噪聲水平兩個因素決定。實船冰載荷測量中冰激響應信號始終存在噪聲信號的干擾，其可采用低通濾波以濾掉高頻噪聲或傳感器合理布置以提升信噪比的方式控制。因此，決定冰載荷求解算法魯棒性的關鍵在于控制反求矩陣的奇異性，而式(1)表示的反卷積型式因采樣時間間隔均較小而導致其規模較大，進而導致矩陣不可避免的較大奇異性。

形函數方法可將荷載時間歷程的識別轉變為識別數目有限的形函數權重的求解，進而極大地降低系數矩陣的維數，提升計算效率和識別算法的魯棒性。如圖1所示將冰載荷p(t)在時間域內離散成m個等時間間隔的時間元Δt，在每個時間元內再選n個時間節點。

圖1 基于線性分布形函數的冰載荷識別正問題Fig.1 Forward problem of ice load identification based on linear distribution shape function method

(2)

αi=[αi1αi2…αin]T

(3)

φi=[φi1(t)φi2(t) …φin(t)]T

(4)

Aiai=pi

(5)

(6)

求解Aiai=pi后可得權重系數矩陣ai，則

(7)

式中，Ni(t)為形函數，即

(8)

若將第i個時間元分為兩個時間節點(即n=2)，則滿足

[Ni1(t)Ni2(t)]

(9)

此時，Ni1(t)和Ni2(t)可采用線性分布形函數，其定義為

(10)

(11)

由有限元形函數理論可知式(10)和式(11)滿足式(12)和式(13)所示的兩條性質，其中式(12)表示在任意時間元內其形函數之和始終等于1，式(13)表示形函數在自身節點j處時的取值為1，而在其他時間節點時(即k≠j時)為零。這兩條性質可有助于檢驗形函數構造的合理性和保障節點載荷與單元內近似冰載荷之間的連續性。

(12)

(13)

(14)

則可將式(14)寫成矩陣形式

ε=Sp

(15)

船體外板結構監測區域的冰壓識別問題對應多源載荷識別問題，其需要將整個監測區域劃分為眾多子區域。船體結構在多源載荷作用下的響應是各個監測子區域上載荷引起響應的線性疊加，因此可將多源載荷識別問題可寫成

(16)

式中:M為需要監測子區域的數目；N為測點的數目；pi為每個監測子區域上的冰載荷時程；εi為測點上的應變時程；Sij為i監測子區域上冰載荷與j測點應變之間的形函數矩陣。

由此，建立了基于線性分布形函數理論的冰載荷識別模型的正問題。

2 共軛梯度最小二乘正則化算法

已知冰激應變響應ε和形函數矩陣S對冰壓p進行求解的反分析過程通常無法直接用最小二乘法或Moore-Penrose逆法，載荷識別時也通常要引入正則化方法[19-21]。適用于求解大維數、非對稱和非正定方程的共軛梯度最小二乘算法(conjugate gradient least squares, CGLS)是一種高效算法[22]，可較好地識別殼結構和鋸齒結構所受的環境載荷。相比于較為常用的截斷奇異值分解(truncated singular value decomposition, TSVD)等直接型正則化算法，CGLS算法是一種迭代型算法，其每次迭代過程中的搜索方向是兩兩共軛的，且其搜索方向是負梯度方向與上一步搜索方向的組合，而截斷奇異值分解方法則需利用正則化參數優選方法對其小奇異值進行截斷以保證求解的適定性。因而，CGLS方法與TSVD方法的相似之處在于兩者均是在已知邊界約束情況下，根據與原不適定求解問題構造相“鄰近”的適定問題以得到原問題的近似解，而區別在于求解步數、正則化算子及其能力提升方法等方面。

CGLS算法首先要對矩陣S進行正規化處理和相關計算量的初始化

STSpk=STp

(17)

p1=0,r1=εerr,q1=STr1

(18)

式中：p1為共軛梯度算法的初始解;r1=εerr為解的殘差的初始值;q1為共軛迭代過程中間量的初始值。

當k>1時，共軛梯度最小二乘迭代算法過程為

(19)

pk=pk-1+αkqk-1

(20)

rk=rk-1-αkSqk-1

(21)

(22)

qk=STrk+βkqk-1

(23)

式中:αk，βk，qk為共軛迭代過程中間量;rk為迭代k次之后解的殘差。

CGLS正則化算法屬于半收斂算法，不合適的迭代步數會導致“過擬合”或“欠擬合”。CGLS正則化算法在進行冰載荷識別時濾波因子在不同求解時間步及迭代步數時均不同，故需要恰當的終止迭代準則以提升算法的正則化能力。本文采用Paige等[23-24]針對最小二乘問題提出的終止迭代法則(式(24))，該終止準則可較為有效判別當前迭代步數是否為最佳迭代終止點

(24)

式中:γ為正常數；rk為解的殘差。

結合形函數理論與共軛梯度最小二乘正則化算法，其總結的冰壓識別技術流程圖如圖2所示，分為冰載荷識別模型構建、結構冰激應變獲取和冰壓穩定識別算法三個部分，其邏輯及遞進關系見圖2。

圖2 基于形函數理論及CGLS算法冰壓識別的流程圖Fig.2 Flow chart of the ice pressure identification using the shape function theory and CGLS algorithm

3 冰載荷識別模型驗證性分析

至此，基于形函數方法和共軛梯度最小二乘正則化算法的冰載荷識別模型建立，該算法的有效性可通過數值分析及試驗驗證進行檢驗。

3.1 數值分析

為從數值角度檢驗冰載荷識別算法，本文建立的結構模型主要參考中國極地科考破冰船“雪龍2號”外板結構和冰載荷承載力研究工作中的典型破冰船外板結構，該外板由肋骨、強肋骨、縱桁和外板組成，具體尺寸信息標記于圖3中。在冰載荷識別過程中，船體結構假設處于線彈性階段。結構阻尼選擇為比例阻尼。彈性模量為206 GPa，泊松比為0.3，結構鋼材密度為7 850 kg/m3。邊界條件參考ABS冰級船規范標準，在模型的上下平面內邊的邊界設置為在x方向自由，左右平面內邊的邊界設置為關于yoz面對稱。有限元模型中間部分采用網格尺寸為50 mm×50 mm較為精細的網格，為提升計算效率，其余部分采用了尺寸為150 mm×150 mm網格。冰載荷監測區域選擇中間區域，對應監測面積為2.8 m×3 m；肋骨冰激應變監測是目前主流技術手段，故應變測量選擇靠近肋骨邊緣且測試方向平行于外板的應變。

圖3 典型冰區船舶外板結構有限元模型Fig.3 FE model of the typical shell structure of ice-going vessel

實船冰載荷測量中，其冰激響應信號始終存在噪聲的干擾，噪聲信號會對識別數值產生較大的干擾。為分析噪聲信號對冰載荷識別模型干擾程度，本文采用加性噪聲的形式進行添加，其實際結構測量的信號可寫作[25-26]

εerr=ε+lnosiestd(ε)rand(-1,1)

(25)

式中:εerr為含有噪聲信號的應變信號；lnosie為一個百分數，表示噪聲水平；std(ε)為ε標準差；rand(-1,1)為(-1,1)之間的隨機數。

圖 4(a)與圖4(c)為在監測點1基于形函數方法識別的冰載荷與施加載荷時程對比，選擇兩種波形的載荷時程，其采樣時間元Δt選擇為0.01 s。同時，為模擬冰載荷時程的隨機特性，在兩個載荷時程在建立時考慮了加載速率、波形特征及峰值差異等影響。分析外板結構沖擊試驗結果可知平行外板方向布放的傳感器靈敏度最高，最能有效反應載荷信號特征，主要原因是船側肋骨在甲板、縱桁約束下組成的梁結構在沖擊載荷作用下發生“整體彎曲變形”效應，相比于“局部擠壓”、“整體剪切變形”等效應更為突出。圖4(b)為工況一(t=0.25 s時)監測區域的應變云圖，其與真實模型試驗的響應特征是相互對應的，即靠近肋骨邊緣應變(平行外板方向)受“整體彎曲效應”影響較大，驗證了模擬響應的可信性。可以看出無噪聲干擾時載荷識別結果與施加載荷時程之間吻合性較好。參考實船冰載荷測試經驗[28-29]，冰激應變的噪聲水平普遍分布在1%～10%。為評估識別算法在極高噪聲水平干擾下的適用性，圖4也增加了20%高噪聲水平干擾時的載荷識別結果。可以看出，當有20%噪聲干擾時識別結果會因其識別模型的連鎖計算格式受到影響，但仍能較為準確描繪出施加載荷的時程特征和峰值特征，識別相對誤差為13.31%。冰載荷監測中的識別算法的求解速度決定了實時監測的報警效率，冰載荷識別模型由型號為Intel Core i7-9700的CPU處理器運行，運行時間為0.36 s；而采用時域反卷積算法的運行時間則為7.46 s，形函數法在求解速度上具有較強的優勢，實時性更加突出。

圖4 形函數法的冰載荷識別結果Fig.4 Ice load identification results using the shape function method

為研究測點位置及測點數目對載荷識別結果的影響，在圖3模型選擇了三個監測點，表1列出了為研究影響因素所采用的監測對比方案說明。

表1 監測識別影響因素評估方案Tab.1 Evaluation scheme of influence factors in the identification

當測點位置選擇距監測區域0.8 m的監測點3進行識別時，20%高噪聲水平干擾下載荷識別相對誤差為15.59%，這主要緣于測點位置的偏遠會加重形函數矩陣奇異性，而運行時間與作用點幾乎一致。當選擇測點1和測點2時，其多源識別問題中的形函數矩陣規模是單源識別問題中形函數矩陣的N2(N=2為測點數量)倍，矩陣奇異性會增強；導致在20%高噪聲水平干擾下多源載荷識別平均相對誤差為16.41%，運行時間也延長至0.48 s。

3.2 試驗驗證

現場測量中很難布置大型加載裝置以進行加載測試，故基于含有肋骨構件的外板結構載荷識別試驗是檢驗模型準確性較為直接有效的方式。本文采用的外板結構由外板、肋骨、強肋骨和縱桁等組成，如圖5所示。模型采用3 mm厚的6061鋁合金板，試驗模型的邊界采用螺栓和邊條組合方式進行約束。模型的彈性模量、材料密度和泊松比分別設為70 GPa,2 700 kg/m3和0.3。圖5右側為對應的測量位置，測量位置的應變片粘貼方式如圖5中的右上角方框內的示意圖所示。應變片的線路連接采用帶溫度補償片的半橋連接方式。動載荷識別試驗中應變信號和力信號同時由DH5922數據采集系統記錄，采樣頻率10 kHz。為保證肋骨與真實結構之間的結構相似性最大化，肋骨縱向端部與縱桁相固定，且監測及測量的位置均遠離邊界。

圖5 載荷識別試驗裝置Fig.5 Experimental facility in the load identification

圖5左上角外板方形區域為對應的載荷測試區域，試驗過程中采用力錘施加載荷，其靈敏度為3.10 mV/N。利用力錘連續向圖5所示載荷測試方形區域施加四次力，圖6中的實線部分即為對應的施加載荷時程，圖6中虛線為識別的載荷時程。冰激應變信號預處理首先圍繞沖擊力作用時刻進行截取，然后采用去零漂及抗混濾波技術進行信號處理，轉錄接近真實的冰激應變信號，操作步驟如圖2右上角所示。然后，按圖2所示整體技術路線圖將冰壓進行識別，獲得測試區域的施加力的識別值。可以看出，因結構自身阻尼和振動特性影響，識別載荷時程會稍遲滯于施加載荷時程并且在脈沖信號尾部出現波動現象。然而，識別載荷仍能較好地描述施加載荷時程特征。另外，實船測試中冰載荷時程具有隨機特性且識別結果通常會受到周邊結構上作用載荷的干擾，故各國學者常選取冰載荷峰值作為統計分析的參數，本算例中識別載荷峰值與施加載荷峰值之間相對誤差均值僅為3.21%。

圖6 施加載荷與識別載荷對比Fig.6 Comparison between the applied loads and identified loads

綜合數值和試驗結果可知，基于形函數的冰載荷識別算法具有求解速度較快、穩定性良好及精度較高的優點，適用于船體外板結構動冰載荷識別分析。

4 結 論

實船冰載荷監測工作在國內尚屬起步狀態，其監測識別算法作為核心技術可直接提升監測的整體水準，進而提高我國冰區艦船的設計建造能力。基于以上研究內容可得出以下結論及研究方向：

(1)本文結合有限元形函數理論、共軛梯度正則化算法和最小二乘問題的終止迭代法則等建立了船體結構動冰載荷快速識別模型，并利用數值和試驗手段對模型識別效果進行了檢驗。結果表明基于形函數載荷識別算法相比于常用的影響系數矩陣法，補充考慮了冰載荷的動力效應；同時相比于Green核函數冰載荷識別算法，求解快速性得到了保障。

(2)本文研究證明了線性形函數方法在描述識別冰載荷時程特征方面具有一定的可行性，而實際中的冰載荷可能會表現出極強的隨機特性和離散特性，因此需要更高階形函數作為基函數以更加準確表征冰載荷時程，提升識別模型普適性。同時，為提升試驗評估的全面性，將在冰載荷識別試驗增加考慮多源載荷識別、真實海冰試樣加載及低溫環境。

(3)現有實船冰載荷測試數據均是利用影響系數矩陣法識別的，國內外均未有學者采用動載荷識別算法對其原有認知的冰載荷時空分布特征、峰值分布擬合情況及峰值參數化表征函數等規律進行對比研究。下一步可利用實測數據進行分析，明確原有冰載荷識別算法的缺陷和修正冰載荷譜等統計規律。