基于cfd方法的長(zhǎng)峰不規(guī)則波數(shù)值模擬

基于cfd方法的長(zhǎng)峰不規(guī)則波數(shù)值模擬

船舶的垂直移動(dòng)運(yùn)動(dòng)主要包括船舶的縱向移動(dòng)和懸掛。嚴(yán)重的縱搖和垂蕩將引起甲板上浪、砰擊和失速等一系列后果,對(duì)船舶在海上的航行安全以及艦載武器的使用影響極大。因此,本文選取船舶的縱向運(yùn)動(dòng)作為研究對(duì)象,基于CFD方法對(duì)船模在長(zhǎng)峰不規(guī)則波中頂浪運(yùn)動(dòng)進(jìn)行數(shù)值模擬,具有迫切的現(xiàn)實(shí)意義。1浮體運(yùn)動(dòng)數(shù)值模擬設(shè)施要實(shí)現(xiàn)對(duì)船模在長(zhǎng)峰不規(guī)則波中頂浪縱向運(yùn)動(dòng)的數(shù)值模擬,首要前提就是構(gòu)建精度滿足耐波性計(jì)算要求的長(zhǎng)峰不規(guī)則波數(shù)值波浪水池。所謂的數(shù)值波浪水池,是指對(duì)非線性波浪水動(dòng)力以及浮體運(yùn)動(dòng)數(shù)值模擬設(shè)施的統(tǒng)稱。它能夠通過(guò)實(shí)驗(yàn)觀測(cè)為各種波浪理論的研究奠定堅(jiān)實(shí)的物理基礎(chǔ),還能為海洋、船舶等工程設(shè)計(jì)提供可靠而高效的試驗(yàn)數(shù)據(jù)。本文基于粘性流理論構(gòu)建長(zhǎng)峰波數(shù)值波浪水池,采用有限體積法(FVM)對(duì)RANS方程和連續(xù)控制方程進(jìn)行離散求解,利用Fluent軟件的二次開發(fā)功能UDF,編寫邊界條件完成數(shù)值造波和阻尼消波功能。選用ITTC海浪譜作為目標(biāo)海浪譜,構(gòu)建長(zhǎng)峰波數(shù)值波浪水池,并對(duì)其計(jì)算精度進(jìn)行誤差計(jì)算。1.1波水池價(jià)值1.1.1雙轉(zhuǎn)臺(tái)波的入射速度數(shù)值造波是指用數(shù)值方法模擬波浪的生成過(guò)程,為數(shù)值模型實(shí)驗(yàn)提供各種形式的波浪環(huán)境條件。本文采用的數(shù)值造波方法是邊界條件造波法,即在數(shù)值波浪水池入口邊界設(shè)置3個(gè)方向的速度。X方向沿水池向下游為正,Z方向向上為正,Y軸與X軸和Z軸符合右手法則。3個(gè)方向的入射速度分別為:???????????????????U=∑i=1NωiAiekizcos(kix?ωit+εi),V=0,W=∑i=1NωiAiekizsin(kix?ωit+εi)。(1)[BFQ]{U=∑i=1ΝωiAiekizcos(kix-ωit+εi),V=0,W=∑i=1ΝωiAiekizsin(kix-ωit+εi)。(1)[BFQ]式中:η為波動(dòng)水面相對(duì)于靜止水面的瞬時(shí)高度;Ai,ki,ωi和εi分別為第i個(gè)組成波的波幅、波數(shù)、圓頻率和初始相位,εi是在(0,2π)范圍內(nèi)的隨機(jī)相位;X軸為波浪傳播方向;U,V,W分別為波浪X軸、Y軸和Z軸的速度分量。1.1.2阻尼消波器計(jì)算當(dāng)波浪到達(dá)水池末端開邊界處,會(huì)引發(fā)水波的二次反射,反射波與入射波疊加,會(huì)導(dǎo)致波場(chǎng)的失真。因此,消除反射波的影響也是數(shù)值波浪水池的一項(xiàng)重要技術(shù)。數(shù)值波浪水池常用的消波技術(shù)主要有輻射邊界條件消波、主動(dòng)消波器消波、阻尼消波3種。本文選取的是阻尼消波。阻尼消波法是指在流場(chǎng)中添加人工粘性,因其對(duì)來(lái)波的頻率和波長(zhǎng)不敏感,可以有效地消除各種頻率和波長(zhǎng)的來(lái)波,因而被廣泛采用。在計(jì)算域出口邊界前設(shè)置1～2倍波長(zhǎng)的阻尼消波段,利用Fluent中的UDF宏DEFINE_SOURCE(mom_source,cell,thread,dS,eqn)編程實(shí)現(xiàn)消波。在阻尼消波段內(nèi),動(dòng)量方程寫為:?u?t+u?u?x+v?u?y=?1ρ?p?x+v(?2u?x2+?2u?y2)?μ(x)u,(2)?u?t+u?u?x+v?u?y=-1ρ?p?x+v(?2u?x2+?2u?y2)-μ(x)u,(2)?v?t+u?v?x+v?v?y=g?1ρ?p?y+v(?2v?x2+?2v?y2)?μ(x)v?v?t+u?v?x+v?v?y=g-1ρ?p?y+v(?2v?x2+?2v?y2)-μ(x)v。(3)其中μ(x)為在阻尼段起點(diǎn)為0的單調(diào)遞增函數(shù),可以取為線性遞增、指數(shù)遞增等形式。取μ=105(x-Xmin_D)/(Xmax_D-Xmin_D)。(4)式中:Xmin_D和Xmax_D分別為消波區(qū)的最小、最大X坐標(biāo)。1.2長(zhǎng)峰不規(guī)則波值波池建設(shè)1.2.1垂直網(wǎng)格化社區(qū)的網(wǎng)格結(jié)構(gòu)三維數(shù)值水槽的網(wǎng)格劃分如圖1所示,本文構(gòu)建的長(zhǎng)峰不規(guī)則波數(shù)值水池長(zhǎng)18m,其中12～18m為消波區(qū),寬3m,深2.5m,自由面以上1.5m,整個(gè)水槽的網(wǎng)格數(shù)為413820。垂直自由面方向的網(wǎng)格尺寸取為有義波高1/5。造波區(qū)沿X軸正方向網(wǎng)格尺寸與垂直于自由面Z方向的最小網(wǎng)格尺寸相同,從自由面到水槽頂部網(wǎng)格按1∶1.1的比例等比分布;從自由面到水池底部網(wǎng)格按1∶1.05的比列等比分布。網(wǎng)格基本上是離自由面越遠(yuǎn)尺寸越大。對(duì)于消波區(qū)的網(wǎng)格劃分,垂直方向劃分與造波區(qū)的一致,水平方向網(wǎng)格以造波區(qū)網(wǎng)格尺度為基準(zhǔn)向右邊界逐漸擴(kuò)大。1.2.2計(jì)算示例描述選取ITTC海浪譜作為目標(biāo)靶譜,對(duì)3種海況下的長(zhǎng)峰不規(guī)則波進(jìn)行數(shù)值模擬。波浪的目標(biāo)參數(shù)見(jiàn)表1。1.2.3長(zhǎng)峰不法波數(shù)值波浪水池波面時(shí)歷曲線長(zhǎng)峰不規(guī)則波數(shù)值模擬瞬時(shí)波面場(chǎng)和局部速度矢量如圖2和圖3所示。圖4和圖5是有義波高分別取H1/3=0.08m,H1/3=0.125m時(shí),長(zhǎng)峰不規(guī)則波數(shù)值波浪水池X=1m,X=6m,X=11m處的波面時(shí)歷曲線對(duì)比圖。從上述數(shù)值造波水池不同位置的波面時(shí)歷曲線對(duì)比可看出,長(zhǎng)峰不規(guī)則波在沿X軸正方向向下游傳播過(guò)程中,波能有一定程度的衰減,這是由于波浪在傳播過(guò)程中,高頻子波衰減導(dǎo)致的。1.2.4目標(biāo)譜的建立對(duì)上面監(jiān)測(cè)得到長(zhǎng)峰不規(guī)則波波面時(shí)歷進(jìn)行譜分析與目標(biāo)譜對(duì)比如圖6所示。分別從譜面積m0、譜峰頻率ωp和有義波高H1/33個(gè)方面對(duì)上述數(shù)值模擬海浪譜進(jìn)行誤差分析,誤差計(jì)算結(jié)果如表2和表3所示。2船舶模型中的行程曲線模型對(duì)長(zhǎng)峰不規(guī)則波中的波浪垂直運(yùn)動(dòng)進(jìn)行了數(shù)值模擬2.1海上模型幾何模型選取具有球鼻首和方位的DTMB5512船模作為研究對(duì)象,該船模是ITTC(國(guó)際船模試驗(yàn)水池會(huì)議)推薦的瘦削型標(biāo)準(zhǔn)船模DTMB5415的全相似幾何模型,以美國(guó)海軍DDG-51型驅(qū)逐艦為模板,船型數(shù)據(jù)詳實(shí),而且與我海軍艦船船型相似(見(jiàn)圖7),適合軍艦參考。表4為該船模及其所對(duì)應(yīng)的實(shí)船尺度的主要數(shù)據(jù)。2.2數(shù)值模擬方案2.2.1垂蕩和縱搖2個(gè)自由度的運(yùn)動(dòng)本文對(duì)DTMB5512船模在長(zhǎng)峰不規(guī)則波中頂浪縱向運(yùn)動(dòng)數(shù)值模擬試驗(yàn)中,只考慮垂蕩和縱搖2個(gè)自由度的運(yùn)動(dòng)。船模CFD耐波性數(shù)值模擬試驗(yàn)對(duì)計(jì)算資源的要求較高,考慮到本文研究所使用的計(jì)算機(jī)配置的實(shí)際情況,以及試驗(yàn)水池網(wǎng)格劃分所帶來(lái)的計(jì)算效率等問(wèn)題,試驗(yàn)計(jì)劃如表5所示。2.2.2耐波性數(shù)值波浪水池網(wǎng)格劃分參照《水面船模耐波性實(shí)驗(yàn)規(guī)程》,將計(jì)算域設(shè)置成長(zhǎng)方體形狀,如圖8所示。船模與計(jì)算域各邊界的位置關(guān)系如下:入口距船首1倍船長(zhǎng),出口距船尾2倍船長(zhǎng),頂部邊界距水線0.5倍船長(zhǎng),底部邊界距水線1倍船長(zhǎng),左、右邊界距船中縱剖面0.5倍船長(zhǎng)。耐波性數(shù)值波浪水池分成5個(gè)區(qū)域進(jìn)行網(wǎng)格劃分,即近船體區(qū)域、近流場(chǎng)區(qū)域、自由面區(qū)域、上下遠(yuǎn)流場(chǎng)區(qū)域及消波區(qū)域、各區(qū)之間互不重疊,且連接處選擇connected方式,如圖9所示。近船體區(qū)網(wǎng)格劃分如圖10所示。2.3udf模擬水池試驗(yàn)過(guò)程根據(jù)船舶六自由度運(yùn)動(dòng)的控制方程(5),當(dāng)波浪作用于船體時(shí),其運(yùn)動(dòng)的速度、角速度以及位置、姿態(tài)等可以通過(guò)控制方程求解、積分得到。對(duì)于流浮耦合運(yùn)動(dòng),波浪作用在船體上的力和力矩使船體產(chǎn)生運(yùn)動(dòng),同時(shí)船體的運(yùn)動(dòng)又對(duì)其周圍流場(chǎng)產(chǎn)生影響。因此本文在數(shù)值模擬中,分段計(jì)算流體與船體運(yùn)動(dòng)的耦合,步驟如下:1)將船模按初始浮態(tài)固定,原點(diǎn)與重心重合;2)對(duì)流場(chǎng)進(jìn)行初始化,設(shè)定初始航速后造波;3)以時(shí)間步長(zhǎng)Δt=0.001s步進(jìn);4)通過(guò)當(dāng)前流場(chǎng)變量迭代求解流場(chǎng)的速度矢量;5)通過(guò)壓力—速度耦合算法獲得壓力場(chǎng);6)求解體積分?jǐn)?shù)方程重構(gòu)自由面;7)讀取船體受力F和力矩M,根據(jù)控制方程mz??=F,I?????=Mmz??=F,Ι?????=Μ求解船體運(yùn)動(dòng)要素,經(jīng)積分得到船體x方向的阻力、垂向受到的垂向力合縱搖角速度及縱搖幅值;8)更新船體位置和浮態(tài);9)返回第3步,求解改變浮態(tài)后各量,并根據(jù)計(jì)算再次改變浮態(tài),按此迭代求解,直到方程組的殘差小于設(shè)定值或迭代次數(shù)達(dá)到設(shè)定值;10)返回第2步并重復(fù)以下步驟,直到設(shè)定的時(shí)間步數(shù)計(jì)算完畢。mz??=?z??∫Lm′dx?z˙∫L(N′?vdm′dx)dx?2ρgz∫Lywdx+???∫Lm′xdx+?˙∫L(N′x?2m′v?xvdm′dx)dx+?∫L(2ρgywx?N′v+v2dm′dx)dx+ζA[∫L2ρgywe?kTcossinkxdxmω∫L(N′?vdm′dx)ekTsincoskxdx?ω2∫Lm′ekTcossinkxdx]cossinωet。(5)[BFQ]mz??=-z??∫Lm′dx-z˙∫L(Ν′-vdm′dx)dx-2ρgz∫Lywdx+???∫Lm′xdx+?˙∫L(Ν′x-2m′v-xvdm′dx)dx+?∫L(2ρgywx-Ν′v+v2dm′dx)dx+ζA[∫L2ρgywe-kΤsincoskxdxmω∫L(Ν′-vdm′dx)ekΤcossinkxdx-ω2∫Lm′ekΤsincoskxdx]sincosωet。(5)[BFQ]通過(guò)以上迭代、循環(huán),可實(shí)現(xiàn)流體與船體運(yùn)動(dòng)的耦合。步驟2按照1.2節(jié)構(gòu)建的數(shù)值波浪水池進(jìn)行造波和消波。試驗(yàn)過(guò)程中可以將事先保存好的穩(wěn)定流場(chǎng)導(dǎo)入耐波性數(shù)值模擬水池,進(jìn)行數(shù)值模擬,這樣可以提高計(jì)算效率,縮短試驗(yàn)時(shí)間。步驟4～6,按前文設(shè)置的數(shù)值方法進(jìn)行計(jì)算。步驟7則通過(guò)UDF編程實(shí)現(xiàn)。步驟9殘差標(biāo)準(zhǔn)可取軟件默認(rèn)設(shè)定值,迭代次數(shù)上限為40步。2.4縱向運(yùn)動(dòng)模擬系統(tǒng)主要解決了以下方面的問(wèn)題對(duì)于DTMB5512船模由于帶球鼻首、具有方尾,線性復(fù)雜,所以實(shí)現(xiàn)其耐波性的數(shù)值模擬,要比一般商用船模困難得多?？刹捎靡韵路椒▽?duì)實(shí)驗(yàn)進(jìn)行改進(jìn):1)改進(jìn)船模貼體網(wǎng)格質(zhì)量、數(shù)量及分區(qū)網(wǎng)格的匹配;2)采用遞增方法實(shí)現(xiàn)船模的縱向運(yùn)動(dòng),即先對(duì)船模進(jìn)行單自由度縱搖,穩(wěn)定之后再增加垂蕩的數(shù)值模擬;3)逐漸增加船模質(zhì)量(將船模質(zhì)量降低到原船模的一半進(jìn)行數(shù)值模擬,當(dāng)殘差穩(wěn)定后再逐漸增大質(zhì)量直至原值);4)減小欠松弛因子,控制單元網(wǎng)格內(nèi)速度的變化量。通過(guò)上述方法,解決了DTMB5512船模復(fù)雜的幾何船形與波浪作用的數(shù)值問(wèn)題。DTMB5512型船模在長(zhǎng)峰不規(guī)則波中縱向運(yùn)動(dòng)數(shù)值模擬的壓力場(chǎng)和速度場(chǎng),如圖11和圖12所示。DTMB5512型船模在長(zhǎng)峰不規(guī)則波中縱搖及垂蕩的時(shí)歷曲線如圖13和圖14所示。3頂浪運(yùn)動(dòng)時(shí)歷曲線在計(jì)算精度驗(yàn)證中的應(yīng)用由于本文對(duì)船模在長(zhǎng)峰不規(guī)則波中頂浪縱向運(yùn)動(dòng)的數(shù)值模擬研究目前在國(guó)內(nèi)外還處于起步階段,未找到具體的水池實(shí)驗(yàn)數(shù)據(jù)?？紤]到基于勢(shì)流理論艦船六自由度計(jì)算,在理論上是成熟的,并在工程應(yīng)用上取得了很多成果,得到了水動(dòng)力學(xué)界認(rèn)可。一般來(lái)說(shuō),SCFD方法由于考慮到流體粘性,計(jì)算精度應(yīng)略高于勢(shì)流理論計(jì)算結(jié)果,但沒(méi)有本質(zhì)上的差異。至于非線性搖蕩則需要另作考慮,本課題局限于線性搖蕩,所以用勢(shì)流理論計(jì)算結(jié)果進(jìn)行驗(yàn)證是可行的。由于勢(shì)流理論切片法計(jì)算過(guò)程相對(duì)繁瑣,可參考文獻(xiàn),本文直接給出DTMB5512型船模在算例波浪環(huán)境下縱搖和垂蕩的響應(yīng)方差。將數(shù)值模擬的DTMB5512型船模搖蕩運(yùn)動(dòng)時(shí)歷曲線運(yùn)用線性譜分析方法進(jìn)行分析,獲得縱搖及垂蕩的搖蕩譜,如圖15和圖16所示。對(duì)上述搖蕩譜進(jìn)行積分即可得到本文數(shù)值模擬DTMB5512船模頂浪運(yùn)動(dòng)縱搖和垂蕩的響應(yīng)方差。對(duì)DTMB5512型船模在長(zhǎng)峰不規(guī)則波中頂浪縱向運(yùn)動(dòng)的數(shù)值結(jié)果與勢(shì)流理論計(jì)算結(jié)果進(jìn)行對(duì)比得出相對(duì)誤差如表6所示。4船模、船模間的布置要求本文運(yùn)用SCFD方法對(duì)艦船在長(zhǎng)峰不規(guī)則波中頂浪運(yùn)動(dòng)進(jìn)行數(shù)值模擬,得到如下可供參考的經(jīng)驗(yàn):1)船舶搖蕩、阻力、操縱、推進(jìn)等課題的數(shù)值研究,技術(shù)細(xì)節(jié)上的一個(gè)主要不同,體現(xiàn)為網(wǎng)格的布設(shè)上。2)船舶耐波性研究的網(wǎng)格,必須將造波和船體網(wǎng)格、動(dòng)網(wǎng)格三者進(jìn)行很好的協(xié)調(diào)、匹配,做到三者的有機(jī)結(jié)合,否則會(huì)出現(xiàn)計(jì)算發(fā)散及非物理現(xiàn)象等不合理現(xiàn)象的發(fā)生。3)船舶多自由度搖蕩試驗(yàn)?zāi)壳暗睦щy主要集中在計(jì)算資源上,由于其計(jì)算量巨大,PC機(jī)以及低端的工作站、服務(wù)器已經(jīng)不能滿足其正常情況下的計(jì)算需要。如果計(jì)算資源等硬件設(shè)施有限,則需要在離散方法、格式,控制方程、調(diào)節(jié)參數(shù)、UDF開發(fā)等“軟件”上下功夫。4)在三維空間內(nèi)完成船模非規(guī)則波中的搖蕩試驗(yàn),其流場(chǎng)及自由面要比船模在規(guī)則波中的復(fù)雜得多,因此在船模貼體網(wǎng)格的布設(shè)上,要求網(wǎng)格質(zhì)量非常高,同時(shí)在interface交界面處,左右網(wǎng)格要尺寸一致,且在船模外表面盡可能多地使用結(jié)構(gòu)性網(wǎng)格,特別是在阻力計(jì)算上,還要盡可能多地布設(shè)邊界層,以提高計(jì)算的精度。5)耐波性流場(chǎng)的高度復(fù)雜性還表現(xiàn)在輸入、輸出及船模的響應(yīng)及其變化率上。具體表現(xiàn)為1個(gè)網(wǎng)格上,至少要輸入3個(gè)方向的線速度u,v,ω,壓力P,參數(shù)k,ε,流體體積分?jǐn)?shù)