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

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

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

模擬無(wú)限流域常用的方法是通過(guò)無(wú)限域截?cái)噙吔?，將無(wú)限域截?cái)喑捎邢抻颍诮財(cái)噙吔缟鲜┘幽軠?zhǔn)確反映無(wú)限水域動(dòng)力特性的邊界條件或單元，無(wú)反射邊界條件和無(wú)限元法就是常用的方法之一，在有限元商用軟件如ABAQUS、ANSYS、LS-DYNA 等都能找到其應(yīng)用方法。采用有限元商用軟件中的無(wú)反射邊界條件或無(wú)限元法時(shí)，為得到合理準(zhǔn)確的計(jì)算結(jié)果，有限域的離散場(chǎng)要求較大，無(wú)限流域的截?cái)噙吔绨霃酵ǔＲ笫欠治鼋Y(jié)構(gòu)半徑的6 倍或以上[1]。這就意味著當(dāng)結(jié)構(gòu)模型較大時(shí)，有限流域的離散場(chǎng)將會(huì)占據(jù)較多的計(jì)算資源，使分析規(guī)模增大。

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

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

對(duì)于圖1 所示的含有無(wú)限水域的水下結(jié)構(gòu)聲固耦合系統(tǒng)，其中無(wú)限水域被截?cái)噙吔鏢 分成有限水域和無(wú)限水域，并且聲學(xué)水域滿足聲學(xué)波動(dòng)方程：

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

對(duì)于無(wú)限水域截?cái)噙吔缫詢?nèi)的水下結(jié)構(gòu)與有限水域的聲固耦合系統(tǒng)，其有限元耦合方程可寫(xiě)為：

2 無(wú)限水域PWA-FEM 總場(chǎng)公式

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

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

對(duì)于圖1 所示的聲固耦合系統(tǒng)，流場(chǎng)經(jīng)由截?cái)噙吔鏢分割成有限域和無(wú)限域，無(wú)限域的動(dòng)力特性通過(guò)截?cái)噙吔鏢上的PWA 方程描述。

在無(wú)限水域截?cái)噙吔鏢上，無(wú)限域散射場(chǎng)壓力Psc和等效加速度a′sc之間關(guān)系由式(8)求導(dǎo)后得出：

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

3 無(wú)限水域SBFEM-FEM 總場(chǎng)公式

采用SBFEM 模擬圖1 所示無(wú)限水域，根據(jù)SBFEM 理論，只需無(wú)限水域截?cái)噙吔鏢，離散模型見(jiàn)圖2，即截?cái)噙吔缟系娜我鈫卧猧 和j 具有相同的相似中心O，整個(gè)無(wú)限水域的動(dòng)力特性通過(guò)截?cái)噙吔缟系腟BFEM 控制方程描述。

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

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

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

4 算例

分別使用PWA-FEM 總場(chǎng)公式和SBFEM-FEM總場(chǎng)公式，對(duì)浸沒(méi)在無(wú)限水域中、受平面波沖擊作用的無(wú)限長(zhǎng)圓筒[17]進(jìn)行線性瞬態(tài)響應(yīng)分析，來(lái)驗(yàn)證總場(chǎng)公式的可行性和合理性，其示意圖見(jiàn)圖3，而其有限元模型見(jiàn)圖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)均為無(wú)量綱量，無(wú)量綱速度為vn/c，無(wú)量綱時(shí)間為ct/r，t為時(shí)間。

