楔形體入水撞擊的自由面分析

楔形體入水撞擊的自由面分析

0在仿真中的應(yīng)用在許多應(yīng)用中，高速流固耦合的影響可以得到體現(xiàn)。被人們熟知的有甲板上浪、入水沖擊、波浪沖擊海洋結(jié)構(gòu)物或者海岸、雷彈入水及船舶出水再入水的砰擊作用。剛性楔形體的常速入水對砰擊的研究具有重要意義。最早前蘇聯(lián)的Dobrovol′skaya基于自相似方法的速度勢理論求解了楔形體常速入水,忽略了流體重力、水的可壓縮性及粘性,得到了物體表面的壓力分布和速度勢。Zhao和Faltinsen利用時域方法也求解了同樣的問題。他們利用更精細(xì)的方法,在高性能計算機(jī)的環(huán)境下利用邊界元方法對此類問題進(jìn)行了求解。Tveitnes對楔形體的常速入水進(jìn)行了試驗研究,并和仿真結(jié)果進(jìn)行了對比。楔形體入水問題還有很多值得解決的問題,例如,如何將2D相似解擴(kuò)展到3D,以及3自由度的楔形體自由入水和復(fù)雜2D結(jié)構(gòu)和3D結(jié)構(gòu)的入水,射流從物體表面分離等都是在計算中需要考慮的問題。本文以柯西積分作為積分方程,利用復(fù)數(shù)變量邊界元法求解了常速入水的相似解。在以往的求解中,為了避免在自由面射流處的震蕩,以及射流變長使網(wǎng)格增多、網(wǎng)格的尺寸減小而引起的數(shù)值發(fā)散問題,都是將射流剪掉,因為射流處壓力很接近大氣壓,所以對物體整體運(yùn)動幾乎沒有影響。將射流剪掉后的處理方法一般是在邊界的缺口補(bǔ)充一個垂直于物面的單元。本文在數(shù)值仿真的過程中不將射流切掉,而是當(dāng)射流變長變薄時利用淺水近似結(jié)合非線性自由面邊界條件進(jìn)行處理。仿真過程中,對不同底升角的楔形體相似解進(jìn)行了研究對比,并對涉及到的一些數(shù)值技術(shù)進(jìn)行了討論。1數(shù)值模擬1.1流體表面為固體邊界條件仿真涉及到的問題是底升角為γ的2D對稱楔形體垂直入水問題,如圖1所示。建立笛卡爾靜止坐標(biāo)系oxy,o點(diǎn)位于靜水面上并和初始時楔形體頂點(diǎn)重合,x正方向沿靜水面向右,y軸垂直向上。假設(shè)流體不可壓縮、無粘、無旋且忽略質(zhì)量和重力。引入速度勢φ和流函數(shù)ψ,由假設(shè)可知φ在流體域內(nèi)滿足拉普拉斯方程在楔形體的表面滿足流體不可穿越固體邊界條件式中:V為物體的運(yùn)動速度;n為物體表面的垂直法向向量,其方向指向流體域外,n=(nx,ny)=(-sinγ,cosγ)。運(yùn)動學(xué)邊界條件規(guī)定自由表面上的水質(zhì)點(diǎn)始終在自由表面上,動力學(xué)邊界條件規(guī)定自由面上的流體壓力始終為大氣壓力。因為沖擊持續(xù)的時間很短,流體重力的影響相在小時間上的展開是2階小量,所以可以忽略。令自由面上的大氣壓為零,得到歐拉形式的自由面邊界條件為控制面上的邊界條件滿足流體不穿越此剛性壁邊界條件1.2邊界條件迭代求解當(dāng)忽略流體重力時,常速入水的楔形體在不同時間段水面隆起的形狀都是相似的,可以利用無量綱變量來進(jìn)行求解,在無量綱相似坐標(biāo)系下入水物體在垂直方向上的位移始終為1。這種問題可以通過迭代法求解,假設(shè)在初始時刻自由液面是水平的且自由面上的速度勢為零,然后通過迭代求出自由面的形狀,直至收斂即可得到相似解。當(dāng)入水速度V是常數(shù)時,將速度勢寫成以下形式式中:ξ=x/Vt;η=y/Vt;?為相似坐標(biāo)系下的速度勢,也滿足拉普拉斯方程。物面邊界條件變?yōu)樵谧杂擅嫔?運(yùn)動學(xué)和動力學(xué)邊界條件變?yōu)楫?dāng)流體緊貼物面飛濺飆升時,剖面的形狀函數(shù)在(ξ,η)內(nèi)是多值的。在求解的過程中,為了避免自由液面上顆粒進(jìn)入物面而引起的數(shù)值錯誤。對右邊物面,將(ξ,η)順時針旋轉(zhuǎn)π/2-γ,得到新的坐標(biāo)(α,β)在新的坐標(biāo)系中,β是沿楔形體表面的。自由面的邊界條件保持不變。此時可通過迭代來求解。對自由面邊界條件直接迭代可能放大誤差導(dǎo)致程序發(fā)散。將自由面邊界條件寫成下面的形式如果將式(11)和式(13)的右邊看成是已知的,則上式就可以看成是1階常微分方程,利用1階常微分方程的解法得到此式即為仿真中迭代所用到的式子,式中,(α0,β0)為控制面和自由面的交點(diǎn)。當(dāng)又細(xì)又長的射流出現(xiàn)時,射流區(qū)域可以結(jié)合淺水近似求解,由淺水理論可知速度勢在射流上是線性變化的,取前2項由自由面邊界條件有?1=cosγ,同樣可以認(rèn)為在射流上β隨α也是線性變化的,結(jié)合自由面邊界條件可以得到速度勢和射流形狀的表達(dá)式式中的系數(shù)為由式(20)可知射流和物面的交點(diǎn)為(αc=cosγ,βc),得到在此點(diǎn)的A=βc,然后由式(19)得即為交點(diǎn)處的速度勢,在每一次的求解過程中(αs,βs)可以通過迭代的方法求出,然后從這一點(diǎn)開始,射流上所有點(diǎn)的新坐標(biāo)點(diǎn)和速度勢就可由淺水近似求出了。2數(shù)值方法2.1分布式akj求解算法?和χ的邊界值問題都需要在相似坐標(biāo)系下利用復(fù)速度勢在復(fù)數(shù)邊界元下求解,由柯西積分方程得式中:z=ξ+i×η;z0為流域內(nèi)的一點(diǎn),積分方向沿流體邊界的正向,將節(jié)點(diǎn)劃分為線性單元之后,可以將邊界的復(fù)速度勢寫成每個單元內(nèi)復(fù)速度勢的線性逼近組合。其中,Nj(z)為線性基函數(shù)將線性基函數(shù)Nj(z)帶入式(22)和式(23),然后讓z0逼近控制體邊界單元的每一個節(jié)點(diǎn)k時,有其中,Akj為矩陣方程的系數(shù)對于Akj表達(dá)式,當(dāng)k分別為j-1,j,j+1時需要利用極限的方法求出,另外,對于表達(dá)式Akj中的復(fù)對數(shù)求解可采用下列方法其中,θ(j+1,j)為z0與zj以及zj+i這2條線所夾的中心角,此時z0即為zk。上述方程可以拆分成下面的形式由邊界條件知,自由面上的速度勢和物面上的流函數(shù)是已知的,移到方程右邊。左邊是待求自由面上的流函數(shù)和物面上的速度勢。對于對稱楔形體的垂直入水,剛性壁面上的流函數(shù)C當(dāng)作零來處理。為了提高矩陣求解質(zhì)量,避免數(shù)值發(fā)散,當(dāng)2≤k≤N1或N2≤k≤N3時取矩陣方程的實部,當(dāng)N1≤k≤N2或N3≤k≤N時取方程的虛部。另外,控制面和自由面以及控制面與物面的交點(diǎn)上的速度勢和流函數(shù)都是已知的,它們都被移動到方程的右邊。每次求解矩陣方程后得到自由面上的流函數(shù)后即可利用式(14)和式(15)得到新的自由面和自由面上的速度勢。再求解矩陣方程,如此循環(huán)直至收斂為止。2.23建立標(biāo)準(zhǔn)條件在迭代時,有的網(wǎng)格單元變長,而有的則變短,且在數(shù)值積分的過程中自由面可能會出現(xiàn)毛刺,使得計算很難進(jìn)行下去。為了使相鄰單元長度的比例始終保持在1左右,需要對自由面網(wǎng)格進(jìn)行重新劃分。首先,以控制面和自由面的交點(diǎn)為起點(diǎn)建立坐標(biāo)和多邊形弧長的3次樣條公式式中:qi為從起點(diǎn)到任一自由面節(jié)點(diǎn)上的多邊形弧長;系數(shù)ai,bi,ci,di可由在每個節(jié)點(diǎn)上坐標(biāo)值、1次導(dǎo)數(shù)、2次導(dǎo)數(shù)以及3次導(dǎo)數(shù)相等求出。其次,基于式(29)和式(30),由高斯積分求出每個單元的弧長及從起點(diǎn)到任一節(jié)點(diǎn)的弧長il。然后建立弧長和自由面節(jié)點(diǎn)坐標(biāo)的3次樣條表達(dá)式最后,根據(jù)式(31)和式(32)得到需要的網(wǎng)格分布,在靠近物面的射流上劃分均勻的網(wǎng)格,在遠(yuǎn)離物面的自由面上網(wǎng)格大小按照弧長等比增大進(jìn)行劃分。在得到自由面上新的坐標(biāo)分布后,可以根據(jù)相同的方法求出速度勢在自由面上的分布。另外,在每一次的自由面更新后需要對物面重新劃分網(wǎng)格,網(wǎng)格的大小和射流上的網(wǎng)格保持在一個數(shù)量級。3底升角的影響首先對收斂性進(jìn)行測試,控制體大小取為(±m(xù)ax(10,10/tanγi),-10)。在物面上劃分均勻網(wǎng)格,網(wǎng)格的大小和自由面上射流部分保持一致。在劃分網(wǎng)格時保持相鄰單元的比例接近1,這樣可以避免數(shù)值發(fā)散?？刂泼嫔系木W(wǎng)格大小保持為物面單元的5倍左右。在測試的過程中,物面上和射流處單元從0.02到0.1分別進(jìn)行測試,發(fā)現(xiàn)當(dāng)網(wǎng)格大小為0.04到0.05內(nèi)時經(jīng)過若干次的迭代之后射流的形狀基本上保持不變,且收斂的形狀也幾乎相同。本文中網(wǎng)格的大小取0.04。為了避免自由面的更新過快,對小的底升角需要使用欠松弛因子。式中:f(i)和f(i+1)分別為第i和第i+1迭代步的值,可為速度勢,也可為坐標(biāo)值;μ為欠松弛因子。對底升角從20°到80°的情況分別進(jìn)行了計算仿真。楔形體表面歸一化壓力p/ρV2的分布為得到的楔形體表面壓力分布和相似解液面隆起分別如圖2和圖3所示。由圖2可知,楔形體入水時的力學(xué)特性和底升角有很大的關(guān)系,當(dāng)?shù)咨呛艽髸r,即楔形體很尖很薄,入水時受到的沖擊壓力很小,同時激起的射流飛濺也很小。而當(dāng)?shù)咨呛苄r,楔形體表面壓力的分布呈現(xiàn)中凹的趨勢,當(dāng)?shù)咨侵饾u減小時表面的壓力急劇增大,射流飛濺也又細(xì)又長。當(dāng)?shù)咨呛苄r入水可能會引起起點(diǎn)效應(yīng),而文中鑒于模型求解復(fù)雜化,沒有考慮這一效應(yīng)。4低通插裝式樁件中底升角的影響本文建立了楔形體多種底升角的入水砰擊模型,利用復(fù)變邊界元發(fā)結(jié)合自由面邊界條件迭代的方法求解了對稱楔形入水的相似解,模擬了垂向入水的力學(xué)行為,對影響砰擊過程的數(shù)值參數(shù),包括單元的劃分、控制體大小的選取以及插值等因素進(jìn)行了