楊超 姜宇 吳志剛

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100083)

水面滑跳是結(jié)構(gòu)體近水面小角度斜撞水過程中發(fā)生的一種特殊運(yùn)動(dòng),日常生活中在河邊進(jìn)行的“打水漂游戲” 就是一種典型的滑跳現(xiàn)象[1]。結(jié)構(gòu)體在滑跳運(yùn)動(dòng)過程中,不斷地以較大的水平速度和較小的垂向速度斜向撞擊水面,損失一部分動(dòng)能的同時(shí),還會(huì)被水面彈起。 近年來,不斷有學(xué)者提出利用這種滑跳運(yùn)動(dòng)來實(shí)現(xiàn)近水面高速機(jī)動(dòng)飛行的新概念跨介質(zhì)飛行器[2-5],也稱為擊水式飛行器。 因此,水面滑跳問題作為跨介質(zhì)飛行器的基礎(chǔ)問題也逐漸受到關(guān)注。

結(jié)構(gòu)體近水面滑跳問題涉及復(fù)雜多相流動(dòng)與結(jié)構(gòu)體大幅運(yùn)動(dòng)的相互作用,尤其是結(jié)構(gòu)體高速斜入水沖擊過程中帶來的受力劇烈變化。 目前,國內(nèi)外對(duì)結(jié)構(gòu)物高速斜入水沖擊載荷問題的研究,主要集中在結(jié)構(gòu)體垂向砰擊水面帶來的載荷,在工程應(yīng)用領(lǐng)域具有一定的局限性。 例如,Zhao和Faltinsen[6]在1993 年針對(duì)二維楔形體的垂直入水進(jìn)行了計(jì)算;樊征和劉瑩瑩[7]采用解析法分析了二維小升角楔形體的垂向入水問題;2008年,許國冬等[8]討論了V 形楔多種入水形式的相似解;Faltinsen 和Chezhian[9]以非軸對(duì)稱的船體模型入水為研究對(duì)象,基于Wagner 理論,研究了三維垂向入水砰擊的數(shù)值方法。 也有部分研究考慮了結(jié)構(gòu)體的彈性效應(yīng)。 例如,1991 年,顧懋祥等[10]采用迭代法計(jì)算了平頭旋轉(zhuǎn)殼體垂直入水時(shí)的水彈性效應(yīng);Takagi[11]在1994 年運(yùn)用Wagner 的平板近似理論計(jì)算了彈性結(jié)構(gòu)的入水沖擊問題;2017—2020 年,張?jiān)狼嗟萚12-13]對(duì)楔形體和弧形體的垂直入水進(jìn)行了試驗(yàn)和商用軟件數(shù)值仿真;2018 年,Mohammad 和Maurizio[14]對(duì)柔性細(xì)長(zhǎng)體的垂直入水沖擊載荷進(jìn)行了理論分析和試驗(yàn)驗(yàn)證;2004—2006 年,Wu 等[15-16]提出了二維楔形體的迭代求解方法,并運(yùn)用二維邊界元法,使用柯西積分計(jì)算了楔形體入水的相似解。 邊界元法是基于邊界積分方程,采用與有限元法[17]類似的分元、離散思想建立起來的。 由于邊界積分方程是定義在邊界上的,邊界元法只要在邊界上劃分單元,相比于區(qū)域解法降低了維度,具備較高的計(jì)算效率和計(jì)算精度[18]。 相對(duì)而言,結(jié)構(gòu)體高速斜入水沖擊問題的研究較少。 2008 年,王永虎、石秀華等[19]以入水沖擊的理論動(dòng)力學(xué)和入水彈道學(xué)理論為基礎(chǔ),完成了魚雷頭部尖拱的斜入水高速?zèng)_擊載荷的分析;2009 年,陳占暉和盧永錦[20]利用流固耦合有限元法研究了立方體斜向沖擊水面的運(yùn)動(dòng)特性;2012 年,鄭金偉和宗智[21]利用LSDYNA 軟件計(jì)算了三維剛性橢圓頭結(jié)構(gòu)的高速斜入水沖擊載荷;2014 年,Scolan[22]完成了三維不規(guī)則外形的剛體結(jié)構(gòu)低速斜入水沖擊載荷試驗(yàn);2021 年,蔡維軒等[23]完成了無旋圓盤的滑跳運(yùn)動(dòng)試驗(yàn);同年,王聰?shù)萚24]對(duì)旋轉(zhuǎn)圓盤的滑跳特性進(jìn)行了理論仿真。 可以看出,考慮結(jié)構(gòu)彈性效應(yīng)影響的三維結(jié)構(gòu)體高速斜入水沖擊載荷數(shù)值求解算法的研究較少。

為了高效求解三維結(jié)構(gòu)體滑跳運(yùn)動(dòng)問題,本文基于邊界元法,建立了一種可以考慮彈性效應(yīng)的三維結(jié)構(gòu)體近水面滑跳運(yùn)動(dòng)的時(shí)域數(shù)值仿真方法。 在三維結(jié)構(gòu)體的滑跳運(yùn)動(dòng)過程中,高速斜入水沖擊載荷由邊界元法求得,利用狀態(tài)方程直接積分法(direct integration of state equation method)求解得到結(jié)構(gòu)節(jié)點(diǎn)位移動(dòng)力響應(yīng),并通過滑跳運(yùn)動(dòng)動(dòng)力學(xué)模型更新結(jié)構(gòu)體的質(zhì)心位置及運(yùn)動(dòng)速度。 這使得本文方法既可以應(yīng)用于三維彈性結(jié)構(gòu)體,也擁有較高的計(jì)算效率。 針對(duì)剛球高速斜入水沖擊現(xiàn)象進(jìn)行了數(shù)值仿真,并與試驗(yàn)值進(jìn)行了對(duì)比,驗(yàn)證了該方法流體載荷求解部分的準(zhǔn)確性。針對(duì)三維球冠體的近水面滑跳運(yùn)動(dòng)進(jìn)行了時(shí)域數(shù)值仿真,得到了對(duì)應(yīng)的運(yùn)動(dòng)特性。 研究了不同的初始參數(shù)(結(jié)構(gòu)質(zhì)量、初始高度、水平拋出速度和半徑大小)和結(jié)構(gòu)彈性效應(yīng)對(duì)球冠體近水面滑跳運(yùn)動(dòng)的影響,總結(jié)了相應(yīng)的規(guī)律。

