回轉體高速傾斜入水的流場特性及結構響應*

高英杰，孫鐵志,2，張桂勇,2,3，尤天慶，殷志宏，宗 智,2,3

（1. 大連理工大學船舶工程學院遼寧省深海浮動結構工程實驗室，遼寧大連116024；2.大連理工大學工業裝備結構分析國家重點實驗室，遼寧大連116024；3.高新船舶與深海開發裝備協同創新中心，上海200240；4.北京宇航系統工程研究所，北京100076）

對結構高速入水的研究，在導彈、魚雷入水、航行器水上回收、船舶砰擊和水上飛機迫降等實際問題上具有重要的工程意義。魚雷等武器多采用薄殼結構形式，在高速入水過程將面臨著結構破壞變形、內部儀器失靈等問題，魚雷是否能進入正常的作戰狀態決定于魚雷的毀傷性能。回轉體入水過程涉及氣-液-固三相耦合作用[1]，在高速入水過程中流固耦合現象更復雜。回轉體在高速入水瞬間受到極大的沖擊載荷，結構可能發生較大變形甚至破壞。因此，在研究過程中需考慮流固耦合作用，分析結構的動響應特征，為結構強度設計提供依據。在對高速入水過程的流固耦合研究中，理論和實驗的方法往往難以實施。隨著計算機和數值技術的發展，數值方法在解決入水問題上不受幾何外形和邊界條件等因素的限制，具有一定的優越性。

對入水問題，已有了大量的理論和實驗研究。1929年，von Karman[2]提出并采用了附加質量和動量守恒的方法研究入水沖擊問題，在理論上求解了二維楔形體入水的沖擊載荷。在此基礎上，1932年，Wagner[3]采用伯努利方程，得到了小攻角楔形體浸水面上的壓力分布。他們為后來入水問題的理論研究奠定了基礎。Gavrilenko等[4]考慮流體的可壓縮性，對邊界條件進行Laplace變化，在理論上分析了不同形狀的軸對稱剛體結構的入水。Worthington 等[5]最早采用高速攝影技術，完整地觀察了小球入水的完整過程，對過程中的自由液面變化和空泡演變過程進行了分析。May 等[6-7]研究了不同結構參數的鋼球入水空泡特性，并分析了不同速度入水下鋼球的載荷參數。Abelson[8]開展了錐頭圓柱體以速度76.2 m/s入水的實驗，分析了入水角度與壓力分布的關系。隨著數值計算能力的發展，數值模擬方法逐步成為解決入水問題的重要研究方法。Panahi[9]基于有限體積方法和動網格技術，分析了結構物入水過程中的運動規律和動力特征。Erfanian 等[10]考慮流固耦合現象，對球鼻艏入水過程進行了數值模擬，提高了數值精度。Neaves等[11]引入空化模型和Tait 狀態并考慮流體的可壓縮性，開展了高速射彈垂直入水過程的數值模擬，較好地分析了流域的壓力場分布情況。

目前，在我國，對于入水問題也有了廣泛的研究。孫輝等[12]通過實驗研究了楔形體入水流固耦合動響應參數特征和不同參數的結構響應規律。張偉等[13]、郭子濤等[14]進行了速度35～160 m/s、不同頭型參數的入水實驗，研究了頭型參數對彈道穩定性的影響并得到了圓柱形結構入水速度衰減的理論公式。陳誠等[15]開展了超空泡航行器的小角度傾斜入水實驗，研究了入水速度和空化器面積對沖擊載荷的影響。Yan 等[16]開展了AUV 航行器高速入水實驗，并基于有限元和SPH 粒子法耦合的方法對航行器的沖擊載荷進行了分析。錢鋮鋮等[17]采用動網格技術對射彈高速入水進行了數值模擬，得到了阻力、壓力和入水速度的關系。張佳悅等[18]采用VOF方法和動網格技術對航行體高速垂直入水過程進行了數值模擬，得到了航行器阻力系數與入水速度的關系及空泡發展規律。何春濤等[19]用有限體積方法結合動網格技術，對結構垂直入水進行了數值模擬，通過與實驗結果對比驗證了數值算法的有效性，進一步分析了空泡閉合特性及空泡演變規律。侯昭等[20]基于雷諾平均方法并采用重疊網格技術，建立了六自由度圓柱傾斜入水的模擬方法，深入討論了空泡形態、彈道特性和受力情況。黃志剛等[21]采用ALE方法，深入分析了不同壁厚回轉體100 m/s入水時沖擊特性和結構強度，指出結構形式是影響入水回轉體結構強度的重要因素。

