畢業(yè)論文范文——基于SPH方法的三維矩形容器液體晃動(dòng)與沖擊數(shù)值模擬

1、畢業(yè)論文基于sph方法的三維矩形容器液體晃動(dòng)與沖擊數(shù)值模擬姓 名：專 業(yè)：航空電子班 級(jí)：完成日期： 指導(dǎo)教師: 摘要:基于sph方法，采用pam-crash軟件對(duì)無(wú)阻尼板和裝有阻尼板的矩形容器在加速度激勵(lì)作用下，容 器內(nèi)液體的晃動(dòng)與沖擊進(jìn)行了三維數(shù)值模擬，并將測(cè)試點(diǎn)的計(jì)算壓力及液體晃蕩運(yùn)動(dòng)模擬與試驗(yàn)進(jìn)行了對(duì)比, 吻合較好。與cfd方法相比,sph方法計(jì)算所得結(jié)果與試驗(yàn)更為接近，從而驗(yàn)證了 sph方法在求解三維液 體大幅晃動(dòng)問(wèn)題方面的準(zhǔn)確性。分析表明：阻尼板的安裝對(duì)液體的晃動(dòng)特性具有顯著的影響,遲滯特性所表 現(xiàn)的壓力雙峰特征在大充液比情況下消失。關(guān)鍵詞:sph ;自由液面；晃動(dòng)；沖擊；三維數(shù)值

2、模擬1引言帶自由液而的液體大幅晃動(dòng)問(wèn)題廣泛存在于航空航天、石油化工和交通運(yùn)輸?shù)阮I(lǐng)域，一 直受到各國(guó)工稈研究人員的關(guān)注。液體大幅晃動(dòng)具有高度的非線性，描述晃動(dòng)的一組方程呈 強(qiáng)利合性，自由液面邊界條件未知。對(duì)于如此復(fù)雜的一個(gè)問(wèn)題，傳統(tǒng)的數(shù)值計(jì)算方法如mac 法、vof法以及常用的有限元法均存在著一些缺陷：mac法對(duì)硬件要求高；vof法需耍高精 度的偏微分方程離散格式；有限元方法在處理液體大變形問(wèn)題時(shí)會(huì)出現(xiàn)網(wǎng)格的扭曲纏繞從而 導(dǎo)致計(jì)算的中止。光滑粒子流體動(dòng)力學(xué)(sph)方法是由lucy和monaghan111提出的。該方法的計(jì)算對(duì)彖是 空間運(yùn)動(dòng)的粒子，采用拉格朗口方法描述。無(wú)網(wǎng)格粒子形式以及拉格朗口

3、性質(zhì)的sph方法避 免了大變形時(shí)網(wǎng)格發(fā)生畸變的問(wèn)題，非常適用于液體大幅晃動(dòng)問(wèn)題的研究。自monaghan121首 次將sph方法應(yīng)用于口由表面水波的流動(dòng)模擬以來(lái)，一些國(guó)內(nèi)外學(xué)者開(kāi)展了對(duì)sph方法的應(yīng) 用研究oiglesias等i?叫采用sph方法模擬二維矩形水槽的縱搖，并與試驗(yàn)進(jìn)行了對(duì)bt；yonghwan kim利用此方法進(jìn)行了載液貨船液艙液體大幅晃動(dòng)及與貨船運(yùn)動(dòng)相耦合的數(shù)值模擬。陳正云 等采用sph方法對(duì)二維液艙內(nèi)液體大幅晃蕩問(wèn)題進(jìn)行了數(shù)值模擬，將計(jì)算得到的波面圖與 試驗(yàn)結(jié)果進(jìn)行了對(duì)比，吻合較好；馬利和梁龍河運(yùn)用光滑粒子流體動(dòng)力學(xué)耦合有限元的方 法進(jìn)行了碰撞和穿透問(wèn)題的模擬。目前，相對(duì)于二