1 問題描述

本文分析的三維結(jié)構(gòu)體近水面滑跳運(yùn)動(dòng)問題可以由圖1 表示,半徑為r0的球冠體結(jié)構(gòu)在質(zhì)心距離水面高度H0的位置上以速度V0水平拋出。圖中采用了如下坐標(biāo)系描述:流場(chǎng)坐標(biāo)系O-xyz,原點(diǎn)固定在流場(chǎng)正中心的O點(diǎn),x軸和y軸固定在與O點(diǎn)同一高度的水平面上,z軸垂直向上。 結(jié)構(gòu)坐標(biāo)系O′-x′y′z′,原點(diǎn)固定在球冠體的質(zhì)心O′點(diǎn),x′軸和y′軸固定在與O′點(diǎn)同一高度的水平面上,z′軸垂直向上。 圖1 中,Ω為流體域,S為流體域邊界。

圖1 球冠體近水面滑跳運(yùn)動(dòng)問題示意圖Fig.1 Schematic diagram of near-water-surface skipping motion of hemispherical structure

2 理論方法

三維結(jié)構(gòu)體高速斜入水沖擊是流體和結(jié)構(gòu)相互作用的強(qiáng)非線性過程。 在很多情況下,流體運(yùn)動(dòng)的速度比聲速小得多,因此,除平頭物體垂直入水外,一般可以忽略流體的可壓縮性[25]。 入水沖擊載荷具有瞬態(tài)極大值的特點(diǎn),與之相比,入水沖擊時(shí)流體的重力因素所造成的影響較小,因此作為簡(jiǎn)化條件,忽略流體重力的影響。

本文建立的基于邊界元法的三維彈性結(jié)構(gòu)體近水面滑跳運(yùn)動(dòng)時(shí)域迭代仿真方法的基本原理如圖2 所示。 三維彈性結(jié)構(gòu)體在近水面滑跳過程中受到斜入水沖擊載荷、摩擦阻力、重力和浮力的共同作用,每個(gè)時(shí)間步內(nèi)流場(chǎng)作用在三維結(jié)構(gòu)體上的斜入水沖擊載荷由邊界元法求解得到,利用狀態(tài)方程直接積分法求解得到結(jié)構(gòu)節(jié)點(diǎn)位移動(dòng)力響應(yīng)。 流場(chǎng)載荷連同重力和浮力代入到滑跳運(yùn)動(dòng)動(dòng)力學(xué)模型中進(jìn)行分析,得到結(jié)構(gòu)體新的質(zhì)心位置及運(yùn)動(dòng)速度。 重新構(gòu)造流場(chǎng)的形狀和邊界條件,并進(jìn)入下一步迭代計(jì)算。

圖2 時(shí)域迭代求解方法基本原理Fig.2 Basic principle of time domain iterative solution method

2.1 邊界元法

根據(jù)流體區(qū)域的簡(jiǎn)化和假設(shè),流體為無旋、無黏、不可壓流體,則三維流體運(yùn)動(dòng)的速度勢(shì)φ滿足如下Laplace 方程:

為了確定方程(1)在空間某個(gè)內(nèi)部域或外部域的解,還必須附加定解條件,即邊界條件。 對(duì)于本文分析的問題,滿足Laplace 方程的第三邊界條件,即混合條件。 流體區(qū)域自由表面上的邊界條件為

對(duì)于三維Laplace 方程,其基本解φ*應(yīng)滿足條件[15]:

若以φ為所求三維Laplace 方程邊值問題的解,而φ*作為輔助解ψ代入Green 第二公式(由Gauss 公式推得):

式(9)中的積分是對(duì)dS(q)積分,即在邊界面對(duì)場(chǎng)點(diǎn)進(jìn)行積分,q代表邊界面上的場(chǎng)點(diǎn)。

為了建立邊界積分方程,只需將基本解的域內(nèi)源點(diǎn)P從流體域內(nèi)趨向于邊界面源點(diǎn)p,于是得到

式中:α為邊界面上源點(diǎn)p處的立體角系數(shù)。

對(duì)于本文分析的三維結(jié)構(gòu)體,進(jìn)行常值三角形單元的邊界劃分,如圖3 所示。 對(duì)于常值三角形單元來說,相鄰單元間未知量本身都不連續(xù)。在每個(gè)邊界元上φ和?φ/?n都是常量,并用該邊界元中點(diǎn)的值來代替。

圖3 常值三角形單元Fig.3 Constant triangular element

p點(diǎn)取遍所有邊界單元的中點(diǎn),就可以得到如下矩陣:

式中:H和G分別為φ和?φ/?n的積分系數(shù)矩陣,僅僅與求解域的幾何特性有關(guān)。

通過求解式(11),就可以得到三維結(jié)構(gòu)體的入水過程中任意邊界元中點(diǎn)的速度勢(shì)。 在求得邊界單元中點(diǎn)處的勢(shì)函數(shù)值后,可以用一種簡(jiǎn)便的數(shù)值微分方法[26-27]計(jì)算浸濕表面上速度分布。設(shè)浸濕表面方程為y=y(x,z),勢(shì)函數(shù)在浸濕表面上的增量可以寫成

式中:u、v、w分別為邊界元中點(diǎn)速度的3 個(gè)分量;l、m、k分別為邊界元中點(diǎn)單位法向量的3 個(gè)分量。

