基于總場公式的水下沖擊結(jié)構(gòu)響應(yīng)分析方法

張 浩，趙梓斌，李上明,2

(1. 中國工程物理研究院總體工程研究所，四川，綿陽 621999；2. 工程材料與結(jié)構(gòu)沖擊振動四川省重點(diǎn)實(shí)驗(yàn)室，四川，綿陽 621999)

對于水下沖擊結(jié)構(gòu)響應(yīng)問題而言，準(zhǔn)確地模擬水下沖擊波作用在結(jié)構(gòu)上的載荷是影響計(jì)算結(jié)果合理性的關(guān)鍵之一。在這類流固耦合問題中，水下結(jié)構(gòu)滿足連續(xù)介質(zhì)條件，使用有限元方法離散，而水下結(jié)構(gòu)周圍的流體域則滿足無限聲學(xué)水域條件，當(dāng)作無粘無旋的理想流體即聲學(xué)波動介質(zhì)處理。通過將有限元方程與模擬無限流域的方程耦合，用于分析無限流域下的結(jié)構(gòu)沖擊響應(yīng)問題。

模擬無限流域常用的方法是通過無限域截斷邊界，將無限域截斷成有限域，在截斷邊界上施加能準(zhǔn)確反映無限水域動力特性的邊界條件或單元，無反射邊界條件和無限元法就是常用的方法之一，在有限元商用軟件如ABAQUS、ANSYS、LS-DYNA 等都能找到其應(yīng)用方法。采用有限元商用軟件中的無反射邊界條件或無限元法時，為得到合理準(zhǔn)確的計(jì)算結(jié)果，有限域的離散場要求較大，無限流域的截斷邊界半徑通常要求是分析結(jié)構(gòu)半徑的6 倍或以上[1]。這就意味著當(dāng)結(jié)構(gòu)模型較大時，有限流域的離散場將會占據(jù)較多的計(jì)算資源，使分析規(guī)模增大。

本文提出的水下沖擊結(jié)構(gòu)響應(yīng)的總場公式是基于分離場技術(shù)[2]，在不考慮空化作用條件下，將總場沖擊波壓力分解為散射場壓力和自由場壓力，自由場壓力由經(jīng)驗(yàn)公式計(jì)算，避免了計(jì)算沖擊波在流體中傳播的過程，降低了分析規(guī)模。將可用于反映無限水域動力特性的PWA 總場公式和SBFEM 總場公式，分別與流固耦合系統(tǒng)的FEM方程耦合，建立PWA-FEM 總場公式和SBFEMFEM 總場公式。PWA 方程是一種高頻端漸近式[3]，可以用于求解早期高頻響應(yīng)，而中后期響應(yīng)與DAA方法[4]相比精確度有所欠缺。但由于其表達(dá)式簡單，易于和FEM 方程耦合，所以有一定的工程應(yīng)用價值[5]。而SBFEM[6]作為一種較新的半解析解方法，在徑向方向即離散面法向方向的解是精確解，嚴(yán)格滿足無窮遠(yuǎn)處邊界條件，已成功應(yīng)用于無限水庫[7]和無限巖土[8]模擬，也可與有限元商業(yè)軟件聯(lián)合使用[9]，能準(zhǔn)確的模擬無限域問題[10]。Zhao 等[11]使用基于SBFEM 的修正高階透射邊界來模擬無限域，證明了總場公式的有效性，許賀等[12]提出了一種基于FEM-SBFEM 的壩-庫水動力耦合簡化分析方法，Yaseri 等[13]使用SBFEM-FEM 模擬土壩-彈性峽谷系統(tǒng)，計(jì)算了不同形狀的無限大彈性峽谷對土壩-彈性峽谷系統(tǒng)固有周期的影響，SBFEM 還可用于土-結(jié)構(gòu)相互作用系統(tǒng)中無限地基的模擬[14-15]。SBFEM 方程對于減小無限域問題的有限域離散半徑有重要意義，降低了分析規(guī)模，并且SBFEM 方程模擬無限域問題時，對邊界形狀要求不高，可適用范圍更廣。

1 聲固耦合系統(tǒng)的FEM 方程

對于圖1 所示的含有無限水域的水下結(jié)構(gòu)聲固耦合系統(tǒng)，其中無限水域被截斷邊界S 分成有限水域和無限水域，并且聲學(xué)水域滿足聲學(xué)波動方程：

