基于廣義多項式混沌法的圓柱殼自由振動不確定性分析

陳爐云,鄧乃旗,楊 念

（1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.中國艦船研究設(shè)計中心，武漢 430064）

0 引 言

理論分析和數(shù)值計算是艦船結(jié)構(gòu)聲振特性預(yù)報的主要方法。對于艦船系統(tǒng)這樣的復(fù)雜對象，即使實施精確的數(shù)字建模也存在一定的偏差，這些偏差不是由于設(shè)計或建造的錯誤因素所導(dǎo)致的，而是由結(jié)構(gòu)固有屬性所決定。由于存在著不確定性因素，結(jié)構(gòu)響應(yīng)參數(shù)僅是真實結(jié)果的近似值，兩者間的偏差可能超出設(shè)計許可，并可能導(dǎo)致重大工程風險。在艦船結(jié)構(gòu)建造過程中，存在著影響其聲振特性的固有的、具有隨機特性的影響參數(shù)，即建造因素不確定性。由于存在建造因素不確定性，使得艦船結(jié)構(gòu)聲振特性參數(shù)實質(zhì)上是一個具有一定概率密度分布的隨機變量，這表明在一定建造工藝水平下，不能消除該不確定性對系統(tǒng)的影響，即使提高工藝水平也只能使不確定性區(qū)間變窄而不能消除。建造因素不確定性包括材料參數(shù)不確定性、幾何尺寸不確定性、邊界條件不確定性、初始條件不確定性、操作環(huán)境不確定性、載荷不確定性等[1-4]。當結(jié)構(gòu)中存在多個不確定性輸入?yún)?shù)時，結(jié)構(gòu)響應(yīng)的預(yù)測值與實際值會存在巨大的偏差[5-7]。

對于結(jié)構(gòu)聲振特性不確定性問題，需評估各參數(shù)的不確定性對最終響應(yīng)的影響，并用概率統(tǒng)計方法對這些影響進行描述以獲得綜合影響程度。針對已知概率密度分布的不確定性輸入，采用數(shù)學模型將輸入不確定性參數(shù)傳導(dǎo)給輸出響應(yīng)，獲得輸出參數(shù)的概率分布特征。不確定性問題的主要計算方法包括攝動法[8]、蒙特卡洛法[9]、可靠性分析法[10]、響應(yīng)面法[11]和譜分解方法[12]等等。譜分解方法在處理高斯分布和非高斯分布問題的計算效率方面有一定優(yōu)勢。為解決非高斯分布和參數(shù)協(xié)方差未知問題，廣義多項式混沌方法（generalized polynomial chaos,GPC）成為當前研究熱點。

本文將GPC 法應(yīng)用于圓柱殼結(jié)構(gòu)自由振動的不確定性分析，將自由振動不確定性研究范圍從隨機載荷不確定性方面擴展到考慮建造因素的不確定性問題。在文中建立自由振動不確定性分析模型，將結(jié)構(gòu)建造因素中的不確定性傳導(dǎo)并反映到圓柱殼結(jié)構(gòu)自由振動特性概率密度分布。研究成果可用于考慮建造質(zhì)量影響的艦艇結(jié)構(gòu)聲振特性的預(yù)報與控制。

1 廣義多項式混沌方法

1.1 廣義多項式混沌方法概述

在不確定性分析中，標準隨機變量的選擇將影響計算收斂性和截斷項數(shù)，選擇最佳匹配多項式可提高計算效率，如表1 所示。在聲振特性不確定性計算中，采用Hermite 多項式擬合不確定性參數(shù)較為合適，其最高階為3次。

表1 不確定參數(shù)對應(yīng)的最佳多項式類型Tab.1 Uncertain parameter types and the corresponding orthogonal polynomials

1.2 廣義多項式混沌方法的隨機投影

選定多項式類型后，定義隨機空間(Ω,γ,P)，其中Ω是總樣本空間、γ是樣本空間中的σ域、P是事件概率。對于隨機空間，空間中任意不確定性變量χ可表示為

式中，Ψi為隨機正交基，xi為待定系數(shù)。隨機正交基Ψi是一組含有隨機向量ξ的多項式，它們滿足正交性條件：

1.3 廣義多項式混沌方法的求解

應(yīng)用廣義多項式混沌方法，將結(jié)構(gòu)中建造非確定性因素作為輸入?yún)?shù)通過結(jié)構(gòu)控制方程傳遞到輸出的響應(yīng)結(jié)果中，開展非高斯分布的隨機問題求解。