圖4 中,qi為控制點(diǎn),Vi為控制點(diǎn)速度向量,Γ為控制點(diǎn)單位法向量,i和j分別為x坐標(biāo)方向和z坐標(biāo)方向的單位向量。

圖4 速度分量求解示意圖Fig.4 Schematic diagram of velocity component solution

在計(jì)算得到每個(gè)邊界單元中心點(diǎn)的速度和速度勢(shì)后,就可以利用下述改進(jìn)的非線性伯努利方程求得每個(gè)邊界單元對(duì)應(yīng)的水動(dòng)力壓強(qiáng)[28]:

式中:pi為邊界單元的壓強(qiáng);ρ為流體密度。

最后,每個(gè)邊界單元的壓強(qiáng)乘以每個(gè)邊界單元的面積,即可得到每個(gè)邊界單元對(duì)結(jié)構(gòu)體的作用力。

2.2 高速斜入水沖擊載荷模型

由于浸沒在水中的體積排擠流體產(chǎn)生液面隆起現(xiàn)象[29],又稱Wagner 型入水沖擊,如圖5 所示。 取決于入水結(jié)構(gòu)體的外形和入水角等初始狀態(tài),下面給出確定液面隆起等效位置He的表達(dá)式,即

式中:θ為入水傾角;t為時(shí)間;η(x,y,z,θ,t)為函數(shù)表達(dá)式。 參考文獻(xiàn)[30],圓球體傾斜入水時(shí),He的取值固定為1.35 倍的初始液面入水深度。考慮自由液面的隆起現(xiàn)象后,會(huì)導(dǎo)致結(jié)構(gòu)體浸濕表面擴(kuò)大,高速斜入水沖擊載荷增大。

結(jié)構(gòu)體高速斜入水沖擊載荷可以分為較大的水平速度帶來的劃水力和較小的垂向速度帶來的垂向入水砰擊載荷2 部分。 當(dāng)結(jié)構(gòu)體高速斜入水時(shí),流體被排開,因此球冠體被浸濕的表面只有與速度矢量V夾角小于90°且低于液面隆起等效位置He的部分[31-32],示意圖如圖5 所示。 同樣,只有這部分浸濕表面會(huì)受到高速斜入水沖擊載荷。

圖5 浸濕表面示意圖Fig.5 Schematic diagram of water soaked surface

對(duì)于結(jié)構(gòu)體浸濕表面的三角形常值單元,計(jì)算高速斜入水沖擊載荷時(shí)的邊界條件為結(jié)構(gòu)體水平運(yùn)動(dòng)速度Vx的法向分量、垂向運(yùn)動(dòng)速度Vz的法向分量與結(jié)構(gòu)體旋轉(zhuǎn)角速度ω帶來的速度Vω的法向分量之和,示意圖如圖6 所示。

圖6 斜入水沖擊載荷模型單元邊界條件Fig.6 Element boundary conditions of oblique water impact load model

當(dāng)使用邊界元法計(jì)算斜入水沖擊載荷時(shí),假設(shè)流體為理想流體。 而實(shí)際流體是有黏性的,因此在斜入水沖擊載荷模型中,結(jié)構(gòu)體與流體相接觸的表面受到摩擦阻力的作用。 摩擦阻力系數(shù)cf可以直接選用帕蘭特-許立汀紊流摩阻計(jì)算公式[3,33],即

式中:Re為雷諾數(shù)。

對(duì)于球冠體來說,結(jié)構(gòu)底部較為平順,因此流體黏性帶來的摩擦阻力主要作用于水平方向,后續(xù)的分析中忽略垂直方向的摩擦阻力。

2.3 結(jié)構(gòu)動(dòng)力響應(yīng)分析

本文采用狀態(tài)方程直接積分法計(jì)算每個(gè)時(shí)間步內(nèi),彈性結(jié)構(gòu)在流場(chǎng)載荷下的彈性動(dòng)力響應(yīng)[34]。 結(jié)構(gòu)體在模態(tài)向量坐標(biāo)下的動(dòng)力平衡方程[35-36]可以表示為

式中:

當(dāng)有限元模型節(jié)點(diǎn)處所受的外力矩陣確定時(shí),就可以通過狀態(tài)方程直接積分法計(jì)算得到三維彈性結(jié)構(gòu)體的位移響應(yīng)[37]。

2.4 滑跳運(yùn)動(dòng)動(dòng)力學(xué)模型

進(jìn)行滑跳運(yùn)動(dòng)的結(jié)構(gòu)體一般為對(duì)稱結(jié)構(gòu),因此只需要考慮結(jié)構(gòu)體的縱向平面內(nèi)運(yùn)動(dòng):水平運(yùn)動(dòng)、垂直運(yùn)動(dòng)和俯仰運(yùn)動(dòng)。 通過動(dòng)量定理和動(dòng)量矩定理[38],可以建立滑跳運(yùn)動(dòng)的全空間運(yùn)動(dòng)方程組。

當(dāng)結(jié)構(gòu)體在空中飛行時(shí),沒有流體載荷的作用,水平方向的運(yùn)動(dòng)方程為

式中:Flzi為浸濕表面單元高速斜入水沖擊載荷的升力分量;Ffz為摩擦阻力的垂直方向分量。

俯仰方向上的運(yùn)動(dòng)方程為

式中:Pc為結(jié)構(gòu)體質(zhì)心位置;Pi為浸濕表面單元中點(diǎn)位置;rx(Pc,Pi) 為結(jié)構(gòu)體質(zhì)心與浸濕表面單元中點(diǎn)x坐標(biāo)的差值;rz(Pc,Pi) 為結(jié)構(gòu)體質(zhì)心與浸濕表面單元中點(diǎn)z坐標(biāo)的差值。

2.5 時(shí)域迭代求解方法