4、維數(shù)值模擬，基于sph方法的液體大幅晃 動(dòng)三維數(shù)值模擬研究比較少。本文釆用sph方法進(jìn)行三維矩形容器液體晃動(dòng)與沖擊的數(shù)值模 擬，進(jìn)一步探討和驗(yàn)證了采用sph方法進(jìn)行帶自由液血的液體大幅晃動(dòng)數(shù)值計(jì)算的準(zhǔn)確性。2 sph數(shù)值方法2.1 sph方法sph方法通過(guò)構(gòu)造一個(gè)近似場(chǎng)函數(shù)/(兀»來(lái)求解一批具有一定質(zhì)量的粒子在任意時(shí)間 的速度和能量，表達(dá)式如下：si式中，兀和片為坐標(biāo)向量，w(x-x',h)為光滑核函數(shù)，力為決定光滑核函數(shù)影響域尺寸的光 滑長(zhǎng)度。一般取monaghan提出的分段三次b樣條函數(shù)作為核函數(shù)，如式(2)所示：2 . 1 3r h r0 r < 32w(r,h)

5、 = a(/< h2-/?)31<r<2 (2)0r>2式屮，r是粒子i和j之間的距離與光滑長(zhǎng)度的比值，6為與維數(shù)相關(guān)的系數(shù)，在一維、二維、三維空間中，a”分別為1/力，15/7刃幾3/2tz7f。通過(guò)對(duì)式(1)離散化，可得粒子i處的粒子近似式如式(3)所示：/(兀)=工一/(xy)-w(x -xj、k)7=1 p)式中，廠和(./= 1,2,., n)為粒子./的密度和質(zhì)量，n為粒子j影響域里的粒子總數(shù), 利用核函數(shù)的可微性，可得函數(shù)空間導(dǎo)數(shù)的粒子近似式：“ 1yiv /(£» =工s vw, - xj,/2) >1 pj從式(3)和(4)可

6、以看出，粒子近似式將函數(shù)及導(dǎo)數(shù)的連續(xù)積分表達(dá)式轉(zhuǎn)換成任意排 列的粒子的離散化求和，從而實(shí)現(xiàn)了空間的無(wú)網(wǎng)格化。2.2控制方程連續(xù)介質(zhì)流體動(dòng)力學(xué)ns方程分別由式(5)和式(6)描述：牛=_殲況(5)dtdu1 jj 一廠= v p h+ f(6)dtppdxf式中，p為液體密度，弘是速度向量，p為壓強(qiáng),“為液體的動(dòng)粘系數(shù)，f為體力。通過(guò) 對(duì)上述方程的轉(zhuǎn)換，可以得到方程的sph表達(dá)式如下：(7)du： a p p：dwy才專噲+才3材f式中，此 tuf, 為人工粘度，用于獲得穩(wěn)定的計(jì)算，具體表達(dá)式如下:nu =(9)a 勺 +0口0；pij0,h.v -x 式中，札=八，1 1 12xjj +02，

7、cii = 2 c + cj) , pij =q(pi + q)，% =-(/2, +hj) , vj =匕一 vj ,xj =xf -x. o q口和0口為常數(shù)，一般取1左右，本文的取值均為1, c為聲速，v為粒子的 速度矢量。對(duì)于不可壓縮流，通過(guò)利用人工壓縮性i可以將其視為可壓縮流，粒子i的壓強(qiáng)門由流 體狀態(tài)方程求得：pi = by -1)(io)a)式中，門為液體初始密度，常數(shù)/ = 7, b為表征粒子聲速的參數(shù)，b = 100p()ghy, h( 為容器內(nèi)液面高度。運(yùn)用sph方法進(jìn)行不可壓縮流計(jì)算時(shí)，時(shí)間步長(zhǎng)與聲速成反比。為了使顯式時(shí)間積分具 有穩(wěn)定性，采用cfl條件進(jìn)行約束】，即要求

8、時(shí)間步長(zhǎng)與粒子間距離的最小值成正比，在 sph中即轉(zhuǎn)化為時(shí)間步長(zhǎng)與最小光滑長(zhǎng)度成比例：(11)har = min()c£2.3固壁邊界條件為了阻止粒子穿透固壁邊界，在sph方法中，通過(guò)在固壁邊界上分布一組虛粒子用于對(duì) 鄰近邊界的實(shí)粒子作用排斥力，從而阻止鄰近邊界的這些粒子非物理穿透邊界，虛粒子和相 鄰粒子之間的作用力沿著兩粒子的中心連線方向。力的大小如式(12)所示：<1 厲(12) 乞1%式中，d、叫和直為可調(diào)參數(shù)。文中£的取值為0ghc, /?,和直一般取12和4。©為當(dāng) 前時(shí)刻虛粒子和鄰近邊界實(shí)粒子的間距，心為截止半徑：若太大，會(huì)造成初始階段，實(shí)粒 子

9、就受到排斥力的作用，擾亂粒子的分布；若太小，粒子可能穿過(guò)邊界也感受不到排斥力 作用，其取值一般與粒子初始間距的大小相近，算例中取其值為粒子的初始間距0.0076mo3數(shù)值算例采用sph方法，不考慮容器內(nèi)液體表面所受張力，分別對(duì)cniig試驗(yàn)"2】中無(wú)阻尼板和裝 有阻尼板的矩形容器在平動(dòng)加速度激勵(lì)作用下，容器內(nèi)液體的晃動(dòng)與沖擊進(jìn)行數(shù)值計(jì)算，并 與試驗(yàn)結(jié)果及cfd計(jì)算結(jié)果進(jìn)行了對(duì)比。試驗(yàn)中在矩形容器的底面和側(cè)壁布置兩個(gè)壓力 測(cè)試點(diǎn)p1和p2,容器的長(zhǎng)、寬、高分別為500mm、400mm、400mm,容器的主視圖及壓力 測(cè)試點(diǎn)布置如圖1所示。充液比為0.5,容器固定于一輛輕型卡車上，車由靜

10、止加速至40km/h 再減速至靜止?fàn)顟B(tài)，小車的實(shí)測(cè)加速度曲線及濾波后用于數(shù)值計(jì)算的曲線如圖2所示，表1 給出了主要的材料屬性參數(shù)。表1材料屬性屬性名稱數(shù)值矩形容器及阻尼板密度/ (kg/m3)7830矩形容器及阻尼板楊氏模量/(pa)2.07ell矩形容器及阻尼板泊松比0.3液體密度/ (kg/m3)998液體楊氏模量/ (pa)2.2e9圖1矩形容器示意圖|ii|5i:0 亠二二j_未濾波1015圖3比較直觀地顯示了在加速度激勵(lì)作用下自由液而晃動(dòng)的三維模擬及其對(duì)應(yīng)的對(duì)稱而 上的粒子速度矢量圖，從圖中可以清楚看出液面的翻卷及波浪的破碎，當(dāng)側(cè)壁處粒子速度方 向垂直指向壁面，且運(yùn)動(dòng)速度比較大時(shí)，會(huì)