由以上研究可以看出，入水沖擊問題受到廣泛關注并有了大量有價值的研究成果。目前，研究重點多集中在低速入水問題，作為一個沖擊時間極短的強非線性問題，高速入水過程中涉及流固耦合作用的實驗和數值研究開展較少。因此，采用流固耦合數值方法，開展高速入水過程中彈塑性材料結構的動響應規律、流體動力特性和空泡演變規律的研究，可以加深對高速入水流場特征的理解，為高速入水結構物設計提供依據，同時為下一步開展高速入水數值計算提供參考。

1 數值計算模型

1.1 流體基本控制方程

在流體求解器中，采用均質平衡流模型中的VOF模型作為多相流模型。VOF方法將氣、汽、液組成的混合介質看作一種可變密度的單一流體，在整個流域內實現速度場，壓力場共享，各相介質成分由各相體積分數α 來描述。該方法可以有效地捕捉自由液面處的變形，在模擬多相流界面之間運動方面應用廣泛。

連續性方程為：

湍流模型采用了應用范圍較廣的SSTk-w模型，該模型綜合了k-w和k-ε 模型的優勢，在近壁面處采用k-w模型而在邊界層外采用k-ε 模型。SSTk-w湍流模型包含了修正的湍流黏性公式，考慮了湍流的剪切應力效應，能更精確地模擬反壓力梯度引起的分離點和分離區大小。其輸運方程為：

式中：k為湍動能，w為湍流耗散率，Gk為層流速度梯度產生的湍動能，Gw為湍動能速度梯度，Yk為波動擴張引起的湍動能耗散，Гk、Гw分別為k、w的有效擴散項，Dw為正交擴散項，Sk、Sw分別為用戶自定義源項。Гk、Гw分別為：

式中：σk、σw分別為k和w的普朗特數，F1和F2為混合函數，μt為湍流黏性系數，a1=0.31為模型參數，S為應變率，對于高雷諾數模型a*=1。

1.2 固體求解器單元和材料模型

魚雷殼體等水下武器一般采用薄殼結構，為了使回轉體在高速入水過程中的結構動響應特性更明顯，本文中在傾斜入水數值模擬算例的固體求解器中采用殼單元模型。固體求解器中采用數值積分方法，分別計算厚度方向每個積分點的應力和應變，并允許非線性材料行為。采取通用殼單元S4R 形式，該單元模型在多數情況下能提供比較精確的模擬結果。

材料本構關系等材料性能數據對數值模擬的精確性和有效性的影響很大，尤其對高速入水沖擊等非線性問題。Johnson-Cook 本構關系是一個經驗性的黏塑性模型，用于高應變率和高溫情況[22]。Johnson-Cook 材料模型形式簡單、使用方便，能夠很好地反映金屬材料的性能，在沖擊動力學上有廣泛的應用。其應力應變關系為：

45鋼的主要材料參數[23]分別為：密度7 850 kg/m3，楊氏模量210 GPa，泊松比0.29；Johnson-Cook 本構關系參數分別為：A=595 MPa，B=580 MPa，C=0.023，n=0.133，m=2.03，Tm=1 492℃，Tt=25℃。

1.3 耦合求解模式

STAR-CCM+和ABAQUS的雙向耦合是一種強流固耦合形式，雙向耦合考慮流體的壓力和剪切力對固體作用引起變形，也同時考慮變形導致流體流動的變化，當固體求解器和流體求解器都采用隱式格式時，該耦合模式允許STAR-CCM+與ABAQUS在一個時間步內多次交換數據。因此，需增加一個內部迭代過程。在內迭代耦合計算時，固體求解器先根據流場中的初始條件計算結構的變形量和位移量估值，耦合交界面處STAR-CCM+流體域界面上的載荷和ABAQUS求解器中耦合界面節點處的載荷需要滿足載荷邊界條件：

式中：ds、df為耦合交界面處固體和流體位移。

為落實監管責任，強化打擊力度，犍為縣局在30個鄉鎮建立工作站，構建起縣、鄉、村、組四級“橫向到邊、縱向到底”的網格化監管體系；通過打造A級餐飲示范店，樹立示范標桿，引導全縣行業看齊；積極推進“兩法銜接”部門聯動工作，建立完善聯席會議長效機制；以專項行動帶動日常監管，嚴格執行“四個最嚴”，三年來共嚴肅查處食品藥品領域各類違法行為680余件，扣押、銷毀不合格產品100余噸，罰沒入庫達600余萬元，移送司法機關9件11人，已判刑7人。