本文建立的基于邊界元法的三維結(jié)構(gòu)體近水面滑跳運(yùn)動(dòng)時(shí)域迭代仿真方法流程如圖7 所示。首先,進(jìn)行分析模型的初始化和分析參數(shù)的設(shè)置,開始時(shí)域迭代計(jì)算。 在每一個(gè)時(shí)間步內(nèi),對(duì)結(jié)構(gòu)體是否觸水進(jìn)行判斷。 如果觸水,由該時(shí)刻下的結(jié)構(gòu)體表面和初始流場(chǎng)進(jìn)行組合構(gòu)造,確定該時(shí)間步的流場(chǎng)形狀及流場(chǎng)邊界條件,并提交邊界元法程序進(jìn)行求解,得到該時(shí)間步內(nèi)流場(chǎng)對(duì)結(jié)構(gòu)體的高速斜入水沖擊載荷。 然后,將流場(chǎng)載荷插值到結(jié)構(gòu)浸濕表面節(jié)點(diǎn),并利用狀態(tài)方程直接積分法得到結(jié)構(gòu)節(jié)點(diǎn)的彈性位移動(dòng)力響應(yīng)。 如果結(jié)構(gòu)體沒有觸水,則該時(shí)間步內(nèi)流場(chǎng)對(duì)結(jié)構(gòu)體的高速斜入水沖擊載荷為0。 最后,根據(jù)滑跳運(yùn)動(dòng)動(dòng)力學(xué)模型,更新結(jié)構(gòu)體的位置與速度信息,并進(jìn)入下一步迭代計(jì)算。

圖7 時(shí)域迭代求解方法流程Fig.7 Flow chart of time domain iterative solution method

3 仿真計(jì)算及分析

3.1 數(shù)值方法驗(yàn)證

1976 年,為了研究高速斜入水沖擊載荷的數(shù)值,Soliman 等[39]完成了鋼球高速斜入水沖擊試驗(yàn),得到了一系列相關(guān)數(shù)據(jù)。

在Soliman 的試驗(yàn)中,高速斜入水沖擊結(jié)構(gòu)為鋼制實(shí)心球體,直徑為25.4 mm,入水速度大小為64.62 m/s。 試驗(yàn)使用安裝在鑄鐵平臺(tái)上的彈藥槍,將炮彈射入裝滿水的有機(jī)玻璃水箱(6 ft長(zhǎng)、1 ft 寬、1 ft 深,1 ft(英尺) =0.304 8 m)。 相關(guān)參數(shù)采用35 mm 攝像機(jī)及頻閃儀進(jìn)行測(cè)量。試驗(yàn)的示意圖如圖8 所示。 圖中的鋼球以入水速度V、入水傾角θ高速?zèng)_擊水面,c代表鋼球質(zhì)心。

圖8 鋼球高速斜入水試驗(yàn)示意圖Fig.8 Schematic diagram of high speed oblique water entry test of steel sphere

為了驗(yàn)證本文提出的三維結(jié)構(gòu)體高速斜入水沖擊載荷時(shí)域迭代求解方法的可行性和準(zhǔn)確性,針對(duì)Soliman 的試驗(yàn)進(jìn)行了數(shù)值仿真。 鋼球具有較大的水平運(yùn)動(dòng)速度及較小的垂向運(yùn)動(dòng)速度。 數(shù)值計(jì)算結(jié)果與試驗(yàn)值的對(duì)比如圖9 和圖10所示。 其中,鋼球離開水面時(shí)運(yùn)動(dòng)速度的數(shù)值計(jì)算結(jié)果與試驗(yàn)值對(duì)比如圖9 所示;鋼球離開水面時(shí)出水傾角的數(shù)值計(jì)算結(jié)果與試驗(yàn)值對(duì)比如圖10 所示。

圖9 離開水面速度數(shù)值計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig.9 Comparison between calculated and experimental values of velocity leaving water surface

圖10 離開水面傾角數(shù)值計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig.10 Comparison between calculated and experimental values of dip angle leaving water surface

通過圖9 和圖10 可知,當(dāng)鋼球入水傾角θ增大時(shí),鋼球在觸水滑跳過程中的動(dòng)能衰減增大,其離開水面時(shí)的速度傾角增大。 此外還可以看出,本文方法的數(shù)值計(jì)算結(jié)果與試驗(yàn)數(shù)值吻合良好,驗(yàn)證了本文提出的三維結(jié)構(gòu)體高速斜入水沖擊載荷時(shí)域迭代求解方法的可行性和準(zhǔn)確性。

3.2 仿真建模

在本文方法中,每一個(gè)時(shí)間步都需要用到結(jié)構(gòu)體的浸濕表面和未受擾動(dòng)的流場(chǎng)來對(duì)下一個(gè)時(shí)間步的流場(chǎng)進(jìn)行重構(gòu),結(jié)構(gòu)動(dòng)力響應(yīng)分析也需要將結(jié)構(gòu)體的模態(tài)分析結(jié)果作為原始數(shù)據(jù)。 因此,需要提前對(duì)結(jié)構(gòu)體進(jìn)行有限元建模和模態(tài)分析,建立未受擾動(dòng)的流場(chǎng),設(shè)置初始參數(shù)。

為了方便構(gòu)成流場(chǎng)邊界,需要預(yù)先對(duì)球冠體進(jìn)行網(wǎng)格劃分,即有限元建模。 球冠體半徑為0.152 4 m,材料密度為1 800 kg/m3,質(zhì)量為1.554 kg。 在結(jié)構(gòu)有限元建模時(shí),球冠體表面均采用三角形單元進(jìn)行網(wǎng)格劃分,構(gòu)建完成的結(jié)構(gòu)有限元模型如圖11 所示。 材料的彈性模量為70 GPa,剪切模量為26 GPa,泊松比為0. 33。通過實(shí)心球冠體的模態(tài)分析可以得到,一階固有頻率高達(dá)3 800 Hz,球冠體的彈性效應(yīng)很小。 因此,實(shí)心球冠體可近似地看作一個(gè)剛性物體。

圖11 球冠體有限元模型Fig.11 Finite element model of spherical crown