11、對(duì)容器壁產(chǎn)生很大的沖擊壓力。圖4給出了容器 運(yùn)動(dòng)過(guò)程中三種典型的液體晃蕩現(xiàn)象一駐波、行進(jìn)波和水躍。圖5為有無(wú)阻尼板情況下，兩 種數(shù)值計(jì)算方法對(duì)晃動(dòng)的模擬與試驗(yàn)的對(duì)比，sph計(jì)算較真實(shí)地模擬出了試驗(yàn)中口由液面隨 著容器加速度變化，從平靜到大幅晃動(dòng)的過(guò)程，該算例屮，相對(duì)于cfd方法，sph方法模 擬的結(jié)果更接近于試驗(yàn)。圖6給出了試驗(yàn)、cfd方法及sph方法所得的測(cè)試點(diǎn)p1處的壓力時(shí)間歷程曲線。從圖 6 (a)中可以看到5s之后測(cè)試點(diǎn)的壓力所呈現(xiàn)的雙峰特征，這是液體晃動(dòng)過(guò)程中遲滯特性的 一種典型表現(xiàn)。其中前一個(gè)峰值是由于液體對(duì)容器壁的沖擊而產(chǎn)生的高脈沖壓力，時(shí)間短， 具有比較明顯的局部特性。這可以從