流體求解器采用網格變形技術實現網格的變形和移動，完成一步耦合計算并繼續更新壓力和剪切力，重復至最大時間或滿足位移收斂條件后進行下一步耦合。在雙向耦合過程中，信息在流體求解器和固體求解器中進行更新和傳遞，并在一定間隔內進行多次交換，稱為耦合時間步。耦合時間步的選擇取決于問題的耦合強度，對于復雜的流固耦合問題，耦合時間步一般與固體及流體計算的時間步相同，以實現收斂和穩定性。過程如圖1所示。

在流固耦合過程中，固體結構被施加形變和剛體運動，因此采用網格變形技術且需適應剛體運動影響。STAR-CCM+采用morphing 網格變形技術，在固體網格形變結束后，流體網格通過映射到固體網格節點上的方式產生變形移動。在有限元分析中，此方法稱為任意拉格朗日-歐拉(ALE)方法。在計算流體力學中，當流體網格與固體網格共形時，流體求解器通過額外求解輸運方程獲得固體的運動信息，通過保角映射來獲取固體的形變信息。

圖1 雙向耦合求解過程Fig.1 Computing processof fluid-structure interaction method

1.4 計算模型設置

為了進一步研究傾斜入水過程中結構動響應情況和流場變化，采用雙向流固耦合方法，開展回轉體傾斜60°入水過程的數值模擬分析。如圖2所示，計算區域分為空氣域和水域，空氣域尺寸為1.5 m×1.0 m×0.5 m，水域尺寸為1.5 m×1.0 m×1.2 m。回轉體入水速度為60 m/s，與水平面夾角為60°，彈體長0.197 m、直徑0.05 m、壁厚3 mm，采用45鋼Johnson-Cook 材料模型。

設計算域四周和底部為速度入口，頂部為壓力出口，避免壁面反射的影響模擬無限水域；設圓柱體表面為壁面，設重疊網格與流體域的界面為重疊網格邊界。如圖3所示，計算時在水面運動區域和圓柱體周圍進行網格加密，其中流體網格數量為3 285 858，固體網格數量為16 432。

圖2 計算域Fig.2 Computational domain

圖3 流體網格模型Fig.3 Mesh of fluid domain

2 計算結果和分析

2.1 方法有效性驗證

為了驗證本次數值模擬方法的有效性，選取文獻[24]中的實驗模型進行計算，并與實驗結果進行對比。實驗系統主要包括發射系統、水箱、高速攝影機、照明裝置等，水箱尺寸為3.3 m×3.3 m×6.0 m。實驗模型為45鋼實心圓柱體，直徑6 mm，實測質量4.88 g，入水速度106.6 m/s。數值模擬中，模型尺寸、質量和入水速度與實驗一致，45鋼的彈塑性階段使用Johnson-Cook 材料模型。通過對比數值模擬和入水實驗結果中的速度衰減、位移和不同時刻的空泡演變規律，驗證流固耦合算法的有效性。

圖4為速度曲線的流固耦合模擬和實驗結果。在回轉體高速撞水瞬間，頭部會受到一個極大的沖擊力，體現為在觸水瞬間回轉體存在一個劇烈的速度衰減，雙向流固耦合模擬可以清楚地捕捉到該現象。同時也可以看出，速度曲線的模擬結果和實驗結果吻合較好。由圖5可以看出，侵蝕位移的數值模擬結果和實驗結果趨近一致，證明了流固耦合數值模擬方法的有效性。

圖4 速度衰減曲線Fig.4 Velocity attenuation curves

圖5 侵蝕位移曲線Fig.5 Erosion displacement curves

圖6給出了回轉體以速度106.8 m/s入水時，不同時刻實驗和數值模擬的空泡形態對比。由圖6可以清楚地觀察到撞水開空泡、空泡發展、空泡閉合、空泡潰滅等過程。在觸水階段，回轉體撞擊水面對水產生沖擊波，空泡開始形成。在0.8 ms時，回轉體完全入水，空泡完全包裹回轉體，水表面收到沖擊產生飛濺現象。在2.0 ms時，隨著回轉體繼續運動，空泡逐漸拉長、體積不斷擴展、直徑將要達到最大，并開始產生回射流。可以看出，數值模擬的空泡演變結果與實驗結果具有較高的一致性，進一步驗證了流固耦合算法的可行有效。