針對(duì)直徑為0.304 8 m 的球冠體結(jié)構(gòu),建立一個(gè)邊長(zhǎng)為1.0 m 的立方體流場(chǎng),初始未受擾動(dòng)的流場(chǎng)如圖12 所示。 由于采用邊界元法進(jìn)行高速斜入水沖擊載荷分析,流場(chǎng)內(nèi)部沒有網(wǎng)格劃分,只是在邊界表面上劃分了共10 682 個(gè)三角形單元。 在滑跳運(yùn)動(dòng)計(jì)算過程中,流場(chǎng)隨結(jié)構(gòu)體質(zhì)心的水平運(yùn)動(dòng)進(jìn)行平移,流場(chǎng)的坐標(biāo)原點(diǎn)始終處于質(zhì)心的正下方。

圖12 初始未受擾動(dòng)的流場(chǎng)Fig.12 Initial undisturbed flow field

對(duì)于標(biāo)準(zhǔn)算例,球冠體在距離水平面1.0 m的高度處,以25.0 m/s 的速度被水平拋出。 時(shí)域仿真總時(shí)長(zhǎng)為1.4 s,迭代步長(zhǎng)為0.4 ms,液體的密度設(shè)置為1 000 kg/m3,重力加速度g為9.8 m/s2。從第一個(gè)時(shí)間步開始,每個(gè)時(shí)間步的流場(chǎng)域上表面都是由流場(chǎng)的自由表面和結(jié)構(gòu)浸濕表面重構(gòu)而成。

3.3 剛性球冠體滑跳運(yùn)動(dòng)仿真

根據(jù)3.2 節(jié)的標(biāo)準(zhǔn)算例參數(shù)設(shè)置,應(yīng)用本文建立的數(shù)值仿真方法對(duì)剛性球冠體進(jìn)行近水面滑跳運(yùn)動(dòng)時(shí)域仿真,仿真總時(shí)長(zhǎng)為1.4 s,仿真結(jié)果如圖13 ～圖19 所示。 其中,圖13 為球冠體近水面滑跳運(yùn)動(dòng)的軌跡仿真結(jié)果,圖14 為球冠體近水面滑跳運(yùn)動(dòng)的垂向速度時(shí)歷曲線,圖16 為球冠體近水面滑跳運(yùn)動(dòng)的水平速度時(shí)歷曲線。

圖13 滑跳運(yùn)動(dòng)軌跡仿真結(jié)果Fig.13 Simulation results of skipping motion trajectory

圖14 滑跳運(yùn)動(dòng)垂向速度時(shí)歷曲線Fig.14 Vertical velocity time history curve of skipping motion

圖15 第1 次觸水垂向速度時(shí)歷曲線Fig.15 Vertical velocity time history curve in the first contact with water

圖16 滑跳運(yùn)動(dòng)水平速度時(shí)歷曲線Fig.16 Horizontal velocity time history curve of skipping motion

根據(jù)圖13 可知,本文建立的基于邊界元法的三維結(jié)構(gòu)體近水面滑跳運(yùn)動(dòng)時(shí)域數(shù)值仿真方法適用于球冠體近水面滑跳運(yùn)動(dòng)的數(shù)值仿真。 球冠體在1.4 s 內(nèi)經(jīng)歷了3 次觸水滑跳運(yùn)動(dòng),每次觸水彈起的質(zhì)心最高位置逐漸降低。 此外還可以看出,第1 次觸水滑跳過程時(shí)長(zhǎng)約24 ms,入水深度為0. 025 m;第2 次觸水滑跳過程時(shí)長(zhǎng)約為28 ms,入水深度為0. 018 m;第3 次觸水滑跳過程時(shí)長(zhǎng)約32 ms,入水深度為0. 014 m。 可以看出,隨著觸水滑跳次數(shù)的增加,觸水滑跳過程的時(shí)長(zhǎng)增加,但入水深度減小。 通過圖14 和圖16 可知,球冠體在結(jié)構(gòu)未觸水的時(shí)間段內(nèi)只受到重力作用,下降速度勻速增大,水平速度不變;當(dāng)球冠體觸水時(shí),受到流場(chǎng)力的作用,在極短的時(shí)間內(nèi)垂向速度由負(fù)轉(zhuǎn)正,水平速度減小。 此外還可以看出,球冠體每次觸水彈起后,垂向速度的最大值都會(huì)發(fā)生衰減。 垂向速度最大值和水平速度的減小,代表了球冠體動(dòng)能在不斷的觸水滑跳過程中被消耗。 滑跳運(yùn)動(dòng)過程中結(jié)構(gòu)體的水平速度不斷下降,當(dāng)觸水彈起高度不足以使結(jié)構(gòu)體脫離水面之后,滑跳運(yùn)動(dòng)結(jié)束。

由于球冠體每次觸水滑跳時(shí)間過短,體現(xiàn)在圖14 和圖16 中,觸水滑跳的垂向速度變化和水平速度變化近似折線形式。 而事實(shí)上,觸水彈跳的過程中速度變化是連續(xù)的。 球冠體第1 次觸水滑跳過程中的垂向速度變化曲線和水平速度變化曲線分別如圖15 和圖17 所示。 通過圖15 可知,第1 次觸水滑跳時(shí)間歷程為0.416 ～0.44 s,時(shí)長(zhǎng)為24 ms。 當(dāng)球冠體剛接觸水面時(shí),垂向速度約為-4.0 m/s,受到流場(chǎng)力的作用,垂向速度迅速變化,約為0. 426 s 時(shí)由負(fù)轉(zhuǎn)正,并繼續(xù)增加至2.2 m/s離開水面。 通過圖17 可知,水平速度由剛接觸水面時(shí)的25 m/s 持續(xù)減小至23.5 m/s 離開水面,衰減率為6%(衰減率λ= (Vin-Vout)/Vin,其中Vin為球冠體第1 次觸水之前的水平速度,Vout為球冠體第1 次觸水之后的水平速度,下同)。