圖1 聲固耦合系統(tǒng)Fig. 1 Acoustic fluid-structure coupling system

對于無限水域截斷邊界以內(nèi)的水下結(jié)構(gòu)與有限水域的聲固耦合系統(tǒng)，其有限元耦合方程可寫為：

2 無限水域PWA-FEM 總場公式

根據(jù)平面波近似理論[3]，壓力和流體質(zhì)點(diǎn)速度的關(guān)系如下：

式中：P為流場壓力； ρ為流體密度；c為流場中聲速；V為流體質(zhì)點(diǎn)速度。

對于圖1 所示的聲固耦合系統(tǒng)，流場經(jīng)由截斷邊界S分割成有限域和無限域，無限域的動力特性通過截斷邊界S上的PWA 方程描述。

在無限水域截斷邊界S上，無限域散射場壓力Psc和等效加速度a′sc之間關(guān)系由式(8)求導(dǎo)后得出：

式(23)為含沖擊波載荷的水下結(jié)構(gòu)聲固系統(tǒng)的全耦合方程，該公式可以用Newmark 法進(jìn)行求解。

3 無限水域SBFEM-FEM 總場公式

采用SBFEM 模擬圖1 所示無限水域，根據(jù)SBFEM 理論，只需無限水域截斷邊界S，離散模型見圖2，即截斷邊界上的任意單元i 和j 具有相同的相似中心O，整個無限水域的動力特性通過截斷邊界上的SBFEM 控制方程描述。

圖2 無限水域SBFEM 離散模型Fig. 2 SBFEM discretization of infinite fluid

式(53)即為含沖擊波載荷的SBFEM-FEM 總場公式。該公式可采用Newmark 法進(jìn)行求解分析無限水域下的結(jié)構(gòu)沖擊響應(yīng)問題。

以上兩種總場公式的推導(dǎo)，均與沖擊波載荷衰減形式無關(guān)，在已知自由場沖擊波壓力和流體加速度的情況下，即可進(jìn)行求解。

4 算例

分別使用PWA-FEM 總場公式和SBFEM-FEM總場公式，對浸沒在無限水域中、受平面波沖擊作用的無限長圓筒[17]進(jìn)行線性瞬態(tài)響應(yīng)分析，來驗(yàn)證總場公式的可行性和合理性，其示意圖見圖3，而其有限元模型見圖4。

圖3 平面沖擊波與水下圓筒Fig. 3 Plane shock wave and submerged cylindrical shell

圖4 水下圓筒與離散水域有限元模型Fig. 4 FE model of submerged cylindrical shell and discretised water domain

圖5 給 出 的是 β=0 (沖 擊波 不 衰減)、 θ=0°、90°和 180°處，圓 筒平均半徑r處 法向速度vn。θ為圓筒外法線與水平軸線的夾角。圖中橫坐標(biāo)和縱坐標(biāo)均為無量綱量，無量綱速度為vn/c，無量綱時間為ct/r，t為時間。

圖5 R/r=4 時PWA-FEM 解與Huang 的解析解和數(shù)值解的比較Fig. 5 Comparison of PWA-FEM solution with Huang's analytical solution and numerical solution when R/r=4

采用PWA-FEM 總場公式對圓筒算例進(jìn)行計(jì)算分析，無限域圓形截斷邊界半徑取圓筒結(jié)構(gòu)半徑的4 倍，見圖4，流場環(huán)向單元分別取40 個、80 個、120 個，流場徑向單元取100 個，圓筒環(huán)向單元與流體域相同，徑向單元取2。圓筒用四節(jié)點(diǎn)平面應(yīng)變單元離散，有限水域用四節(jié)點(diǎn)平面聲學(xué)單元離散，無限水域用截斷邊界上的PWA 總場公式模擬。圖5給出的是基于PWA-FEM 總場公式的數(shù)值解、Huang[19]的級數(shù)解析解和采用無反射邊界條件計(jì)算的數(shù)值解。數(shù)值解采用的是三維模型，水下圓筒用C3D8R 單元離散，有限流域用A3D8R 單元離散，徑向模型網(wǎng)格參數(shù)與PWAFEM 總場公式設(shè)置相同，軸向單元數(shù)設(shè)置為1，軸向厚度設(shè)為0.01，無限流域通過在無限流域截斷邊界設(shè)置Planar 型無反射邊界條件模擬，流固耦合面用Tie 連接，分別采用total wave 公式和scattered wave公式計(jì)算。