應(yīng)用GPC 法求解不確定性問題的步驟為：（1）將不確定性參數(shù)表示成多項式混沌展開式形式，求解該展開式系數(shù)；（2）用廣義多項式混沌表達式表達不確定性輸入?yún)?shù)和所求解響應(yīng)結(jié)果；（3）將輸入?yún)?shù)不確定性通過結(jié)構(gòu)的隨機控制方程傳導(dǎo)至輸出參數(shù)中，求得響應(yīng)表達式中的未知系數(shù)。GPC法求解過程如圖1所示。

圖1 基于GPC的求解過程示意圖Fig.1 GPC-based solving flow chart

設(shè)結(jié)構(gòu)聲振方程滿足線性條件，將圖1所示的求解步驟應(yīng)用于結(jié)構(gòu)聲振特性不確定性分析中，需獲得結(jié)構(gòu)中輸入?yún)?shù)（如材料密度、板殼厚度、肋骨型號等）的概率密度分布。

通過不確定性傳導(dǎo)功能，將結(jié)構(gòu)參數(shù)不確定性特性傳導(dǎo)至輸出參數(shù)（即結(jié)構(gòu)聲振特性），獲得結(jié)構(gòu)聲振特性概率密度分布，求解過程如圖2所示。

圖2 基于GPC的結(jié)構(gòu)聲振特性不確定性求解Fig.2 Dynamical solution based on GPC

2 不確定性輸入?yún)?shù)

在應(yīng)用GPC進行結(jié)構(gòu)聲振特性不確定性求解中，需求解式（4）中待定系數(shù)xi。設(shè)變量χ和ξ處于同一隨機空間中且具有相同的概率分布，采用Galerkin投影法描述待定系數(shù)xi：

式中，a、b是標準隨機變量上、下限。設(shè)式中χ(ξ)為未知函數(shù)，需進行隨機空間轉(zhuǎn)換以獲得不確定性變量χ關(guān)于參數(shù)ξ的表達式χ=F(ξ)，變換方程如圖3所示。

圖3 不同隨機分布之間的變換Fig.3 Conversion between uncertain parameters

在圖3 中，不確定性變量χ在隨機空間Ω1中，概率分布為ρ1(χ)，區(qū)間為[c,d]；隨機參數(shù)ξ在隨機空間Ω2中，概率分布為ρ2(ξ)，區(qū)間為[a,b]。不確定性變量χ和標準隨機變量ξ間的變換滿足：

對于大多數(shù)隨機分布函數(shù)，式(7)中隨機分布函數(shù)F(ξ)并不是一個顯性表達式，需引入標準平均分布為中間變量實現(xiàn)對待定系數(shù)xi的轉(zhuǎn)換方程，將式（7）代入式（6）中，可得

建立不確定性變量χ和標準隨機變量ξ對于標準平均分布μ(0,1)的空間變換方程：

3 不確定性傳導(dǎo)方程

3.1 不確定性傳導(dǎo)

在結(jié)構(gòu)聲振特性不確定性分析中，研究輸入?yún)?shù)不確定性對聲振特性的影響，將輸入?yún)?shù)（結(jié)構(gòu)材料、結(jié)構(gòu)幾何尺寸等）的不確定性特性通過運動控制方程有效地傳導(dǎo)到輸出參數(shù)（結(jié)構(gòu)響應(yīng)、聲學響應(yīng)等）概率密度分布中。

在應(yīng)用有限元法進行確定性固有頻率系數(shù)λri分析中，涉及確定性剛度矩陣系數(shù)ki、確定性質(zhì)量矩陣系數(shù)mi和確定性振型向量系數(shù)φri。由于涉及大量積分運算且非常耗時，需對傳導(dǎo)方程進行改進。

3.2 配點法概述

對于大多數(shù)結(jié)構(gòu)，其變量χ的概率密度分布ρ1(χ)無法通過明確的函數(shù)來表示，采用基于配點法（collection method，CM）結(jié)合試驗數(shù)據(jù)獲得ρ1(χ)的概率密度分布具有一定優(yōu)勢。配點法基本原理是基于最小二乘法原理，使兩種預(yù)測結(jié)果間的誤差達到最小，以確定相關(guān)的不確定性參數(shù)。

求解式（13），獲得廣義多項式混沌系數(shù)xi，進而獲得χ的概率密度。配點法可直接避開通過式（8）中復(fù)雜積分運算求解系數(shù)的過程，減小計算復(fù)雜度和計算量。

3.3 配點法求解響應(yīng)混沌系數(shù)

為提高計算效率，將配點法和廣義多項式混沌法進行結(jié)合，用于作為輸入?yún)?shù)的不確定性變量混沌系數(shù)及作為輸出的結(jié)構(gòu)響應(yīng)表達式中混沌系數(shù)的求解，獲得結(jié)構(gòu)響應(yīng)概率密度分布，其求解過程如圖4所示。

圖4 配點法求解響應(yīng)混沌系數(shù)Fig.4 Flow chart for obtaining response coefficients based on collocation method