圖17 第1 次觸水水平速度時(shí)歷曲線Fig.17 Horizontal velocity time history curve in the first contact with water

球冠體第1 次觸水過程中的高速斜入水沖擊載荷仿真結(jié)果如圖18 和圖19 所示。 其中,圖18為球冠體在第1 次觸水過程中的高速斜入水沖擊載荷垂直方向分量時(shí)歷圖,圖19 為球冠體在第1次觸水過程中的高速斜入水沖擊載荷水平方向分量時(shí)歷圖。

圖18 第1 次觸水垂向受力時(shí)歷曲線Fig.18 Vertical stress time history curve in the first contact with water

圖19 第1 次觸水水平受力時(shí)歷曲線Fig.19 Horizontal stress time history curve in the first contact with water

通過圖18 和圖19 可以看出,在球冠體的第1 次觸水過程中,高速斜入水沖擊載荷提供了阻力和升力,都是很快達(dá)到峰值并逐漸減小。 分開來看,垂直方向受到的載荷方向向上,起到升力的作用,數(shù)值相對(duì)較大,在0. 422 s 時(shí)達(dá)到峰值887.2 N,在其作用下,球冠體下降速度迅速減為0 并加速上升,直至離開水面;水平方向持續(xù)受到阻力的作用,數(shù)值相對(duì)較小,0.422 s 時(shí)達(dá)到峰值-216.0 N,在其作用下水平運(yùn)動(dòng)速度不斷減小。

對(duì)于球冠體而言,直徑為0.304 8 m,第1 次觸水滑跳時(shí)運(yùn)動(dòng)速度為25 m/s,流體密度為1 000 kg/m3。 經(jīng) 過 計(jì) 算, 此 時(shí) 的 雷 諾 數(shù) 高 達(dá)7.5 ×106。 由圖19 可知,在球冠體第1 次觸水滑跳的過程中,考慮流體黏性的高速斜入水沖擊載荷水平方向分量峰值為-216.0 N,未考慮流體黏性的高速斜入水沖擊載荷水平方向分量峰值為-205.5 N,影響約為5. 1%。 由此可見,在高速斜入水沖擊仿真過程中,流體黏性帶來的影響較小,但是也不可忽略。 本文后續(xù)的仿真中,都已考慮流體黏性帶來的摩擦阻力影響。

3.4 變參分析

3.4.1 結(jié)構(gòu)質(zhì)量對(duì)滑跳運(yùn)動(dòng)的影響

為了得到結(jié)構(gòu)質(zhì)量對(duì)球冠體滑跳運(yùn)動(dòng)的影響規(guī)律,在不改變其他分析參數(shù)的情況下,改變球冠體的密度,使得結(jié)構(gòu)質(zhì)量m0分別為1. 0,1.554,3.0 kg,仿真總時(shí)長(zhǎng)為1.3 s。 不同結(jié)構(gòu)質(zhì)量的球冠體滑跳軌跡對(duì)比如圖20 所示,不同結(jié)構(gòu)質(zhì)量的滑跳運(yùn)動(dòng)垂向速度時(shí)歷對(duì)比如圖21 所示,不同結(jié)構(gòu)質(zhì)量的滑跳運(yùn)動(dòng)水平速度時(shí)歷對(duì)比如圖22所示。

圖20 不同結(jié)構(gòu)質(zhì)量的滑跳軌跡對(duì)比Fig.20 Comparison of skipping motion trajectory with different structural quality

圖21 不同結(jié)構(gòu)質(zhì)量的滑跳垂向速度時(shí)歷對(duì)比Fig.21 Time history comparison of skipping motionvertical velocity with different structural quality

圖22 不同結(jié)構(gòu)質(zhì)量的滑跳水平速度時(shí)歷對(duì)比Fig.22 Time history comparison of skipping motion horizontal velocity with different structural quality

由圖20 和圖21 可知,當(dāng)結(jié)構(gòu)質(zhì)量增大時(shí),結(jié)構(gòu)體每次觸水滑跳入水深度增加,彈起速度和高度也會(huì)增大,導(dǎo)致每次滑跳落點(diǎn)之間的水平距離加大。 在相同的時(shí)間和距離內(nèi),隨著結(jié)構(gòu)質(zhì)量的增加,觸水滑跳的次數(shù)減少。 由圖22 可知,質(zhì)量為1.0 kg 時(shí),第1 次觸水滑跳后水平速度降為23.64 m/s,衰減率為5. 44%;質(zhì)量為1. 554 kg時(shí),第1 次觸水滑跳后水平速度降為23. 5 m/s,衰減率為6%;質(zhì)量為3.0 kg 時(shí),第1 次觸水滑跳后水平速度降為23. 22 m/s,衰減率為7. 12%。當(dāng)結(jié)構(gòu)質(zhì)量增大時(shí),每次觸水彈跳水平速度的下降量增大。

3.4.2 初始高度對(duì)滑跳運(yùn)動(dòng)的影響

為了得到初始高度對(duì)球冠體滑跳運(yùn)動(dòng)的影響規(guī)律,在不改變其他分析參數(shù)的情況下,改變初始高度H0分別為0. 5,1. 0,1. 5 m,仿真總時(shí)長(zhǎng)為1.3 s。 不同初始高度的滑跳軌跡對(duì)比如圖23 所示,不同初始高度的滑跳垂向速度時(shí)歷對(duì)比如圖24所示,不同初始高度的滑跳水平速度時(shí)歷對(duì)比如圖25 所示。

圖23 不同初始高度的滑跳軌跡對(duì)比Fig.23 Comparison of skipping motion trajectories with different initial heights

圖24 不同初始高度的滑跳垂向速度時(shí)歷對(duì)比Fig.24 Time history comparison of skipping motion vertical velocity with different initial heights

