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

(1. 中國船舶科學(xué)研究中心 水動力學(xué)科研部，江蘇 無錫 214082; 2. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024)

極地豐富自然資源、便利水運(yùn)交通和深刻地緣影響是刺激海洋強(qiáng)國研發(fā)建造冰區(qū)船的內(nèi)在動力[1]，俄美加日韓及北歐諸國在冰區(qū)船研發(fā)領(lǐng)域均具有極強(qiáng)的話語權(quán)，其中俄羅斯擁有世界上數(shù)量最多、噸位最大且推進(jìn)性能最好的破冰船隊(duì)伍，目前仍在建造多艘重型破冰船以鞏固其北極地位?！吨袊谋睒O政策》[2]指出“技術(shù)裝備是認(rèn)知、利用和保護(hù)北極的基礎(chǔ)”，然而我國破冰船在總體船型、推進(jìn)系統(tǒng)、冰帶構(gòu)件優(yōu)化等方面均存在較大的技術(shù)短板[3]。冰載荷作為冰區(qū)船設(shè)計(jì)及安全評估中的關(guān)鍵輸入始終是兼具科學(xué)及工程意義的研究熱點(diǎn)，其可直接用于船型優(yōu)化[4]、結(jié)構(gòu)安全評估[5]和疲勞分析[6]等。

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

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

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

1 冰載荷識別的形函數(shù)理論

在連續(xù)時域內(nèi)，船體結(jié)構(gòu)由冰載荷引起的冰激響應(yīng)可采用杜哈梅爾積分的形式表示[16-17]，即

(1)

式中:g(t)為對應(yīng)響應(yīng)與沖擊載荷關(guān)系的Green核函數(shù)；p(τ)為冰載荷；ε(t)為冰激應(yīng)變。

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

形函數(shù)方法可將荷載時間歷程的識別轉(zhuǎn)變?yōu)樽R別數(shù)目有限的形函數(shù)權(quán)重的求解，進(jìn)而極大地降低系數(shù)矩陣的維數(shù)，提升計(jì)算效率和識別算法的魯棒性。如圖1所示將冰載荷p(t)在時間域內(nèi)離散成m個等時間間隔的時間元Δt，在每個時間元內(nèi)再選n個時間節(jié)點(diǎn)。

圖1 基于線性分布形函數(shù)的冰載荷識別正問題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后可得權(quán)重系數(shù)矩陣ai，則

(7)

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

(8)

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

[Ni1(t)Ni2(t)]

(9)

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

(10)

(11)

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

(12)

(13)

(14)

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

ε=Sp

(15)

船體外板結(jié)構(gòu)監(jiān)測區(qū)域的冰壓識別問題對應(yīng)多源載荷識別問題，其需要將整個監(jiān)測區(qū)域劃分為眾多子區(qū)域。船體結(jié)構(gòu)在多源載荷作用下的響應(yīng)是各個監(jiān)測子區(qū)域上載荷引起響應(yīng)的線性疊加，因此可將多源載荷識別問題可寫成

(16)

式中:M為需要監(jiān)測子區(qū)域的數(shù)目；N為測點(diǎn)的數(shù)目；pi為每個監(jiān)測子區(qū)域上的冰載荷時程；εi為測點(diǎn)上的應(yīng)變時程；Sij為i監(jiān)測子區(qū)域上冰載荷與j測點(diǎn)應(yīng)變之間的形函數(shù)矩陣。

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

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

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

CGLS算法首先要對矩陣S進(jìn)行正規(guī)化處理和相關(guān)計(jì)算量的初始化

STSpk=STp

(17)

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

(18)

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

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

(24)

式中:γ為正常數(shù)；rk為解的殘差。

結(jié)合形函數(shù)理論與共軛梯度最小二乘正則化算法，其總結(jié)的冰壓識別技術(shù)流程圖如圖2所示，分為冰載荷識別模型構(gòu)建、結(jié)構(gòu)冰激應(yīng)變獲取和冰壓穩(wěn)定識別算法三個部分，其邏輯及遞進(jìn)關(guān)系見圖2。

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

3 冰載荷識別模型驗(yàn)證性分析

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

3.1 數(shù)值分析

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

圖3 典型冰區(qū)船舶外板結(jié)構(gòu)有限元模型Fig.3 FE model of the typical shell structure of ice-going vessel

實(shí)船冰載荷測量中，其冰激響應(yīng)信號始終存在噪聲的干擾，噪聲信號會對識別數(shù)值產(chǎn)生較大的干擾。為分析噪聲信號對冰載荷識別模型干擾程度，本文采用加性噪聲的形式進(jìn)行添加，其實(shí)際結(jié)構(gòu)測量的信號可寫作[25-26]

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

(25)

式中:εerr為含有噪聲信號的應(yīng)變信號；lnosie為一個百分?jǐn)?shù)，表示噪聲水平；std(ε)為ε標(biāo)準(zhǔn)差；rand(-1,1)為(-1,1)之間的隨機(jī)數(shù)。

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

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

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

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

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

3.2 試驗(yàn)驗(yàn)證

現(xiàn)場測量中很難布置大型加載裝置以進(jìn)行加載測試，故基于含有肋骨構(gòu)件的外板結(jié)構(gòu)載荷識別試驗(yàn)是檢驗(yàn)?zāi)Ｐ蜏?zhǔn)確性較為直接有效的方式。本文采用的外板結(jié)構(gòu)由外板、肋骨、強(qiáng)肋骨和縱桁等組成，如圖5所示。模型采用3 mm厚的6061鋁合金板，試驗(yàn)?zāi)Ｐ偷倪吔绮捎寐菟ê瓦厳l組合方式進(jìn)行約束。模型的彈性模量、材料密度和泊松比分別設(shè)為70 GPa,2 700 kg/m3和0.3。圖5右側(cè)為對應(yīng)的測量位置，測量位置的應(yīng)變片粘貼方式如圖5中的右上角方框內(nèi)的示意圖所示。應(yīng)變片的線路連接采用帶溫度補(bǔ)償片的半橋連接方式。動載荷識別試驗(yàn)中應(yīng)變信號和力信號同時由DH5922數(shù)據(jù)采集系統(tǒng)記錄，采樣頻率10 kHz。為保證肋骨與真實(shí)結(jié)構(gòu)之間的結(jié)構(gòu)相似性最大化，肋骨縱向端部與縱桁相固定，且監(jiān)測及測量的位置均遠(yuǎn)離邊界。

圖5 載荷識別試驗(yàn)裝置Fig.5 Experimental facility in the load identification

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

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

綜合數(shù)值和試驗(yàn)結(jié)果可知，基于形函數(shù)的冰載荷識別算法具有求解速度較快、穩(wěn)定性良好及精度較高的優(yōu)點(diǎn)，適用于船體外板結(jié)構(gòu)動冰載荷識別分析。

4 結(jié) 論

實(shí)船冰載荷監(jiān)測工作在國內(nèi)尚屬起步狀態(tài)，其監(jiān)測識別算法作為核心技術(shù)可直接提升監(jiān)測的整體水準(zhǔn)，進(jìn)而提高我國冰區(qū)艦船的設(shè)計(jì)建造能力?；谝陨涎芯績?nèi)容可得出以下結(jié)論及研究方向：

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

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

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