圖5 R/r=4 時(shí)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 總場(chǎng)公式對(duì)圓筒算例進(jìn)行計(jì)算分析，無(wú)限域圓形截?cái)噙吔绨霃饺A筒結(jié)構(gòu)半徑的4 倍，見(jiàn)圖4，流場(chǎng)環(huán)向單元分別取40 個(gè)、80 個(gè)、120 個(gè)，流場(chǎng)徑向單元取100 個(gè)，圓筒環(huán)向單元與流體域相同，徑向單元取2。圓筒用四節(jié)點(diǎn)平面應(yīng)變單元離散，有限水域用四節(jié)點(diǎn)平面聲學(xué)單元離散，無(wú)限水域用截?cái)噙吔缟系腜WA 總場(chǎng)公式模擬。圖5給出的是基于PWA-FEM 總場(chǎng)公式的數(shù)值解、Huang[19]的級(jí)數(shù)解析解和采用無(wú)反射邊界條件計(jì)算的數(shù)值解。數(shù)值解采用的是三維模型，水下圓筒用C3D8R 單元離散，有限流域用A3D8R 單元離散，徑向模型網(wǎng)格參數(shù)與PWAFEM 總場(chǎng)公式設(shè)置相同，軸向單元數(shù)設(shè)置為1，軸向厚度設(shè)為0.01，無(wú)限流域通過(guò)在無(wú)限流域截?cái)噙吔缭O(shè)置Planar 型無(wú)反射邊界條件模擬，流固耦合面用Tie 連接，分別采用total wave 公式和scattered wave公式計(jì)算。

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

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

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

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

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

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

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

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

圖9 給出了SBFEM-FEM 總場(chǎng)公式對(duì)該算例進(jìn)行模擬計(jì)算的結(jié)果。與PWA-FEM 總場(chǎng)公式、常規(guī)有限元方法相比，用SBFEM 模擬無(wú)限水域時(shí)，能大大減小無(wú)限水域截?cái)噙吔绲陌霃?，減少有限域單元數(shù)量和計(jì)算規(guī)模。

圖9 R/r=1.3 時(shí)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 時(shí)，SBFEM-FEM 總場(chǎng)公式的計(jì)算結(jié)果與解析解相比，出現(xiàn)了部分偏差，而在環(huán)向單元數(shù)K逐漸增加到80 和120 后，SBFEM-FEM 總場(chǎng)公式的計(jì)算結(jié)果收斂于解析解和數(shù)值解，說(shuō)明了SBFEM-FEM總場(chǎng)公式模擬無(wú)限水域的合理性和可靠性，且收斂速度要快于PWA-FEM 總場(chǎng)公式；此外，在環(huán)向單元數(shù)K=120 的情況下，SBFEM-FEM 總場(chǎng)公式計(jì)算結(jié)果與PWA-FEM 總場(chǎng)公式計(jì)算結(jié)果基本一致時(shí)，所需要的無(wú)限流域截?cái)噙吔绨霃矫黠@小于PWA-FEM 總場(chǎng)公式，說(shuō)明了SBFEM-FEM 總場(chǎng)公式可有效減小無(wú)限流域截?cái)噙吔绨霃剑瑴p小有限域離散范圍和單元網(wǎng)格數(shù)量，減小了分析規(guī)模。

現(xiàn)將數(shù)值算例的圓形截?cái)噙吔绺臑榫匦谓財(cái)噙吔鐣r(shí)，系統(tǒng)模型如圖10 所示，以便考察PWAFEM 和SBFEM-FEM 總場(chǎng)公式對(duì)其他邊界形狀的適應(yīng)性。

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

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

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

5 結(jié)論

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

(1) PWA-FEM 和SBFEM-FEM 總場(chǎng)公式均可用來(lái)模擬沖擊波作用下的含無(wú)限水域的結(jié)構(gòu)沖擊響應(yīng)問(wèn)題，該公式將沖擊波載荷等效成截?cái)噙吔缟系牡刃Ч?jié)點(diǎn)力，簡(jiǎn)化了沖擊波傳播的模擬過(guò)程；

(2) SBFEM-FEM 總場(chǎng)公式因其繼承了SBFEM半解析的特點(diǎn)，即在徑向方向(指向無(wú)窮遠(yuǎn)的方向)的解精確滿足能量吸收性，故模擬無(wú)限流域問(wèn)題時(shí)，可取較小的截?cái)喟霃剑?/p>

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

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

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