2.2 速度60 m/s傾斜入水

圖6 不同時刻的空泡形態Fig.6 Cavity features at different times

圖9為不同時刻回轉體底部的應力分布，可以清楚地觀察到應力主要集中在回轉體底部中心和圓周邊緣。在觸水瞬間，圓周邊緣部分與水接觸，受到極大的沖擊載荷，應力集中在觸水圓周邊緣處；隨后，應力向底部中心處轉移，擴展至單元中心處和整個圓周邊緣；最終，應力始終集中在底部中心處和圓周處，保持穩定狀態。圖10～11為回轉體底部中心單元、半徑中心單元和圓周處單元的有效應力、應變曲線，可以看出：底部中心應力的峰值為171 MPa，大于半徑中心處和圓周處的；底部中心應力衰減頻率大于半徑中心和圓周處的，在18.8 ms左右應力發生波動產生脈沖峰值，三者應變情況與應力情況趨勢保持一致。

圖7 回轉體軸向、底部徑向的單元分布Fig.7 Element distribution of revolution body

圖8 回轉體的應力沿軸向分布曲線Fig.8 Equivalent stress curves of elementsalong axis

圖12～13為不同時刻空泡輪廓圖和流場壓力云圖，可進一步分析回轉體入水過程中結構響應的脈動情況。在約17.0 ms時，空泡發生表面閉合，在約2.0 ms時，尾部空泡收縮使空泡內部產生一個短暫的高壓區；由于回轉體阻力主要為壓差阻力，空泡內部的高壓區使回轉體首尾壓差降低，回轉體加速度呈下降趨勢；空泡內部壓力升高，導致回轉體整體受壓變大，應力、應變相應增加。在19.5 ms時，空泡尾部收縮完成，空泡內部高壓區不明顯，回轉體整體加速度變大且應力應變減小。在回轉體的有效應力曲線上，這個現象體現為17.0～19.0 ms出現的脈動峰值。

圖11 回轉體底部的應變沿徑向分布曲線Fig. 11 Strain curvesof bottom elementsalong radius

圖12 不同時刻的空泡演變Fig.12 Cavity featuresat different times

圖13 不同時刻的應力分布Fig.13 Stress distributionsat different times

采用流體固體求解器協同模擬方法，可以清楚地觀察到回轉體傾斜入水過程的完整空泡演變過程，如圖14所示。回轉體經歷撞水開空泡、空泡形成、空泡發展、空泡頸縮和空泡閉合的過程。初始時刻，回轉體撞水后，水面發生流動分離向回轉體兩側運動，撞水部分右側的水面隆起，回轉體軸向空泡呈現不對稱形態，軸向右側空泡大于左側，隨著入水時間增加，空泡不對稱性變弱；在4 ms時，空泡完全包裹回轉體；在12 ms前，空泡沿徑向和回轉體軸向擴張，空泡充分發展；在17 ms左右，空表面閉合，產生回射流；然后，空泡尾部開始收縮，空泡即將開始潰滅。

圖14 速度60 m/s 傾斜入水的空泡演變Fig. 14 Cavity features for oblique water-entry at v0=60 m/s

2.3 入水速度的影響

為了研究不同入水速度對于回轉體的載荷特性和流場演變的影響，基于已建立的流固耦合方法，同時開展了入水速度60、80、100 m/s下的數值模擬。在流固耦合數值模擬前，先比較剛體和變形體底部中心壓力和速度的變化情況。圖15為剛體和變形體底部中心點壓力曲線。表1為不同入水速度下圓柱底部中心點的壓力峰值，相比剛體結果，流固耦合方法的回轉體底面中心點壓力曲線的峰值較小，且隨著時間會有明顯的脈動現象。這是因為，彈性體圓柱入水時其底部觸水部分會發生變形來起到緩沖作用。因為緩沖變形呈反復狀態，彈性回轉體的壓力曲線產生脈動，一段時間后，隨著變形量的減少，脈動幅值也逐漸減小。由圖16可見：初始時刻，彈性回轉體速度和剛體回轉體速度變化基本相同；隨后，彈性回轉體在入水過程中產生變形來緩沖受到的沖擊力，反映在速度曲線上，彈性回轉體速度衰減將慢于剛體回轉體。

圖15 入水速度80 m/s時回轉體底部中心點壓力曲線Fig.15 Pressure curve of bottom central point at v0=80 m/s