從圖4 可知，應(yīng)用配點法結(jié)合GPC 法求解非確定性問題時，可避免應(yīng)用確定性剛度矩陣系數(shù)、確定性質(zhì)量矩陣系數(shù)、確定性振型向量系數(shù)求解確定性固有頻率系數(shù)中的推導(dǎo)和計算?？芍苯永贸墒煊邢拊虡I(yè)軟件進行求解，且可依托商業(yè)軟件的計算優(yōu)勢簡化計算過程，使計算效率最大化。

4 數(shù)值分析

4.1 模型概述

以圓柱殼模型為例，采用Galerkin 投影方法確定不確定性參數(shù)，分別采用GPC 法和蒙特卡洛法計算固有頻率，驗證GPC法的正確性。

具有兩端簡支圓柱殼結(jié)構(gòu)的幾何參數(shù)如表2所示。結(jié)構(gòu)材料為鋼，力學特性為：密度ρ=7860 kg/m3，彈性模量Em=2.1×1011Pa，泊松比μ=0.3。

表2 圓柱殼幾何參數(shù)Tab.2 Cylindrical shell geometry parameters

4.2 彈性模量概率密度

彈性模量概率密度分布圖如圖5所示。

圖5 彈性模量概率密度分布Fig.5 PDF of elastic moduli

4.3 固有頻率概率計算

基于Hermite 多項式，結(jié)合GPC 法和配點法計算圓柱殼結(jié)構(gòu)固有頻率不確定性，計算步驟為：（1）定義圓柱殼結(jié)構(gòu)幾何參數(shù)、材料參數(shù)及彈性模量概率密度分布函數(shù)；（2）求解圓柱殼結(jié)構(gòu)的確定性固有頻率解；（3）將不確定性參數(shù)表示成多項式混沌展開式的形式；（4）推導(dǎo)含建造因素的圓柱殼結(jié)構(gòu)自由振動方程；（5）引入隨機變量，將建造不確定性參數(shù)定義為不確定性變量；（6）基于Hilbert 空間，引入正交基Ψi；（7）應(yīng)用逆冪法（inverse power method）求解自由振動方程的特征值；（8）設(shè)定方程收斂標準，進行收斂性判斷；（9）應(yīng)用Galerkin投影法求解多項式混沌表達式中系數(shù)xi；（10）進行結(jié)構(gòu)固有頻率非確定性數(shù)值計算；（11）將計算結(jié)果與蒙特卡洛法進行對比。

圖6 所示為用GPC 和不同模擬次數(shù)的蒙特卡洛法計算第一階固有頻率值對比，結(jié)構(gòu)固有頻率的變化滿足正態(tài)分布。

圖6 GPC與蒙特卡洛法計算對比Fig.6 Comparison between GPC and different order MC

圖7 分別為用GPC 法、5000 次模擬蒙特卡洛法計算的圓柱殼結(jié)構(gòu)前4 階固有頻率的概率密度分布，圖中豎直虛線表示為圖5 中所定義的彈性模量期望值對應(yīng)的固有頻率。

由圖7 可知，不同階結(jié)構(gòu)固有頻率的變化滿足正態(tài)分布。采用GPC 法計算結(jié)果同蒙特卡洛法計算結(jié)果具有較好的一致性，GPC 法在求解動力學問題方面的有效性得到了驗證。

圖7 GPC法同5000次模擬蒙特卡洛方法的固有頻率對比Fig.7 Natural frequency comparison between GPC method and 5000 MC

4.4 計算效率驗證

為驗證GPC 法計算效率，對比GPC 法和不同模擬次數(shù)蒙特卡洛法的計算時間，如表3所示。

由表3 可知，在計算精度基本相似情況下，采用5000 次模擬蒙特卡洛法需137 247.6 s，而采用GPC 法僅需要184.2 s，計算時間遠小于蒙特卡洛法。圖7 和表3 表明，GPC 法在求解實際大型結(jié)構(gòu)工程不確定性振動問題時具有明顯的優(yōu)勢和潛力。

表3 GPC法同蒙特卡洛法的計算時間對比Tab.3 Consuming time of GPC and Monte Carlo

5 結(jié) 語

本文針對艦艇建造過程中的聲振特性不確定性問題，推導(dǎo)了GPC 法對不確定性參數(shù)的求解公式及不確定性傳導(dǎo)公式。以圓柱殼結(jié)構(gòu)為例，將GPC 法與有限元數(shù)值方法相結(jié)合開展了結(jié)構(gòu)動力學的非確定性分析。通過數(shù)值計算，驗證了所提出方法的可行性和有效性。研究成果可應(yīng)用于艦艇結(jié)構(gòu)聲振特性的精確預(yù)報，并可為艦艇聲學質(zhì)量控制提供策略。