從圖5 可以看到，在有限域取結(jié)構(gòu)半徑4 倍的條件下，PWA-FEM 總場公式的計(jì)算結(jié)果與解析解、總波解和散波解吻合較好，隨著流場環(huán)向單元數(shù)的增加，PWA-FEM 總場公式的計(jì)算結(jié)果逐漸收斂于解析解、總波解和散波解。PWA-FEM 總場公式在環(huán)向單元較少時無法準(zhǔn)確模擬無限流域，而隨著環(huán)向單元的增加，模擬的準(zhǔn)確度也隨之增加，并逼近于解析解、總波解和散波解。這說明了PWA-FEM 總場公式模擬無限流域有一定的合理性和可信度。

圖5 中的PWA-FEM,K=40,R/r=4 計(jì)算曲線出現(xiàn)了較大波動，為探究該現(xiàn)象發(fā)生的內(nèi)在機(jī)理原因，提取了圓筒算例的流固耦合面上 θ=0處的聲學(xué)流體總壓進(jìn)行分析，計(jì)算結(jié)果見圖6。

圖6 R/r=4 時，流固耦合面上 θ=0處的聲學(xué)流體總壓Fig. 6 Acoustic fluid total pressure of fluid-structure interaction at θ =0 when R/r=4

從圖6 可以看出，在無限域截斷半徑不變的條件下，環(huán)向單元數(shù)K=40 時，流固耦合面上θ=0處的聲學(xué)流體總壓是振蕩變化的，但總體呈穩(wěn)定變化趨勢。圖6 的PWA-FEM,K=40,R/r=4 計(jì)算曲線出現(xiàn)誤差的原因是PWA-FEM 總場公式的環(huán)向單元數(shù)較少，流體有限元網(wǎng)格過大，隨著網(wǎng)格加密，計(jì)算結(jié)果逐步收斂。

圖7 是在環(huán)向單元數(shù)取120 時，取不同無限流域截斷半徑條件下的計(jì)算結(jié)果。從計(jì)算結(jié)果可看出，在環(huán)向單元數(shù)不變的條件下，隨著有限流域半徑的增加，計(jì)算結(jié)果也逐漸收斂于解析解、總波解和散波解。這說明了使用PWA-FEM 總場公式模擬無限流域時，仍需要較大的有限域才能準(zhǔn)確的模擬無限域。

圖7 中可以明顯看出，PWA-FEM,R/r=2.2,K=120 的計(jì)算曲線在無量綱時間ct/r=4～5 范圍內(nèi)出現(xiàn)較大誤差，為進(jìn)一步探究該誤差出現(xiàn)的內(nèi)在機(jī)理原因，提取了圓筒算例的流固耦合面上 θ=π處的聲學(xué)流體總壓進(jìn)行分析，計(jì)算結(jié)果見圖8。

圖7 K=120 時PWA-FEM 解與Huang 的解析解和數(shù)值解的比較Fig. 7 Comparison of PWA-FEM solution with Huang's analytical solution and numerical solution when K=120

圖8 K=120 時，流固耦合面上 θ=π處的聲學(xué)流體總壓Fig. 8 acoustic fluid total pressure of fluid-structure interaction at θ=π when K=120

圖8 是在環(huán)向單元數(shù)等于120，取不同無限流域截斷半徑條件下，流固耦合面上 θ=π處的聲學(xué)流體總壓。從圖中可以看出，在無限流域截斷半徑較小時，PWA-FEM 求解的壓力在后期計(jì)算中出現(xiàn)了較大波動，而在無限流域截斷半徑較大時，PWA-FEM 求解的總壓與ABAQUS 總波公式計(jì)算的結(jié)果非常吻合。出現(xiàn)這種現(xiàn)象的原因是有限流域較小時，PWA-FEM 總場公式對邊界的吸收特性未能準(zhǔn)確的模擬，導(dǎo)致了部分散射波在邊界反射重新進(jìn)入計(jì)算域，進(jìn)而使流體總壓出現(xiàn)較大波動，所以圖7 中PWA-FEM,R/r=2.2,K=120 計(jì)算曲線在無量綱時間ct/r=4～5 范圍內(nèi)會出現(xiàn)較大誤差。