圖25 不同初始高度的滑跳水平速度時(shí)歷對(duì)比Fig.25 Time history comparison of skipping motion horizontal velocity with different initial heights

通過圖23 和圖24 可知,初始拋出高度增加,結(jié)構(gòu)體每次觸水滑跳入水深度增加,觸水彈起的高度和彈起速度也隨之增大,導(dǎo)致每次滑跳落點(diǎn)之間的水平距離加大,在同樣的時(shí)間和距離內(nèi)觸水滑跳的次數(shù)變少。 通過圖25 可知,初始高度為0. 5 m 時(shí),第1 次觸水滑跳后水平速度降為24.12 m/s,衰減率為3.52%;初始高度為1.0 m時(shí),第1 次觸水滑跳后水平速度降為23. 5 m/s,衰減率為6%;初始高度為1.5 m 時(shí),第1 次觸水滑跳后水平速度降為23. 0 m/s,衰減率為8%。初始拋出高度的增大,會(huì)導(dǎo)致每次觸水彈跳水平速度的下降量增加,動(dòng)能損耗增多。

3.4.3 水平拋出速度對(duì)滑跳運(yùn)動(dòng)的影響

為了得到水平拋出速度對(duì)球冠體滑跳運(yùn)動(dòng)的影響規(guī)律,在不改變其他分析參數(shù)的情況下,改變球冠體的水平拋出速度V0分別為20,25,30 m/s,仿真總時(shí)長(zhǎng)為1.3 s。 不同水平拋出速度的滑跳軌跡對(duì)比如圖26 所示,不同水平拋出速度的滑跳運(yùn)動(dòng)垂向速度時(shí)歷對(duì)比如圖27 所示。

通過圖26 和圖27 可以看出,隨著水平拋出速度的增大,相同時(shí)間內(nèi)球冠體的水平位移增加,每次觸水滑跳的彈起高度和彈起速度也會(huì)增加,相同時(shí)間和距離內(nèi)觸水滑跳的次數(shù)減少。

圖26 不同水平拋出速度的滑跳軌跡對(duì)比Fig.26 Comparison of skipping motion trajectories with different horizontal throw speeds

圖27 不同水平拋出速度的滑跳垂向速度時(shí)歷對(duì)比Fig.27 Time history comparison of skipping motion vertical velocity with different horizontal throw speeds

3.4.4 半徑大小對(duì)滑跳運(yùn)動(dòng)的影響

為了得到半徑大小對(duì)球冠體滑跳運(yùn)動(dòng)的影響規(guī)律,在不改變其他分析參數(shù)的情況下,改變球冠體的半徑大小r0分別為0.1,0.152 4,0.2 m,仿真總時(shí)長(zhǎng)為1.3 s。 不同半徑大小的滑跳軌跡對(duì)比如圖28 所示,不同半徑大小的滑跳運(yùn)動(dòng)垂向速度時(shí)歷對(duì)比如圖29 所示,不同半徑大小的滑跳運(yùn)動(dòng)水平速度時(shí)歷對(duì)比如圖30 所示。

圖28 不同半徑大小的滑跳軌跡對(duì)比Fig.28 Comparison of skipping motion trajectories with different radius sizes

圖29 不同半徑大小的滑跳垂向速度時(shí)歷對(duì)比Fig.29 Time history comparison of skipping motion vertical velocity with different radius sizes

圖30 不同半徑大小的滑跳水平速度時(shí)歷對(duì)比Fig.30 Time history comparison of skipping motion horizontal velocity with different radius sizes

通過圖28 和圖29 所示,隨著球冠體的半徑減小,結(jié)構(gòu)每次觸水彈跳的入水深度增加,每次觸水滑跳的高度和彈起速度也隨之增加,導(dǎo)致每次觸水滑跳運(yùn)動(dòng)之間飛行的距離增大,在同樣的時(shí)間或者距離內(nèi)觸水彈起的次數(shù)減少。 此外,通過圖30 可知,球冠體半徑為0.1 m 時(shí),第1 次觸水滑跳后水平速度降為22. 83 m/s,衰減率為8.68%;球冠體半徑為0.152 4 m 時(shí),第1 次觸水滑跳后水平速度降為23.5 m/s,衰減率為6%;球冠體半徑為0.2 m 時(shí),第1 次觸水滑跳后水平速度降為23.84 m/s,衰減率為4.64%,隨著球冠體半徑的減小,結(jié)構(gòu)每次觸水彈跳水平速度的下降量增大,動(dòng)能損耗增多。

3.5 彈性效應(yīng)的影響

為了得到結(jié)構(gòu)體彈性效應(yīng)對(duì)球冠體滑跳運(yùn)動(dòng)的影響規(guī)律,在不改變其他分析參數(shù)的情況下,考慮結(jié)構(gòu)體的彈性效應(yīng)并進(jìn)行時(shí)域仿真分析。 在考慮結(jié)構(gòu)彈性效應(yīng)的影響時(shí),發(fā)現(xiàn)實(shí)心高剛度材料的球冠體幾乎沒有彈性效應(yīng)。 為了放大結(jié)構(gòu)的彈性效應(yīng),將實(shí)心球冠體改為空心,殼體的厚度為1.0 mm,材料剛度同時(shí)降低,彈性模量改為0.7 MPa,泊松比為0.33。 進(jìn)行模態(tài)分析后得到,空心球冠體的一階彈性模態(tài)頻率為12.32 Hz,二階彈性模態(tài)頻率為14.91 Hz,仿真總時(shí)長(zhǎng)為1.1 s。 有無彈性效應(yīng)時(shí)結(jié)構(gòu)體第1 次觸水所受垂向流場(chǎng)載荷時(shí)歷對(duì)比如圖31 所示,有無彈性效應(yīng)時(shí)結(jié)構(gòu)體第1 次觸水所受水平方向流場(chǎng)載荷時(shí)歷對(duì)比如圖32所示,有無彈性效應(yīng)的滑跳軌跡對(duì)比如圖33 所示,有無彈性效應(yīng)的滑跳運(yùn)動(dòng)垂向速度時(shí)歷對(duì)比如圖34所示,有無彈性效應(yīng)的滑跳運(yùn)動(dòng)水平速度時(shí)歷對(duì)比如圖35 所示。 考慮彈性效應(yīng)時(shí),球冠體最低點(diǎn)的彈性變形時(shí)歷圖如圖36 所示。

