趙九龍,閆發(fā)鎖,樊 磊
(1.哈爾濱工程大學(xué),哈爾濱150001;2.交通運(yùn)輸部 上海打撈局,上海 200090;3.德州農(nóng)工大學(xué),德克薩斯州 77840,美國(guó))
基于改進(jìn)Wagner方法的軸對(duì)稱(chēng)體入水“后時(shí)段”砰擊載荷數(shù)值計(jì)算
趙九龍1,2,閆發(fā)鎖1,樊 磊3
(1.哈爾濱工程大學(xué),哈爾濱150001;2.交通運(yùn)輸部 上海打撈局,上海 200090;3.德州農(nóng)工大學(xué),德克薩斯州 77840,美國(guó))
文章基于三維邊界元方法研究了三維軸對(duì)稱(chēng)體入水砰擊載荷的數(shù)值算法。算法從三維力學(xué)模型出發(fā),繼承Wagner自由液面抬升理論,引入浸深因子Cw以確定自由液面抬升高度,將自由液面線(xiàn)性化處理,同時(shí)考慮網(wǎng)格運(yùn)動(dòng),在自由液面附近對(duì)網(wǎng)格進(jìn)行截?cái)嘀貥?gòu),以確保水下濕面積的精準(zhǔn)。算法中使用考慮加權(quán)運(yùn)動(dòng)項(xiàng)的非線(xiàn)性伯努利方程計(jì)算得到入水結(jié)構(gòu)的表面壓力,進(jìn)一步積分可得入水結(jié)構(gòu)的總體受力;另外,算法中引入虛擬面元的概念,將砰擊載荷計(jì)算時(shí)歷延長(zhǎng)至液面高于墜物制高點(diǎn)之后,擴(kuò)大了傳統(tǒng)入水載荷計(jì)算的時(shí)歷范圍,并且文中引用Alaoui半球入水試驗(yàn),對(duì)算法的正確性及適用性進(jìn)行了驗(yàn)證。
軸對(duì)稱(chēng)體;砰擊載荷;液面升高;網(wǎng)格截?cái)?;虛擬面元
結(jié)構(gòu)物在短時(shí)間內(nèi)迅速進(jìn)入液體,隨之產(chǎn)生液體飛濺的現(xiàn)象稱(chēng)之為砰擊。然而,砰擊涵蓋的面非常廣泛,廣義的砰擊概念是指結(jié)構(gòu)物與流體發(fā)生的碰撞現(xiàn)象,一般都涉及到固、液、氣三種介質(zhì)的相互耦合,是一種瞬態(tài)的強(qiáng)非線(xiàn)性、強(qiáng)耦合和極其復(fù)雜的物理現(xiàn)象。在工程實(shí)際中,砰擊問(wèn)題非常易見(jiàn),例如:水上飛機(jī)水面降落、船舶在惡劣海況中的底部砰擊、艦船上的救生艇落水、航天運(yùn)載器的水面回收、魚(yú)雷空投入水等問(wèn)題都屬于該范疇。
砰擊現(xiàn)象對(duì)應(yīng)的的實(shí)際力學(xué)模型要涉及到固、液、氣三種介質(zhì)的耦合,而且具體情況復(fù)雜多變,目前為止仍未能提出一種貼切真實(shí)情況、普適性強(qiáng)的理論模型。但是為了便于對(duì)問(wèn)題展開(kāi)研究,通常使用“結(jié)構(gòu)物入水沖擊”模型作為具體的力學(xué)模型來(lái)分析。
許多工程結(jié)構(gòu)物在高速入水時(shí),所面臨的挑戰(zhàn)是入水過(guò)程瞬時(shí)所受到的巨大水動(dòng)壓力,這往往會(huì)導(dǎo)致入水結(jié)構(gòu)的局部大變形甚至破斷。近年來(lái)人們?cè)谟?jì)算入水壓力方面做了許多工作,但是大部分都是以二維視角展開(kāi)研究,研究基于的入水模型也多見(jiàn)于二維楔形體,數(shù)值方面做得比較好的有Wu等[1-2]研究出了二維V形楔的迭代求解方法,并且運(yùn)用二維邊界元法,使用柯西積分計(jì)算了楔形體入水的相似解,許國(guó)冬[3]討論了V形楔多種入水形式的相似解。但是由于實(shí)際工程應(yīng)用中人們使用的結(jié)構(gòu)形式多樣,二維入水算法存在著一定的應(yīng)用局限性,因此對(duì)于三維入水算法的開(kāi)發(fā)有著迫切的需求。
目前還沒(méi)有非常成熟的三維入水算法供于應(yīng)用。近年來(lái)的Korobkin和Scolan[4-5]從逆推的Wagner問(wèn)題及線(xiàn)性Wagner問(wèn)題兩方面出發(fā)研究了三維入水砰擊理論,解決了流場(chǎng)流動(dòng)、壓力分布及流域形狀三方面的問(wèn)題。Faltisen和Chezhian[6]以非軸對(duì)稱(chēng)的船體模型入水為研究對(duì)象,基于邊界元方法、發(fā)展的Wagner理論,研究了三維入水砰擊的數(shù)值方法。
基于此,本文也在三維入水方面進(jìn)行了研究,主要研究了基于三維邊界元方法的砰擊理論算法及其在砰擊載荷計(jì)算中的應(yīng)用。數(shù)值方面的研究從三維力學(xué)模型出發(fā),基于三維邊界元方法,繼承Wagner自由液面抬升理論,引入了“浸深因子”Cw[7]以確定自由液面抬升的高度,在這一高度處建立新的自由液面并稱(chēng)之為“等效自由液面”。算法中考慮到網(wǎng)格運(yùn)動(dòng),對(duì)入水過(guò)程中每一時(shí)間步時(shí)刻對(duì)應(yīng)的網(wǎng)格編號(hào)、節(jié)點(diǎn)坐標(biāo)等進(jìn)行追蹤計(jì)算,由此得出網(wǎng)格與等效自由液面的位置關(guān)系,進(jìn)一步可精確追蹤到每一時(shí)間步時(shí)刻對(duì)應(yīng)的濕網(wǎng)格數(shù)據(jù),同時(shí)對(duì)等效自由液面附近的網(wǎng)格根據(jù)實(shí)際情況進(jìn)行截?cái)嘀厣?,以確保水下濕面積的精確;通過(guò)時(shí)間步的循環(huán)計(jì)算,最終可得到入水結(jié)構(gòu)的表面“壓力”和總體“受力”兩種載荷時(shí)歷。另外,提出了虛擬面元的概念,使算法可以應(yīng)用于墜物浸水之后的砰擊時(shí)段(液面高于入水物體的最高點(diǎn),即所謂的入水“后時(shí)段”,將傳統(tǒng)的砰擊載荷計(jì)算時(shí)歷延長(zhǎng)。
1.1 改進(jìn)Wagner方法在砰擊載荷計(jì)算中的應(yīng)用原理
Wagner在Von-Karman[7]線(xiàn)性砰擊理論的基礎(chǔ)上,提出了自由液面的抬升效應(yīng),采用平板擬合法,相對(duì)精確了浸濕半寬、附加質(zhì)量、沖擊力及壓力分布的求解。在求解過(guò)程中所采用的相當(dāng)平板是對(duì)具體物面的近似,而且采用線(xiàn)性伯努利方程對(duì)壓力分布進(jìn)行求解,整個(gè)過(guò)程雖然計(jì)算量小,但是計(jì)算精度還有待進(jìn)一步提高。本文在原始Wagner方法的基礎(chǔ)上進(jìn)行了改進(jìn),引入了邊界元方法[8],考慮具體浸水物面,采用非線(xiàn)性伯努利方程對(duì)壓力進(jìn)行求解。
談及流體載荷的計(jì)算,近年來(lái)出現(xiàn)了多種計(jì)算方法,比較常用的有攝動(dòng)法、有限元法、有限差分法、邊界元法。由于邊界元方法將空間的問(wèn)題轉(zhuǎn)化為邊界表面的積分問(wèn)題,將問(wèn)題減小一維,大大減少了未知量的數(shù)目,另一方面它能夠處理無(wú)窮遠(yuǎn)處的邊界條件,對(duì)外部問(wèn)題有很好的適用性,正是由于這兩點(diǎn)才使得邊界元方法在流體載荷計(jì)算中得以迅速發(fā)展。本文所使用的邊界元方法適用于求解無(wú)界流場(chǎng)中的三維無(wú)升力繞流問(wèn)題,分布奇點(diǎn)法的所有要點(diǎn)都包含在內(nèi)。速度勢(shì)可以用分布源和分布偶極來(lái)描述,本文算法僅僅涉及到分布源的概念。
控制方程及邊界條件如公式(1)、(2):
將濕物面劃分為N個(gè)面元可以得到離散化的方程組:
式中:aij是j面元對(duì)i面元的影響系數(shù):
通過(guò)求解方程組(3)便可以求解出各面元中心點(diǎn)處的分布源密度,進(jìn)一步可得到個(gè)面元中心點(diǎn)的速度勢(shì):
緊接著可以求出各面元中心點(diǎn)處的誘導(dǎo)速度:
再進(jìn)一步可根據(jù)各時(shí)刻各點(diǎn)的速度勢(shì)數(shù)值,對(duì)時(shí)間進(jìn)行微分計(jì)算便可得到這樣就可以使用非線(xiàn)性的伯努利方程:
計(jì)算得到各時(shí)刻對(duì)應(yīng)的繞流壓力。
從上面的推導(dǎo)可以計(jì)算出無(wú)界流中運(yùn)動(dòng)物體的表面繞流壓力值,但是人們通常研究的入水砰擊現(xiàn)象發(fā)生在半無(wú)界流場(chǎng)中,而且正是自由液面的附近。為了解決理論的適用問(wèn)題,考慮線(xiàn)性的自由液面,由于自由液面上的速度勢(shì)φ=0,即可將自由液面下的濕表面沿著自由液面進(jìn)行上下對(duì)稱(chēng)處理,形成一個(gè)沿自由液面對(duì)稱(chēng)的閉合曲面s+s′,這樣就將問(wèn)題等同于考慮在無(wú)界流場(chǎng)中運(yùn)動(dòng)的s+s′,使問(wèn)題適用于上述理論。
1.2 自由液面處理方法
本文算法在對(duì)自由液面進(jìn)行線(xiàn)性化處理的時(shí)候,根據(jù)Wagner理論,考慮到了液面抬升,如圖1所示的入水模型。通常情況下,物體入水過(guò)程中,自由液面與物體相交點(diǎn)處的位置均要高于原始的自由液面,流體沿物面有向上爬升的趨勢(shì),嚴(yán)格來(lái)講計(jì)算過(guò)程中必須要追蹤自由液面的運(yùn)動(dòng)形式,這樣一定會(huì)使求解更加精確,但是對(duì)三維入水模型來(lái)講,自由液面的追蹤非常的困難,而且鮮有適用性比較廣泛的方法,例如Faltinsen和Chezhian[6],在三維船體入水砰擊的自由液面追蹤中也不能夠很精確地描述自由液面的具體形狀及變化,他們所采用的方法是沿著船體與自由液面交界的一圈均勻標(biāo)記若干點(diǎn),之后分別以這些標(biāo)記點(diǎn)為起點(diǎn),平行于原始的自由液面畫(huà)出射線(xiàn),最終由做出的許多射線(xiàn)近似地?cái)M合出一個(gè)抬升后的自由液面。該方法目前來(lái)講無(wú)疑是一種比較準(zhǔn)確的追蹤方法,基于此,在自由液面的處理過(guò)程中,本文程序也沿用了這種“射線(xiàn)延伸”的近似,但是較之有所簡(jiǎn)化。
從工程實(shí)用角度,為了簡(jiǎn)化流固交界附近數(shù)值處理的難度和計(jì)算量。本文算法采用線(xiàn)性化的等效自由液面,使用經(jīng)驗(yàn)公式換算出物面與自由液面每一時(shí)刻的觸點(diǎn)高度,以確定該觸點(diǎn)位置,并以之為起點(diǎn)平行于水平面引出射線(xiàn),認(rèn)為該射線(xiàn)所在的水平面便是新的自由液面,在這里稱(chēng)之為“等效自由液面”。為了確定等效自由液面的位置,同時(shí)引入了浸深因子Cw的概念,如圖2所示。
圖1 Faltinsen關(guān)于自由液面的追蹤Fig.1 The free surface tracking method by Faltinsen
圖2 本文采用的自由液面處理方法Fig.2 The free surface caculation by this paper
如果知道Cw值,便很容易可以確定等效自由液面的位置。本文參照了美國(guó)水面武器研究中心白橡樹(shù)實(shí)驗(yàn)室的一份研究報(bào)告[9],其中詳細(xì)介紹了軸對(duì)稱(chēng)物體垂向入水與任意物體傾斜入水兩種類(lèi)別的Cw值選取辦法。
表1 Cw選取方法Tab.1 Calculation of Cw
顯然使用該方法可以使砰擊載荷計(jì)算得到很大程度上簡(jiǎn)化,提高了計(jì)算效率。
1.3 面元網(wǎng)格的運(yùn)動(dòng)及重構(gòu)
本文將邊界元模型入水運(yùn)動(dòng)的過(guò)程劃分為若干時(shí)間步進(jìn)行分析,對(duì)于每一時(shí)間步時(shí)刻都要統(tǒng)計(jì)出等效自由液面以下的面元及其節(jié)點(diǎn)信息,這是計(jì)算誘導(dǎo)速度系數(shù)的基本前提。很顯然,水面以下遠(yuǎn)離等效自由液面的網(wǎng)格保持原有的完整性,因此不必對(duì)其幾何形狀進(jìn)行改變;而在等效自由液面附近,沿著物體與自由液面的交界線(xiàn)有一圈網(wǎng)格被等效自由液面截?cái)啵尚碌乃倪呅巍⑽暹呅魏腿切蚊嬖?;但是考慮到本文算法中要求提供的面元必須是長(zhǎng)寬比比值適中的四邊形面元,而且為了提高程序的精確性,不能單純地將這部分網(wǎng)格進(jìn)行忽略,因此在程序中對(duì)這部分網(wǎng)格進(jìn)行了截?cái)嘀厣幚?,具體的處理方法如下圖3所示。
另外,考慮到網(wǎng)格運(yùn)動(dòng),于是便不能直接使用常規(guī)的非線(xiàn)性伯努利方程來(lái)計(jì)算載荷,需要對(duì)其進(jìn)行改進(jìn),添加運(yùn)動(dòng)項(xiàng)目。根據(jù)泰勒公式變換得到:
圖3 自由液面交界處的網(wǎng)格截?cái)嗉爸貥?gòu)方法Fig.3 The gridding truncation method nearby the boundary of the free surface and the object
圖4 墜物砰擊后時(shí)段邊界元模型示意圖Fig.4 The BEM of the later entry time
1.4 砰擊后時(shí)段算法適用性研究
按照上述的方法僅僅可以計(jì)算墜物未完全浸水的砰擊時(shí)段,而對(duì)于墜物制高點(diǎn)低于等效液面之后的時(shí)段則無(wú)法進(jìn)行計(jì)算。在物體完全浸水之后,短時(shí)間內(nèi)物體的上表面事實(shí)上還沒(méi)有被液體覆蓋,因此在計(jì)算中物體的濕表面s與投影面s′將無(wú)法組成封閉面,本文所引用的無(wú)界流面元法則不能適用,在后續(xù)求解方程組的過(guò)程中會(huì)出現(xiàn)無(wú)解的情況。
如圖4所示,T1時(shí)刻之后s與s′分離,之后便是本文所指的“后砰擊”時(shí)段。
對(duì)此,本文提出了添加虛擬面元的方法,使入水“后時(shí)段”的砰擊載荷計(jì)算得以實(shí)現(xiàn),下圖5是本文所設(shè)計(jì)的三種入水邊界元模型。圖中a是實(shí)際的入水邊界元模型,b是添加的虛擬邊界元模型,每個(gè)面都有兩種顏色,其中棕色代表物體的觸水濕表面。第一種模型根據(jù)實(shí)際情況對(duì)a圖邊界元模型加了一個(gè)上蓋,并且劃分邊界元網(wǎng)格,即b圖所示圓形結(jié)構(gòu)。這樣在T1時(shí)刻之后,a、b共同組成了一個(gè)閉合的邊界元模型,經(jīng)與抬升后的自由液面對(duì)稱(chēng)便可以形成兩個(gè)閉合的邊界元模型,符合理論要求。如果使用這樣的水動(dòng)力模型進(jìn)行計(jì)算的話(huà),實(shí)際上就是認(rèn)為在半球入水之后液體立即將半球體吞沒(méi),即半球上蓋之上被流體覆蓋。第二種模型,實(shí)際上就是模擬了在T1時(shí)刻之后a的上緣與抬升后的自由液面間的液面形狀,這一段液面上一定存在著分布源,它對(duì)流場(chǎng)性質(zhì)有重要的影響,本文計(jì)算時(shí)將其視作濕表面的一部分,可以計(jì)算出其上的分布源強(qiáng)度,從而對(duì)a表面面元速度勢(shì)進(jìn)行修正,最后在計(jì)算試件垂向壓力時(shí)并不是對(duì)a和b上的所有面元進(jìn)行加權(quán)計(jì)算,僅僅考慮a上的面元壓力,b僅僅是作為一個(gè)虛擬的濕表面。第三種模型類(lèi)似于第二種,只是假想T1時(shí)刻后a上緣的自由液面成圓筒狀,因此將b設(shè)計(jì)成了圓筒狀的一層面元。
圖5 3種虛擬邊界元模型的添加實(shí)例Fig.5 3 adding methods for the virtual surface elements
2.1 Alaoui剛性半球定常速垂向入水[15]
本文使用Abaqus/CAE程序建立了相應(yīng)的邊界元模型。半球的縱深H=77.34 mm,根據(jù)Nisewanger提供的半球入水自由液面升高公式計(jì)算可得在T1時(shí)刻半球恰恰完全浸沒(méi)于水(半球底面恰與等效自由液面平齊)。對(duì)于18 m/s的情況來(lái)講,T1=2.61 ms;對(duì)于20 m/s的情況,T1=2.41 ms。試驗(yàn)分別給出了6 ms和5 ms的時(shí)歷曲線(xiàn),因此可以斷定T1之后的數(shù)值是完全浸沒(méi)之后的數(shù)值(試件最高點(diǎn)低于等效自由液面),但是由于入水情況的瞬時(shí)性,半球周?chē)囊后w被拍擊而濺向四周,在浸沒(méi)之后的短時(shí)間內(nèi)是不會(huì)將水下物體覆蓋,從以往的試驗(yàn)中也可以觀(guān)察到這種現(xiàn)象。這就提出了新的問(wèn)題:T1時(shí)刻之后的濕表面的構(gòu)成問(wèn)題,即要探尋一種合適的水動(dòng)力模型(邊界元模型)來(lái)精確T1時(shí)刻之后的壓力計(jì)算。
根據(jù)圖5提出的三種添加邊界元模型的方法進(jìn)行砰擊載荷求解,計(jì)算結(jié)果如圖6-8所示。
圖6 18 m/s半球入水受力時(shí)歷Fig.6 The force-time curve of the hemisphere for 18 m/s
圖7 20 m/s半球入水受力時(shí)歷Fig.7 The force-time curve of the hemisphere for 20 m/s
可見(jiàn)使用邊界元模型2的虛擬面元形式計(jì)算結(jié)果與實(shí)驗(yàn)匹配良好,所以對(duì)于半球入水“后時(shí)段”的砰擊載荷計(jì)算可以采用該添加虛擬面元的方法進(jìn)行計(jì)算。
(1)本文基于三維邊界元方法,編寫(xiě)了可用于計(jì)算軸對(duì)稱(chēng)結(jié)構(gòu)入水砰擊載荷的計(jì)算程序,在自由液面處理方面繼承Wagner液面抬升理論,引入“浸深因子”Cw以確定液面抬升高度,簡(jiǎn)化了自由液面處理方式,提高了程序計(jì)算效率。
(2)算法中考慮了網(wǎng)格運(yùn)動(dòng),在自由液面附近對(duì)網(wǎng)格進(jìn)行截?cái)嘀貥?gòu)處理,精準(zhǔn)了每一時(shí)間步對(duì)應(yīng)的水下濕面積。
(3)使用該程序算法計(jì)算了Alaoui圓球體入水載荷,載荷計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)匹配良好,驗(yàn)證了本文計(jì)算程序的可行性與準(zhǔn)確性;并從另一方面說(shuō)明了線(xiàn)性自由液面對(duì)于砰擊載荷計(jì)算有較好的適用性。
(4)提出了入水“后時(shí)段”的砰擊載荷計(jì)算方法,擴(kuò)大了傳統(tǒng)的砰擊載荷時(shí)歷計(jì)算范圍。
圖8 半球阻力系數(shù)時(shí)歷Fig.8 The resistance-coefficients-time curve of the hemisphere
[1]Wu G X,Sun H,He Y S.Numerical simulation and experimental study of water entry of a wedge in free fall motion[J]. Journal of Fluids and Structures,2004,19:277-289.
[2]Wu G X.Numerical simulation of water entry of twin wedge[J].Journal of Fluids and Structures,2006,22:99-108.
[3]Xu G D,Duan W Y,Wu G X.Numerical simulation of oblique water entry of an asymmetrical wedge[J].Ocean Engineering,2008,35(16):1597-1603.
[4]Scolan Y M,Korobkin A A.Three-dimensional theory of water impact:Part 1.Inverse Wagner problem[J].Journal of Fluid Mechanics,2001,440:293-326.
[5]Korobkin A A,Scolan Y M.Three-dimensional theory of water impact:Part 2.Linearized Wagner problem[J].Journal of Fluid Mechanics,2006,549:343-373.
[6]Faltinsen O M,Chezhan M.A generalized wagner method for three-dimensional slamming[J].Journal of Ship Research, 2005,49(4):279-287.
[7]Von Karman T.The impact on seaplane floats during landing[R].NACA TN 321.Oct.1929.
[8]戴遺山,段文洋.船舶在波浪中運(yùn)動(dòng)的勢(shì)流理論[M].北京:國(guó)防工業(yè)出版社,2008:11-35.
[9]Prediction of impact pressure,forces,and moments during vertical and oblique water entry[R].NSWC.15,January,1977.
[10]Nisewanger C R.Experimental determination of pressure distribution on a sphere during water entry[J].NAVWEPS 7808, 1961.
[11]Faltinsen O M,Landrini M,Greco M.Slamming in marine applications[C].Journal of Engineering Mathematics,2004,48: 187-217.
[12]Krobkin A A,Pukhnachou V V.Initial stage of water impact[J].Annual Review Fluid Mechanics,1988,20:159-185.
[13]Takagi K,Dobashi J.Influence of trapped air on the slamming of a ship[J].Journal of Ship Research,2003,47(3):187-193.
[14]Sharov V F.Ship bottom impact upon wave[J].Sudostroyeniye,1958,4:5-9.
[15]Aboulghit EI Malki Alaoui,Alain Neme,Nicolas Jacques.Numerical and experimental studies of simple geometries in slamming[J].International Journal of Offshore and Polar Engineering,2011,21(3):216-224.
Numerical calculation of the slamming load for axisymmetric geometries during the later entry time based on generalized Wagner theory
ZHAO Jiu-long1,2,YAN Fa-suo1,FAN Lei3
(1.Harbin Engineering University,Harbin 150001,China;2.Shanghai Salvage Co./Ocean C&I,Shanghai 200129,China; 3.Texas A&M University,College Station,Texas 77840,USA)
A numerical code used for calculating the slamming load for axisymmetric geometries based on 3D BEM was mainly introduced.The 3D mechanical model and Wagner’s theory which considering the uplift free surface were both used in the numerical code.Cw which was called wetting factor was introduced to ensure the exact height of the uplift free surface.At the same time,the grid’s movement was also considered and along the interface between the geometry and the free surface one group of grids were modified for a precise wet area.The Nonlinear Bernoulli equation considering the grid’s movement was used for calculating the hydrodynamic pressure and the force could be obtained through integration.Besides,the conception of Virtual surface element was brought here so that the slamming load when the top point of the falling object was under the free surface could be achieved.That is to say,the traditional calculated time process was extended.At the same time,Alaoui slamming test was used to test the code’s validity.
axisymmetric geometries;slamming load;uplift of the free surface; grid’s modification;virtual surface element
TV131.2
A
10.3969/j.issn.1007-7294.2016.07.007
1007-7294(2016)11-1412-08
2016-04-07
趙九龍(1988-),男,工程師,E-mail:tcjiulong@sina.com;閆發(fā)鎖(1977-),男,教授。