圖9 給出了SBFEM-FEM 總場公式對該算例進(jìn)行模擬計(jì)算的結(jié)果。與PWA-FEM 總場公式、常規(guī)有限元方法相比，用SBFEM 模擬無限水域時，能大大減小無限水域截斷邊界的半徑，減少有限域單元數(shù)量和計(jì)算規(guī)模。

圖9 R/r=1.3 時SBFEM-FEM 解與Huang 的解析解、數(shù)值解和PWA-FEM 解的比較Fig. 9 Comparison of SBFEM-FEM solution with Huang's analytical solution , numerical solution and PWA-FEM solution when R/r=4

從圖9 可以看出，在環(huán)向單元數(shù)K=40 時，SBFEM-FEM 總場公式的計(jì)算結(jié)果與解析解相比，出現(xiàn)了部分偏差，而在環(huán)向單元數(shù)K逐漸增加到80 和120 后，SBFEM-FEM 總場公式的計(jì)算結(jié)果收斂于解析解和數(shù)值解，說明了SBFEM-FEM總場公式模擬無限水域的合理性和可靠性，且收斂速度要快于PWA-FEM 總場公式；此外，在環(huán)向單元數(shù)K=120 的情況下，SBFEM-FEM 總場公式計(jì)算結(jié)果與PWA-FEM 總場公式計(jì)算結(jié)果基本一致時，所需要的無限流域截斷邊界半徑明顯小于PWA-FEM 總場公式，說明了SBFEM-FEM 總場公式可有效減小無限流域截斷邊界半徑，減小有限域離散范圍和單元網(wǎng)格數(shù)量，減小了分析規(guī)模。

現(xiàn)將數(shù)值算例的圓形截斷邊界改為矩形截斷邊界時，系統(tǒng)模型如圖10 所示，以便考察PWAFEM 和SBFEM-FEM 總場公式對其他邊界形狀的適應(yīng)性。

圖10 水下圓筒與矩形邊界Fig. 10 Underwater cylinder and rectangular boundary

圖11 是矩形截斷邊界條件下，PWA-FEM 總場公式和SBFEM 總場公式的計(jì)算結(jié)果。在截斷邊界為方形、環(huán)向單元數(shù)K=120 的條件下，SBFEMFEM 總場公式的計(jì)算結(jié)果一定程度上優(yōu)于PWAFEM 總場公式，與圓形邊界相比結(jié)果相差不大?？傮w看，文中給出的兩類總場公式均能適應(yīng)不同邊界形狀。

圖11 不同邊界形狀的計(jì)算結(jié)果對比Fig. 11 Comparison of calculation results of different boundary shapes

5 結(jié)論

本文基于分離場技術(shù)提出的總場公式理論具有以下特點(diǎn)：

(1) PWA-FEM 和SBFEM-FEM 總場公式均可用來模擬沖擊波作用下的含無限水域的結(jié)構(gòu)沖擊響應(yīng)問題，該公式將沖擊波載荷等效成截斷邊界上的等效節(jié)點(diǎn)力，簡化了沖擊波傳播的模擬過程；

(2) SBFEM-FEM 總場公式因其繼承了SBFEM半解析的特點(diǎn)，即在徑向方向(指向無窮遠(yuǎn)的方向)的解精確滿足能量吸收性，故模擬無限流域問題時，可取較小的截斷半徑；

(3) PWA-FEM 總場公式形式簡單，易編程實(shí)現(xiàn)，計(jì)算效率較高，但其無限水域的離散域相對較大，建議取水下結(jié)構(gòu)長度的3 倍左右；

(4) 當(dāng)前兩類方法均可適用于任意形狀邊界的無反射特性的模擬；

(5) PWA-FEM 和SBFEM-FEM 總場公式在推導(dǎo)和求解過程中與計(jì)算結(jié)構(gòu)規(guī)模大小無關(guān)，在本文中僅用較為簡單的小規(guī)模算例進(jìn)行驗(yàn)證，而對于更大規(guī)模的真實(shí)結(jié)構(gòu)則暫時未進(jìn)行驗(yàn)證。