圖31 有無彈性效應(yīng)的第1 次觸水垂向受力時(shí)歷對(duì)比Fig.31 Time history comparison of vertical stress in the first contact with water with or without elastic effect

圖32 有無彈性效應(yīng)的第1 次觸水水平受力時(shí)歷對(duì)比Fig.32 Time history comparison of horizontal stress in the first contact with water with or without elastic effect

圖33 有無彈性效應(yīng)的滑跳軌跡對(duì)比Fig.33 Comparison of skipping motion trajectories with or without elastic effect

圖34 有無彈性效應(yīng)的滑跳垂向速度時(shí)歷對(duì)比Fig.34 Time history comparison of skipping motion vertical velocity with or without elastic effect

圖35 有無彈性效應(yīng)的滑跳水平速度時(shí)歷對(duì)比Fig.35 Time history comparison of skipping motion horizontal velocity with or without elastic effect

圖36 球冠體最低點(diǎn)的彈性變形時(shí)歷曲線Fig.36 Time history curves of elastic deformation at the lowest point of spherical crown

通過圖31 和圖32 可知,考慮結(jié)構(gòu)體的彈性效應(yīng)之后,觸水滑跳過程中流體作用在結(jié)構(gòu)上的垂向載荷和水平方向載荷都有所減小,這也影響了結(jié)構(gòu)體的滑跳運(yùn)動(dòng)。 垂向受力峰值由887.2 N降為792.7 N;水平方向受力峰值由-205.5 N 變?yōu)?165.3 N。 通過圖33 和圖34 可知,考慮彈性效應(yīng)之后,球冠體每次觸水滑跳的彈起高度降低,彈起速度也降低。 同時(shí),每次觸水彈起飛行的距離也會(huì)降低,在同樣的時(shí)間或者距離內(nèi)觸水彈起的次數(shù)會(huì)增多。 通過圖35 可以知道,考慮彈性效應(yīng)之后, 第1 次觸水滑跳后水平速度降為23.88 m/s,衰減率為4. 48%。 這說明考慮彈性效應(yīng)后,球冠體每次觸水滑跳水平速度衰減率降低。

通過圖36 可知,球冠體最低點(diǎn)的垂向彈性形變?yōu)橐恢毕蛏蠅嚎s,水平方向彈性形變?yōu)橄蚝笞冃巍?這與高速斜入水沖擊載荷的作用方向一致。此外,垂直方向的彈性形變大于水平方向的彈性形變。

4 結(jié) 論

1) 本文基于邊界元法,建立了一種三維結(jié)構(gòu)體近水面滑跳運(yùn)動(dòng)的時(shí)域數(shù)值仿真方法。 該方法既可以應(yīng)用于任意外形的三維結(jié)構(gòu)體,也擁有較高的計(jì)算效率。 針對(duì)三維剛性球冠體的近水面滑跳運(yùn)動(dòng)進(jìn)行了時(shí)域數(shù)值仿真,得到了對(duì)應(yīng)的運(yùn)動(dòng)特性,驗(yàn)證了該方法的可行性和準(zhǔn)確性。 在此基礎(chǔ)上,可以進(jìn)一步應(yīng)用于擊水式飛行器的滑跳彈道仿真分析。

2) 在球冠體的每次觸水過程中,高速斜入水沖擊載荷提供的阻力和升力,均快速達(dá)到峰值并逐漸減小,且垂向受力相較于水平方向受力會(huì)先一步達(dá)到峰值。 垂向受力的數(shù)值相對(duì)較大,水平方向的受力數(shù)值相對(duì)較小。

3) 初始參數(shù)對(duì)球冠體滑跳運(yùn)動(dòng)有較大的影響。 隨著結(jié)構(gòu)質(zhì)量的增加,每次觸水滑跳的彈起速度和高度增大,每次滑跳落點(diǎn)之間的水平距離加大,在相同的時(shí)間和距離內(nèi)觸水滑跳的次數(shù)減少,每次觸水彈跳水平速度的衰減量增大。 初始高度的增加,會(huì)導(dǎo)致每次觸水滑跳的彈起高度和彈起速度增大,每次滑跳落點(diǎn)之間的水平飛行距離加大,在同樣的時(shí)間和距離內(nèi)觸水滑跳的次數(shù)變少,每次觸水彈跳水平速度的衰減量增大。 隨著水平拋出速度的增大,每次觸水滑跳的彈起速度和高度增大,每次滑跳落點(diǎn)之間的水平距離加大,在相同的時(shí)間和距離內(nèi)觸水滑跳的次數(shù)減少。球冠體半徑的減小,會(huì)導(dǎo)致每次觸水滑跳的彈起高度和彈起速度增加,每次觸水彈跳運(yùn)動(dòng)之間飛行的距離增大,在同樣的時(shí)間或者距離內(nèi)觸水彈起的次數(shù)減少,每次觸水彈跳水平速度的衰減量增大。

4) 彈性效應(yīng)對(duì)球冠體滑跳運(yùn)動(dòng)有一定的影響,考慮彈性效應(yīng)后,會(huì)減小高速斜入水沖擊載荷,導(dǎo)致球冠體彈起速度和高度降低,且滑跳過程中水平速度衰減率減小。 因此,對(duì)彈性效應(yīng)較明顯的結(jié)構(gòu)體進(jìn)行滑跳運(yùn)動(dòng)數(shù)值仿真時(shí),應(yīng)當(dāng)考慮彈性效應(yīng)的影響。