表1 不同入水速度的壓力峰值Table 1 Peak pressures at different initial water entry velocities

圖16 速度80 m/s入水時速度曲線Fig.16 Velocity curves at v0=80 m/s

圖17為采用流固耦合方法的回轉體以不同速度入水的歸一化速度衰減曲線。在初始時刻，由于受到撞水瞬間產生的巨大沖擊力，回轉體速度衰減劇烈且速度越大衰減程度越劇烈；隨后，由于沖擊力迅速降低，在短時間內回轉體速度衰減程度放緩且衰減程度仍與入水速度有關，入水速度越大衰減程度越大；最后，在阻力、慣性力和重力的聯合作用下，速度呈類似線性下降，速度衰減狀態趨于穩定且不同入水速度的衰減程度近乎一致。

圖18為不同入水速度下回轉體底部中心單元的等效應力曲線。回轉體以不同速度入水，其底部中心單元的等效應力變化趨勢保持一致，入水瞬間回轉體應力急劇升高，形成瞬時峰值，隨后應力峰迅速下降，最后呈震蕩衰減狀態。應力變化趨勢與入水速度有關，等效應力峰值與入水速度平方相關，入水速度越大其峰值出現時間越早，且應力衰減震蕩現象越明顯，衰減幅度越大。回轉體100 m/s入水時，底部中心單元等效應力峰值為596.4 MPa，略微超出材料的屈服應力595 MPa，回轉體底部產生了小幅度的塑性應變。由此可見，在入水速度超過100 m/s時，該回轉體模型會導致結構產生塑性形變甚至失效，所以在設計回轉體模型時，還應該考慮入水瞬間的沖擊力作用下的結構強度問題。

圖19為不同入水速度下流場不同時刻的空泡形態曲線。在4 ms時，空泡處于敞開階段，空泡后部分呈外凸狀態，空泡完全包裹回轉體；在8 ms時，空泡充分發展，且入水速度越高空泡長度越長，半徑也越大；在12 ms時，空泡發展完畢，空泡徑向尺寸開始收縮，隨著入水位移加深，在水壓力和空泡表面張力作用下，入水速度80、100 m/s的回轉體開始出現頸縮現象，60 m/s的回轉體頸縮現象明顯；在16 ms時，空泡下半部分直徑較小，上半部分直徑較大，因為回轉體速度已經大幅度降低，回轉體底部附近流體獲得的動能較低，可見60 m/s的回轉體即將發生表面閉合，而80、100 m/s的回轉體表面閉合的時間更晚。因為入水速度較高時，在徑向水表面會獲得較高的擴張速度，直徑也會擴展更大，表面閉合在壓差作用下向軸向收緊所需要的時間也更長。

圖17 歸一化的速度衰減曲線Fig.17 Normalized velocity attenuation curves

圖18 回轉體中心單元的有效應力曲線Fig.18 Equivalent stress curves of central elements

圖19 空泡形態對比Fig.19 Comparisons of cavity features

3 結論

采用雙向流固耦合方法，對回轉體傾斜入水過程中結構載荷特性和流場演變進行了數值模擬，得到了如下結論。

（1）回轉體傾斜高速入水過程中，載荷先集中在觸水部分邊緣處，后集中于底部中心處，因空泡表面閉合后收縮，會產生載荷脈動峰值。

（2）采用流固耦合方法，轉體入水載荷結果小于剛體計算的，因彈性體圓柱入水時其底部觸水部分會發生變形起到緩沖作用，而緩沖變形的反復，使彈性回轉體的載荷曲線產生波動。

（3）回轉體撞水后，軸向空泡先呈現不對稱形態，隨著入水時間延長，空泡不對稱性變弱；入水速度60 m/s下，空泡先發生表面閉合，回轉體入水越快，空泡表面閉合越晚。

（4）回轉體入水載荷變化與入水速度有關，入水速度越大其峰值出現越早，且震蕩越明顯，速度超過100 m/s時回轉體會產生塑性形變，在結構設計上需考慮回轉體的結構強度。

本文中，入水速度遠低于聲速，未來將對亞聲速和超聲速下的入水問題展開流固耦合數值模擬，分析流體可壓縮性對入水過程的影響；同時，可進一步考慮結構的模態、材料參數和結構參數等對高速入水過程中結構動響應的影響及結構的損傷破壞情況。期望為結構物的強度設計提供參考。