斜浪中船舶運動的三維數(shù)值模擬

斜浪中船舶運動的三維數(shù)值模擬

船舶抗波性是指船舶在一定條件下的運行性能，它關(guān)系到船舶的航行安全和海上武器的使用效率。對船舶耐波性的預(yù)報是船舶水動力領(lǐng)域的一項重要研究內(nèi)容,因此,建立一套有效的船舶耐波性預(yù)報方法具有重要意義。目前工程上一般采用基于勢流理論的切片法對船舶耐波性進(jìn)行預(yù)報,然而第24屆ITTC耐波性委員會最終報告和建議指出:“耐波性數(shù)值計算還很不成熟。目前很多計算工具是基于勢流理論求解的,考慮粘性的船舶耐波性水動力計算的非常少,而且多數(shù)為微幅運動?！眲萘鞣椒ê雎粤肆黧w的粘性作用,且對具有強(qiáng)非線性特征問題的處理能力欠缺,使其應(yīng)用有一定的局限性。近年來,隨著數(shù)值計算技術(shù)的發(fā)展及計算機(jī)硬件水平的提升,計算流體動力學(xué)(computationalfluiddynamics,CFD)方法在船舶水動力研究領(lǐng)域取得了豐碩成果,并逐步應(yīng)用于船舶快速性、操縱性、耐波性等性能的預(yù)報,且具有費用低、無觸點流場測量、可獲得詳細(xì)的流場信息等優(yōu)勢,其應(yīng)用前景相當(dāng)廣闊,是今后發(fā)展的重要方向之一。雖然CFD方法在船舶耐波性研究中已經(jīng)取得了一些進(jìn)展,但從Gothenburg2010CFDWorkshop上所發(fā)表的論文來看,目前的研究多針對船舶頂浪縱搖和垂蕩兩自由度運動。國內(nèi)亦有學(xué)者進(jìn)行了簡單船型的耐波性CFD研究,但總體來說,研究進(jìn)展相對滯后。實際上,船舶在海上多為斜浪航行,且伴有多自由度耦合搖蕩運動,因此,研究實際船型在復(fù)雜環(huán)境條件下的響應(yīng),是未來船舶CFD研究的發(fā)展趨勢。為了擴(kuò)展CFD方法在船舶耐波性研究中的應(yīng)用,本文建立了三維粘性流數(shù)值波浪水池,在此基礎(chǔ)上,對帶球鼻艏和方艉的DTMB5512船模在規(guī)則波中斜浪航行的垂蕩、縱搖及橫搖三自由度耦合運動進(jìn)行了數(shù)值模擬研究。1含實施含流數(shù)值波浪水池控制方程本文的數(shù)值仿真是在計算機(jī)上建立的數(shù)值波浪水池中進(jìn)行的。粘性流數(shù)值波浪水池的控制方程包括連續(xù)性方程、動量方程、湍流模型方程;船舶耐波性數(shù)值模擬必須考慮水和空氣之間的自由面,自由面的模擬采用VOF方法。1.1u3000流體速度及流體密度1、3船舶流場可視為不可壓縮粘性流場,其連續(xù)性方程和動量方程為?ρ?t+?(ρui)?xi=0i=1?2?3(1)?ρ?t+?(ρui)?xi=0i=1?2?3(1)?(ρui)?t+?(ρuiuj)?xj=??xj?(ρui)?t+?(ρuiuj)?xj=??xj(μ?ui?xj-ρˉu′iu′j)-?p?xi+ρfii=1?2?3(2)(μ?ui?xj?ρu′iu′jˉˉˉˉˉˉˉˉ)??p?xi+ρfii=1?2?3(2)其中:ui為流體時均速度分量;p為流體壓強(qiáng);fi為流體體積力分量;ρ為流體密度;u為流體的動力粘性系數(shù);u′i為相對于時均流速的湍流脈動速度分量;-ρˉu′iu′j?ρu′iu′jˉˉˉˉˉˉ為雷諾應(yīng)力。1.2采用船舶計算流體力學(xué)中應(yīng)用較為廣泛的RNGk-ε二方程湍流模型封閉上述方程組。該模型是由重正化群理論推導(dǎo)得出,其湍流動能和湍流耗散率方程為?(ρk)?t+?(ρkui)?xi=??xj?(ρk)?t+?(ρkui)?xi=??xj[αkμeff?k?xj][αkμeff?k?xj]+Gk-ρε(3)?(ρε)?t+?(ρεui)?xi=??xj[αεμeff?ε?xj]+C*1εkGk-C2ερε2k(4)其中:μeff=μ+μt?μt=ρCμk2ε?Cμ=0.0845?αk=αε=1.39?C*1ε=C1ε-η(1-η/η0)1+βη3?C1ε=1.42?C2ε=1.68?η=(2Eij?Eji)1/2kε?Eij=12(?ui?xj+?uj?xi)?η0=4.378?β=0.012。1.3網(wǎng)格單元的體積分?jǐn)?shù)cq自由液面的模擬采用VOF(volumeoffluid)方法,通過計算每個網(wǎng)格單元的體積分?jǐn)?shù)Cq確定自由面位置。體積分?jǐn)?shù)Cq表示單元內(nèi)第q項流體體積與總體積的比例,其輸運方程為{?Cq?t+u?Cq?x+v?Cq?y+w?Cq?z=0n∑q=1Cq=1q=1?2(5)其中:q=1為空氣項,q=2為水項。2波水池建設(shè)2.1規(guī)則波的波高和速度方程采用邊界條件造波法生成波浪,即通過給定造波邊界處波高和波速的解析解,實現(xiàn)數(shù)值造波。若船模前進(jìn)的速度為U,浪向角為χ,則規(guī)則波的波高和速度方程分別為η=Acos(6){u=Aωekzcos[k(xcosχ+ysinχ-Utcosχ-ωt]cosχ+Uv=Aωekzcos[k(xcosχ+ysinχ-Utcosχ-ωt]sinχw=Aωekzsin[k(xcosχ+ysinχ-Utcosχ-ωt](7)其中:A為波幅;k為波數(shù);ω為波浪的圓頻率。2.2阻尼消波作用采用阻尼消波方法消波。阻尼消波方法是指在消波段區(qū)域的動量方程中加上阻尼源項來消波,其對來波的頻率和波長不敏感,可以有效地消除各種頻率和波長的來波。本文在計算域出口邊界前設(shè)置1倍波長的阻尼消波區(qū),同時在消波區(qū)內(nèi)沿波浪傳播方向網(wǎng)格逐漸變疏,這樣可以增大人工阻尼,也起到了消波作用。在阻尼消波區(qū)內(nèi),動量方程為{?u?t+u?u?x+w?u?z=-1ρ?p?x+v(?2u?x2+?2u?z2)-μ(x)u?w?t+u?w?x+w?w?z=g-1ρ?p?z+v(?2w?x2+?2w?z2)-μ(x)w(8)其中:μ(x,y)為在阻尼段起點為零的單調(diào)遞增函數(shù),可以取為線性遞增、指數(shù)遞增等形式,本文取μ(x?y)=105√x2+y2x2max+y2max(9)其中:xmax和ymax分別為消波區(qū)沿x軸和y軸的最大坐標(biāo)值。3船舶運動的數(shù)學(xué)模型3.1根據(jù)接觸問題的金融支付規(guī)則為了計算船舶在波浪中的搖蕩運動,需定義三個坐標(biāo)系統(tǒng):第一個是固定坐標(biāo)系,其空間位置不變,流體的速度、壓力等在該坐標(biāo)系下求解;第二個是平移坐標(biāo)系,其隨船體作平移運動;第三個是隨船坐標(biāo)系,其隨船體作旋轉(zhuǎn)運動,船體的搖蕩運動在該坐標(biāo)系下求解。設(shè)Λ=(φθψ)為船舶旋轉(zhuǎn)的歐拉角,對應(yīng)的分別為船舶的橫搖角、縱搖角和首搖角。如果只考慮角位置關(guān)系,向量在固定坐標(biāo)系與隨船坐標(biāo)系之間的轉(zhuǎn)換矩陣為[˙φ˙η˙ψ]=[1sinφtanθcosφtanθ0cosφ-sinφ0sinφ/cosθcosφ/cosθ]?[pqr](10)設(shè)Ω=(pqr)為角速度在運動坐標(biāo)系下的投影,則角速度在隨船坐標(biāo)系下的投影與歐拉角的導(dǎo)數(shù)之間的轉(zhuǎn)換關(guān)系為S=[cosφcosθcosψsinθsinφ-cosψsinθcosφ+sinψcosφsinψsinφsinψcosθsinψsinθsinφ+sinsψsinθcosφ-cosψcosφcosψsinφ-sinθcosθsinφcosθcosφ](11)船體可視為剛體,船舶搖蕩運動遵循剛體運動的動量定理和動量矩定理。當(dāng)重心位于原點時,在隨船坐標(biāo)系下,船舶空間運動一般方程為{m(˙u-vr+wq)=Xm(˙v-wp+ur)=Ym(˙w-uq+vp)=ΖΙx˙p+(Ιz-Ιy)qr=ΚΙy˙q+(Ιx-Ιz)rp=ΜΙz˙r+(Ιy-Ιx)pq=Ν(12)其中:F=(XYZ)為船體所受外力,L=(KMN)為船體所受外力合力對原點的力矩。3.2數(shù)值模擬的步驟船舶搖蕩運動數(shù)值模擬的關(guān)鍵是如何實現(xiàn)流場與船體運動的耦合求解,本文數(shù)值模擬的步驟如圖1所示。通過上述數(shù)值計算步驟,即可完成船舶粘性流場與船體搖蕩運動的耦合求解。4船舶搖動運動值的模擬4.1自由度耦合運動數(shù)值模擬選用DTMB5512船模作為數(shù)值計算對象,對船模中速(Fr=0.280)、斜頂浪航行(浪舷角χ=45°)時船模三自由度耦合運動進(jìn)行數(shù)值模擬。DTMB5512船模的主尺度數(shù)據(jù)如表1所示,數(shù)值模擬的方案如表2所示。4.2搖蕩運動數(shù)值模擬方法船舶斜浪航行時的計算域設(shè)置如圖2所示,波浪沿x方向和y方向均有速度,需要將計算域的一個側(cè)邊界也設(shè)置為造波邊界,并給定波高和速度。由于波浪亦在y方向傳播,因此需要將消波區(qū)的范圍在y方向上擴(kuò)展;否則波浪傳到側(cè)邊界時,將會產(chǎn)生反射波,影響工作區(qū)流場。圖3為計算域劃分示意圖。船舶搖蕩運動的模擬必須依靠動網(wǎng)格技術(shù)來實現(xiàn)。本文采用網(wǎng)格整體移動方法與滑移網(wǎng)格方法相結(jié)合的動網(wǎng)格方案,對船舶三自由度搖蕩運動進(jìn)行數(shù)值模擬。該方法思路是:根據(jù)船體所受的力和力矩計算出船體的線速度和角速度,使計算域所有網(wǎng)格,隨船體一起進(jìn)行縱搖和垂蕩運動,同時采用滑移網(wǎng)格方法,使船體附近的圓柱體區(qū)域隨船體進(jìn)行橫搖運動,這樣可以保持計算域內(nèi)網(wǎng)格的空間拓?fù)浣Y(jié)構(gòu)不變,在數(shù)值計算過程中,有效避免因網(wǎng)格質(zhì)量下降而導(dǎo)致數(shù)值計算發(fā)散問題。采用有限體積法對微分方程進(jìn)行離散,采用SIMPLE算法對流場進(jìn)行解算。壓力插值采用PRESTO格式,動量、湍流動能和湍流耗散率均采用二階迎風(fēng)格式離散,體積分?jǐn)?shù)采用ModifiedHRIC格式離散。4.3船模斜浪行車的發(fā)生率為了驗證數(shù)值造波的精確性,本文在不帶船模的數(shù)值水池中進(jìn)行了波浪的數(shù)值模擬(去除船行波影響),斜向流場的瞬時波面如圖4所示,在船中位置布設(shè)波高儀記錄波面時歷。圖5為fe=1.35Hz時的測試結(jié)果。從圖中可以看出,CFD仿真得到的波浪時歷與理論值吻合良好,證明了本文方法生成的斜向規(guī)則波具有較高的精度,可為船舶耐波性數(shù)值計算提供高精度的數(shù)值波浪環(huán)境。圖6為船模斜浪航行的瞬時波高圖。圖7～9為在遭遇頻率fe=1.2時,船模的縱搖、垂蕩及橫搖時歷。從圖中可以看出,船模斜浪航行時,縱搖和垂蕩運動穩(wěn)定比較快,而橫搖運動穩(wěn)定相對較慢,在數(shù)個周期后達(dá)到穩(wěn)定。船舶斜浪航行垂蕩、縱搖及橫搖運動頻率響應(yīng)函數(shù)的CFD計算結(jié)果與線性切片理論計算結(jié)果的對比如圖10～12所示。從圖中可以看出,本文計算得到的垂蕩和縱搖運動的頻率響應(yīng)函數(shù)結(jié)果略偏小,而橫搖的頻率響應(yīng)函數(shù)結(jié)果略偏大。造成這種差異的原因可能是由于線性切片理論沒有考慮流體的粘性影響,也沒有考慮縱向運動(縱搖和垂蕩)與橫向運動(橫搖)之間的耦合影響,而CFD方法可以計入上述因素的影響?？傮w看來,本文數(shù)值模擬得到的頻率響應(yīng)函數(shù)與線性勢流理論計算結(jié)果基本符合,證明了本文方法研究船舶斜浪航行的