基于移動(dòng)粒子半隱式方法的液柱倒塌計(jì)算

基于移動(dòng)粒子半隱式方法的液柱倒塌計(jì)算

自由地表流動(dòng)是自然界的一種常見(jiàn)現(xiàn)象。由于存在運(yùn)動(dòng)邊界，邊界位置只能在初始階段確定，然后需要作為計(jì)算結(jié)果的一部分進(jìn)行求解。因此，自由地表流動(dòng)的數(shù)值模擬一直是計(jì)算數(shù)學(xué)領(lǐng)域的研究重點(diǎn)，通常需要特殊處理。在傳統(tǒng)方法中，mac方法和vof（vluid）方法被廣泛使用。近年來(lái)，一種新的數(shù)值解算方法（mesh自由方法）引起了國(guó)際計(jì)算界的關(guān)注。在傳統(tǒng)的網(wǎng)格方法中，由于網(wǎng)絡(luò)切割難以實(shí)施，網(wǎng)格突變和網(wǎng)絡(luò)移動(dòng)等問(wèn)題的出現(xiàn)，而使用無(wú)網(wǎng)格方法，解算方法具有明顯的優(yōu)勢(shì)。因此，使用無(wú)網(wǎng)格方法進(jìn)行模擬自由表面流動(dòng)相對(duì)簡(jiǎn)單?，F(xiàn)在，許多無(wú)網(wǎng)格方法已經(jīng)研究了自由表面流動(dòng)。在這項(xiàng)工作中，使用了一種半隱形無(wú)網(wǎng)格方法（pms）方法來(lái)模擬和分析自由表面流動(dòng)。1求解流體的運(yùn)動(dòng)學(xué)模型MPS方法是日本東京大學(xué)Koshizuka等學(xué)者于20世紀(jì)90年代提出的一種求解不可壓縮流動(dòng)的無(wú)網(wǎng)格數(shù)值解法.其基本思想是在求解區(qū)域內(nèi)部及其邊界上布置一系列離散的粒子來(lái)表示宏觀的流體,流體運(yùn)動(dòng)控制方程中的各項(xiàng)微分算子以粒子間相互作用的形式表達(dá).每個(gè)粒子都有相應(yīng)的位置、質(zhì)量、動(dòng)量及能量等信息,這些信息通過(guò)引入核函數(shù)由其鄰點(diǎn)粒子信息近似獲得.在Lagrange坐標(biāo)系下追蹤各個(gè)粒子的運(yùn)動(dòng)變化規(guī)律,從而獲得整個(gè)流場(chǎng)的流動(dòng)信息.1.1u3000流體動(dòng)力黏度方程由于MPS方法求解問(wèn)題的拉哥朗日特性,故給出拉哥朗日形式下的流體運(yùn)動(dòng)控制方程:質(zhì)量守恒方程Dρ/Dt+ρ??u=0(1)Dρ/Dt+ρ??u=0(1)動(dòng)量守恒方程ρDuDt=-?p+μ?2u+ρf(2)ρDuDt=??p+μ?2u+ρf(2)式中:u為流體速度矢量;ρ為流體密度;p為壓力;μ為流體動(dòng)力黏性系數(shù);f為作用在流體單位質(zhì)量上的外力矢量.方程左端均以拉哥朗日時(shí)間微分形式表達(dá),對(duì)流項(xiàng)包含于其中.1.2粒子作用半徑re在MPS方法中,流體運(yùn)動(dòng)控制方程中的各項(xiàng)微分算子均以粒子相互作用的形式表達(dá).每個(gè)粒子通過(guò)核函數(shù)與其鄰點(diǎn)粒子發(fā)生相互作用,在MPS方法中所使用的核函數(shù)為w(r)={rer-1r＜re0r≥re(3)式中:r是粒子i和j之間的距離,r=|rj-ri|;re稱為粒子作用半徑.核函數(shù)具有緊支集特性,對(duì)于粒子間距離r≥re的粒子相互作用可以忽略,即粒子僅與位于其粒子作用半徑re范圍內(nèi)有限數(shù)目的鄰點(diǎn)粒子發(fā)生相互作用.定義i粒子的粒子數(shù)密度為?n?i=∑j≠iw(|rj-ri|)(4)則單位容積中的粒子數(shù)為?ρn?i=?n?i∫Vw(r)dv(5)流體密度等于單位容積中的粒子數(shù)與粒子質(zhì)量的乘積.假設(shè)各個(gè)粒子的質(zhì)量相等,則由式(5)可以證明流體密度與粒子數(shù)密度成正比.對(duì)于不可壓縮流體,流體密度保持定常相當(dāng)于其粒子數(shù)密度維持一個(gè)常數(shù),即n0=∑j≠iw(|r0j-r0i|)(6)設(shè)φi=φ(ri),表示位于ri粒子的物理參數(shù)標(biāo)量值,兩相鄰粒子i、j之間的梯度矢量可以表示為(φj-φi)(rj-ri)/|rj-ri|2,則i粒子的梯度矢量就可以通過(guò)其與其他粒子的梯度矢量加權(quán)平均得到,即??φ?i=dn0∑j≠iφj-φi|rj-ri|2(rj-ri)w(|rj-ri|)(7)式中:d為求解問(wèn)題的空間維數(shù).拉普拉斯算子?2φ的模擬通過(guò)擴(kuò)散方程由粒子與其鄰點(diǎn)粒子之間物理量遷移獲得,以粒子作用形式表示為??2φ?i=2dλn0∑j≠i(φj-φi)w(|rj-ri|)(8)式中λ=∫Vw(r)r2dv∫Vw(r)dv(9)λ的引入是使?2φ用粒子作用模型模擬得到的結(jié)果與擴(kuò)散方程解析解相一致.上述各模型中的粒子作用半徑re并不要求一致.1.3u3000固壁面粒子壓力公式無(wú)網(wǎng)格粒子方法目前尚處于高速發(fā)展期,各類方法中仍存在不少問(wèn)題亟待解決,其中之一就是程序的穩(wěn)定性問(wèn)題.MPS方法應(yīng)用于不可壓縮流體,在模擬固壁面時(shí)固壁靜止不動(dòng),當(dāng)流體與固壁發(fā)生沖擊碰撞時(shí),固壁面粒子壓力在瞬時(shí)可能會(huì)達(dá)到一個(gè)極大值.計(jì)算表明,這會(huì)給程序帶來(lái)不穩(wěn)定.為了解決這個(gè)問(wèn)題,本文中引入可壓縮因子α,并假設(shè)流體密度與壓力成線性關(guān)系,即ρn+1iρ0=1+αpn+1i(10)由于流體密度與粒子數(shù)密度成正比,故上式可寫為nn+1in0=1+αpn+1i(11)可壓縮因子并無(wú)特殊的物理意義,只是為了提高M(jìn)PS方法計(jì)算的穩(wěn)定性而引入的一個(gè)參數(shù).在本文計(jì)算中取α=10-7.當(dāng)程序中引入可壓縮因子以后,計(jì)算的穩(wěn)定性大大提高.1.4u3000壓力泊松方程的求解MPS方法采用半隱式時(shí)間推進(jìn)算法.在每個(gè)時(shí)層內(nèi),首先顯式計(jì)算動(dòng)量方程中的黏性項(xiàng)及源項(xiàng),得到粒子速度及位置的估算值,然后隱式求解連續(xù)方程獲得粒子速度及位置的修正值,最終得到粒子在下一時(shí)層的速度及位置.通過(guò)對(duì)各個(gè)粒子在各個(gè)時(shí)層運(yùn)動(dòng)規(guī)律的追蹤,獲得整個(gè)流場(chǎng)流體的流動(dòng)信息.具體算法如下所述.(1)給定粒子的初始分布,包括給定粒子的速度及位置的初場(chǎng):u0i,r0j.(2)在時(shí)間步長(zhǎng)Δt之后,顯式計(jì)算動(dòng)量方程中的黏性項(xiàng)及源項(xiàng),得到粒子速度及位置的估算值,即u*i=uni+Δt(ν?2uni+f)(12)r*i=rnj+Δtu*i(13)式中:ν為流體運(yùn)動(dòng)黏性系數(shù).(3)根據(jù)估算得到的粒子位置r*j計(jì)算此時(shí)的粒子數(shù)密度〈n*〉i,這個(gè)粒子數(shù)密度與流體粒子數(shù)密度常數(shù)n0有偏差,因此對(duì)各個(gè)粒子的粒子數(shù)密度進(jìn)行如下修正n*+n′=n0(14)同時(shí)由連續(xù)性方程得-??u′i=1Δtρn+1i-ρ*iρ0=1Δtnn+1i-?n*?in0(15)式中:粒子速度修正值由動(dòng)量方程隱式計(jì)算得到u′=-Δtρ?pn+1(16)聯(lián)立式(11)、式(14)～式(16),得到以下壓力泊松方程??2pn+1?i=-ρΔt2{?n*?i-n0n0-αpn+1i}(17)上式左端項(xiàng)以拉普拉斯模型替代,可以離散得到一組代數(shù)方程組,求解該方程組即可得到各粒子的壓力值pi.(4)由計(jì)算得到壓力值后,再根據(jù)式(16)計(jì)算各粒子的速度修正值,即u′=-Δtρ?pn+1(5)得到粒子在下一時(shí)層的速度及位置分布un+1i=u*i+u′i(18)rn+1i=r*i+u′iΔt(19)由上述算法可以看出,MPS方法具有物理概念清晰、算法簡(jiǎn)單易于實(shí)現(xiàn)的特點(diǎn).1.5自由表面上粒子數(shù)密度的計(jì)算傳統(tǒng)的數(shù)值方法求解流體力學(xué)問(wèn)題時(shí)對(duì)邊界的處理往往比較困難,但MPS方法對(duì)邊界的處理則相對(duì)比較簡(jiǎn)單.使用MPS方法可以很容易判定自由表面.由于在自由表面的外部區(qū)域沒(méi)有粒子存在,因此自由表面上粒子的粒子數(shù)密度較小.計(jì)算中凡滿足?n*?i＜βn0(20)條件的粒子即可認(rèn)為是自由表面粒子,式中β為一個(gè)小于1的常數(shù).對(duì)于固體壁面則可以很簡(jiǎn)便地以固定粒子表示,這些固定粒子速度總為0.為使流體粒子不進(jìn)入固壁區(qū)域,固壁粒子中靠近流體的第一層粒子需與流體粒子一起參與壓力計(jì)算,固壁粒子的層數(shù)則須滿足計(jì)算粒子團(tuán)密度所需的層數(shù).2流體物性參數(shù)液柱倒塌算例是一個(gè)模擬流體自由表面流動(dòng)的經(jīng)典算例,其易于建模,且存在實(shí)驗(yàn)數(shù)據(jù),因此被許多算法采用以考核程序的有效性.在VOF及SPH等方法中均以此作為算例對(duì)自由面流動(dòng)進(jìn)行了模擬計(jì)算.Koshizuka等在文獻(xiàn)中使用MPS方法對(duì)液柱倒塌進(jìn)行了成功的模擬計(jì)算,本文采用引入可壓縮因子后改進(jìn)的MPS方法程序?qū)σ褐顾憷M(jìn)行了計(jì)算.計(jì)算模型取于文獻(xiàn)中,如圖1所示,在初始時(shí)刻,矩形液柱以可移動(dòng)的擋板支持,靜止于容器的左端.液柱底部寬度L=0.146m,高h(yuǎn)=0.292m.將擋板迅速抽出,液柱在重力的作用下開(kāi)始倒塌,液體流向容器右壁并與之碰撞,在碰撞速度較大的情況下,將會(huì)有液滴濺出,隨后液滴又會(huì)下落至容器中.計(jì)算中取流體物性參數(shù)為:流體密度ρ=1.0×103kg·m3;動(dòng)力黏性系數(shù)μ=0.5kg/(m·s);重力加速度g=9.8m/s2.以MPS方法進(jìn)行數(shù)值模擬,液柱以648個(gè)粒子表示,粒子均勻布置于求解區(qū)域,則初始布置時(shí)相鄰兩粒子之間距離l0=0.008m.固壁面以3層共474個(gè)粒子表示,粒子初始布置如圖2所示.計(jì)算中各粒子作用模型中的粒子作用半徑分別取為粒子數(shù)密度及梯度模型中re=2.1l0,拉普拉斯模型中re=4.0l0.自由表面系數(shù)取β=0.97.從液柱倒塌開(kāi)始時(shí)刻t=0計(jì)算到t=1.0s時(shí),在IntelCeleronCPU2.00GHz主頻的計(jì)算機(jī)上計(jì)算時(shí)間約為1.5h.圖3中給出了計(jì)算中某幾個(gè)時(shí)間點(diǎn)的模擬結(jié)果.從圖3可以看出,液柱倒塌以后液柱在容器底部向右端流動(dòng),在0.3s時(shí)與右壁面發(fā)生碰撞,流體沿右壁面上升,并在0.6s左右開(kāi)始下落,0.8s時(shí)下落流體與容器中的流體發(fā)生沖擊,產(chǎn)生波動(dòng)向左端壁面推進(jìn)并在1.0s左右與左端壁面發(fā)生碰撞.使用MPS方法對(duì)液柱倒塌算例進(jìn)行了成功的模擬,方法簡(jiǎn)單易行,對(duì)自由面的判定可以通過(guò)計(jì)算粒子團(tuán)密度簡(jiǎn)單地獲得.文獻(xiàn)中給出了由實(shí)驗(yàn)得到的液柱倒塌實(shí)驗(yàn)中與右壁面碰撞之前液柱前沿至左端固壁的距離z與時(shí)間的關(guān)系.與文獻(xiàn)中形式一致,對(duì)MPS方法計(jì)算結(jié)果進(jìn)行無(wú)量綱化處理,圖4中示出了MPS方法計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果之間的比較.從圖4可以看出,MPS方法數(shù)值解與實(shí)驗(yàn)結(jié)果吻合得較好,圖中數(shù)值解與實(shí)驗(yàn)結(jié)果之間的偏差可以認(rèn)為是實(shí)驗(yàn)中流體與固壁之間存在摩擦阻力,而計(jì)算中忽略該阻力影響的緣故.文獻(xiàn)中由SOLA-VOF方法得到的數(shù)值解也示于圖4中,可以看出,MPS方法得到的結(jié)果與SOLA-VOF方法基本吻合且更接近于實(shí)驗(yàn)結(jié)果.考慮到MPS算法只需計(jì)算粒子數(shù)