12、圖2中的加速度曲線分析出：5-10s這段時(shí)間內(nèi)加速度 變化大，容器壁主要承受水躍(圖4 (a)或者行進(jìn)波(圖4 (b)晃動(dòng)現(xiàn)彖引起的瞬時(shí)沖擊 壓力。后一個(gè)峰值也是由于沖擊產(chǎn)生的脈沖壓力，但是時(shí)間相對(duì)較長(zhǎng)，而且沖擊壓力峰值變 小，這主要是由于快速并且連續(xù)的液體壓力作用在沒(méi)有被完全淹沒(méi)的容器壁(圖4 (c),壓 力特性主要表現(xiàn)為連續(xù)性，而不像前一個(gè)峰值中所明顯表現(xiàn)出的瞬時(shí)性。當(dāng)給容器加裝阻尼 板后，液體對(duì)容器壁的沖擊力度降低，由于阻尼板的而積較大，阻尼板的安裝一定程度上可 以等效為增加了容器的充液比，在大充液比情況下，遲滯現(xiàn)象消失，因此圖6 (b)中沒(méi)有明 顯的雙峰特性，且壓力峰值減小，晃動(dòng)現(xiàn)象更

13、多的表現(xiàn)為駐波運(yùn)動(dòng)(圖4 (d)。另外，在無(wú) 阻尼板工況下，當(dāng)加速度減至0并基本趨于穩(wěn)定之后，p1處的壓力值開(kāi)始在初始靜態(tài)壓力值 處附近振蕩，并呈現(xiàn)出有阻尼的逐漸衰減的趨勢(shì)；在有阻尼板工況下，p1處壓力值則很快收 斂于初始靜態(tài)壓力值。通過(guò)將圖6中的sph計(jì)算結(jié)果與試驗(yàn)及cfd計(jì)算結(jié)果對(duì)比可以發(fā)現(xiàn)，sph方法計(jì)算所 得到的曲線較cfd方法不論在變化趨勢(shì)、壓力作用周期還是壓力峰值方面都與試驗(yàn)吻合更 好。計(jì)算峰值稍低于試驗(yàn)所測(cè)峰值，這是由于前面所提及的對(duì)小車的加速度激勵(lì)曲線進(jìn)行了 濾波，濾去了一部分波峰的緣故。對(duì)于8-10s時(shí)間段，兩張圖屮的計(jì)算曲線與試驗(yàn)吻合均不 是很理想，從數(shù)值計(jì)算所得曲線的斜率

14、可以看出壓力變化比較突然，這可能是由于試驗(yàn)中低 頻加速度激勵(lì)測(cè)量誤差所致i®,該時(shí)間段加速度測(cè)量值變化程度要比真實(shí)情況劇烈，導(dǎo)致數(shù) 值計(jì)算中的液體沖擊幅度較大，而在真實(shí)試驗(yàn)中，此時(shí)間段內(nèi)容器壁承受的壓力主要表現(xiàn)為 持續(xù)性的壓力，因而壓力曲線變化比較平緩。(a)無(wú)阻尼板(b)有阻尼板圖3自山液而晃動(dòng)的三維模擬及對(duì)應(yīng)的對(duì)稱而粒子速度矢屋圖瞬時(shí)壓力(a)水躍(無(wú)阻尼板)瞬時(shí)壓力(b)行進(jìn)波(無(wú)阻尼板)連續(xù)壓力(c)駐波(無(wú)阻尼板)連續(xù)壓力(d)駐波(有阻尼板)圖4典型晃蕩運(yùn)動(dòng)試驗(yàn)sph計(jì)算cfd計(jì)算sph計(jì)算試驗(yàn)cfd計(jì)算sobunhs8hmuthk 8huutth2 6生utth 40

15、tih 禺 70§h啜tuuiu(a)無(wú)阻尼板(b)冇阻尼板圖5液體晃動(dòng)的數(shù)值模擬與試驗(yàn)對(duì)比3000-12500-20001500-1000-500015時(shí)間/s試驗(yàn)-o-cfdil /sph計(jì)算£=二二(a)無(wú)阻尼板(b)有阻尼板圖6 sph數(shù)值計(jì)算與試驗(yàn)及cfd數(shù)值計(jì)算壓力曲線對(duì)比4結(jié)論采用sph方法，對(duì)矩形容器在平動(dòng)加速度激勵(lì)作用下自由液面的晃動(dòng)進(jìn)行了三維數(shù)值模 擬，通過(guò)與文獻(xiàn)屮試驗(yàn)數(shù)據(jù)的對(duì)比，發(fā)現(xiàn)計(jì)算所得到的壓力曲線在變化趨勢(shì)、壓力作用周期 和峰值方面都與試驗(yàn)吻合較好，此外，文中算例采用的sph方法計(jì)算所得結(jié)果比cfd方法 更接近于試驗(yàn)；整個(gè)晃動(dòng)過(guò)程模擬也與試驗(yàn)相近

16、，并且比較真實(shí)地模擬出了自由液血大幅晃 動(dòng)下采用有限元方法難以模擬的碎波現(xiàn)象。因此，sph方法進(jìn)行帶自由液面的液體大幅晃動(dòng) 計(jì)算是合理的、準(zhǔn)確的。本文為求解類似關(guān)于液體晃動(dòng)的問(wèn)題如飛機(jī)油箱燃油晃動(dòng)、航天器 推進(jìn)劑晃動(dòng)提供了有效的方法和手段。在本文的計(jì)算中，筆者為提高計(jì)算精度，所用粒子數(shù)達(dá)到18萬(wàn)，盡管釆用了 4個(gè)cpu 進(jìn)行并行計(jì)算，但是發(fā)現(xiàn)并行效率并不高，因此仍花費(fèi)了大量的計(jì)算吋間，如何提高并行效 率將是下一步所要探討的問(wèn)題。另外，文中進(jìn)行了容器增加阻尼板后的液體晃動(dòng)計(jì)算，阻尼板的數(shù)量、形狀、彈性和安 裝位置等因素對(duì)容器內(nèi)液體的晃動(dòng)及容器壁所受的沖擊壓力具有比較人的影響，關(guān)于容器防 晃阻尼板

17、優(yōu)化這方面的研究，今后將采用sph方法結(jié)合優(yōu)化算法開(kāi)展進(jìn)一步的工作。參考文獻(xiàn)：1 r a. gingold, j j monaghan. smoothed particle hydrodynamics-theory and application to non-spherical stars j. royal astronomical society, monthly notices, 1977, 181: 375-389.2 j j monaghan. simulating free surface flow with sph j journal of computational physi

18、cs, 1994, 110: 399-406.3 a. souto-iglcsias, l perez-rojas, r zamora rodriguez simulation of anti-roll tanks and sloshing lypc problems with smoothed particle hydrodynamics j ocean engineering, 2004, 31: 1169-1192.4 a. souto-iglesias, l delorme, l perez-rojas cl al. liquid moment amplitude assessment

19、 in sloshing typo problems with smooth particle hydrodynamics j ocean engineering, 2006, 33: 1462-1484.5 kim. yonghwan. experimental and numerical analyses of sloshing flows j journal of engineering mathematics, 2007, 58: 191-210.6 陳正云，朱仁慶，祁江濤.基于sph法的二維液體大幅晃蕩數(shù)值模擬j,船海工程，200& 37(2): 44-47.7 馬利，陶偉明，郭乙木等.sph耦合有限元方法的水射流彈塑性碰撞模擬j.浙江大學(xué)學(xué)報(bào)(工學(xué)版)， 200& 42(2): 259-263.8 梁龍河，董師舜，王政.剛體與sph耦合的三維數(shù)值方法研究j,高壓物理學(xué)報(bào)，200& 22(3): 275-280.9 l. i). libersky, a. g petschek, t c carney et al. high strain lagrangian hydrod