蘇鐵熊 馬理強(qiáng)? 劉謀斌 常建忠

1 引言

液滴沖擊固壁面現(xiàn)象廣泛存在于工業(yè)生產(chǎn)的各個(gè)領(lǐng)域,如噴墨打印、噴涂印刷、輪機(jī)葉片腐蝕、燃油噴射霧化等[1,2].液滴沖擊固壁面是一種典型的不可壓縮自由表面流動(dòng)問(wèn)題,包含復(fù)雜的流動(dòng)現(xiàn)象和流動(dòng)過(guò)程,如液滴的鋪展、破碎、反彈、飛濺現(xiàn)象等,影響因素涉及液滴的密度、直徑、黏性、沖擊速度,以及固壁面的物理特性如溫度、表面粗糙度、壁面的浸潤(rùn)性等.因此,針對(duì)液滴沖擊固壁面問(wèn)題的機(jī)理研究對(duì)環(huán)境工程、微納米工程以及醫(yī)藥工程等有著十分重要的作用[3,4].

國(guó)內(nèi)外的許多研究人員通過(guò)理論分析、實(shí)驗(yàn)觀察和數(shù)值模擬等方法對(duì)液滴沖擊固壁面問(wèn)題進(jìn)行了研究.Worthington[5]首先系統(tǒng)地研究了液滴沖擊問(wèn)題.Roisman等[6]研究了液滴沖擊固壁面過(guò)程中接觸角的影響.Qiang等[7]采用光滑粒子動(dòng)力學(xué)方法研究了液滴在氣固交界面變形移動(dòng)問(wèn)題.Bussmann等[8]研究了液滴沖擊固壁面的產(chǎn)生的飛濺現(xiàn)象.Eggers等[9]研究了較大沖擊速度下的液滴沖擊固壁面問(wèn)題.Ellis等[10]研究了表面粗糙度對(duì)液滴沖擊的影響.Ma等[11]采用改進(jìn)的光滑粒子動(dòng)力學(xué)方法研究了液滴沖擊問(wèn)題.Yang等[12]研究了速度對(duì)液滴撞擊超疏水壁面行為特性的影響.Sikalo等[13]采用實(shí)驗(yàn)和數(shù)值模擬研究了液滴沖擊鋪展過(guò)程中接觸角的影響,Liu和Chang[4]采用粒子方法研究了多相流動(dòng)過(guò)程中接觸角和壁面浸潤(rùn)效應(yīng)的影響.

光滑粒子動(dòng)力學(xué)(smoothed particle hydrodynamics,SPH)方法[14,15]是一種自適應(yīng)拉格朗日型無(wú)網(wǎng)格粒子方法.在SPH方法中,SPH使用粒子離散及代表所模擬的介質(zhì),并且基于粒子體系估算和近似介質(zhì)運(yùn)動(dòng)的控制方程[16,17].作為一種拉格朗日方法,通過(guò)跟蹤粒子的運(yùn)動(dòng)確定物質(zhì)的運(yùn)動(dòng),不需要復(fù)雜的算法來(lái)追蹤諸如自由表面、移動(dòng)邊界及運(yùn)動(dòng)界面等運(yùn)動(dòng)特征,能夠自然地描述介質(zhì)的運(yùn)動(dòng)過(guò)程,避免計(jì)算對(duì)流或輸運(yùn),因此,特別適合模擬涉及自由表面流動(dòng)、運(yùn)動(dòng)界面、變形邊界和大變形問(wèn)題[18-21].

傳統(tǒng)的SPH方法精度和穩(wěn)定性較差,尤其是在邊界區(qū)域以及粒子分布不均勻的區(qū)域.為提高傳統(tǒng)SPH方法的精度和穩(wěn)定性,對(duì)密度和核梯度近似形式進(jìn)行了修正,此方法能夠保持粒子相互作用過(guò)程中總動(dòng)量始終守恒.采用黎曼解法的SPH流體控制方程,液滴與固壁面相互作用過(guò)程中,考慮到液滴的尺度和支持域內(nèi)粒子較少,本文采用了一種新型的粒子間相互作用力(IIF)模型來(lái)模擬表面張力的影響,避免了界面的尖角以及邊界處等粒子缺失和分布不均勻的地方曲率計(jì)算誤差較大的問(wèn)題.應(yīng)用改進(jìn)的SPH方法對(duì)液滴沖擊固壁面問(wèn)題進(jìn)行了數(shù)值模擬.

2 光滑粒子動(dòng)力學(xué)

2.1 SPH方法核近似與粒子近似

SPH方法使用一組任意分布的粒子來(lái)代表所模擬的介質(zhì)系統(tǒng),基于粒子體系離散控制流體運(yùn)動(dòng)的偏微分方程,并通過(guò)適當(dāng)?shù)淖儞Q對(duì)所獲得的SPH運(yùn)動(dòng)方程進(jìn)行求解,從而仿真流體系統(tǒng)的動(dòng)力學(xué)特征.利用SPH方法對(duì)偏微分方程的近似包括兩步：核近似和粒子近似.對(duì)于傳統(tǒng)的SPH方法,在某點(diǎn)(或某個(gè)粒子)i上,對(duì)任意函數(shù) f(x)的核近似(〈f(x)〉)可由下式定義

其中,f為三維坐標(biāo)向量x的函數(shù),Ω為包含x的積分區(qū)域,W是光滑函數(shù),h定義了光滑函數(shù)W的影響區(qū)域,稱為光滑長(zhǎng)度.

f(x)的離散形式的粒子近似是對(duì)相關(guān)粒子i支持域內(nèi)所有粒子進(jìn)行加權(quán)求和(如圖1)

式中N為粒子i的支持域內(nèi)所有粒子總數(shù),ρj和mj為粒子 j的密度和質(zhì)量,Wij為粒子 j對(duì)粒子i產(chǎn)生影響的光滑函數(shù).(2)式表明粒子i處的任一函數(shù)值可通過(guò)應(yīng)用光滑函數(shù)對(duì)其緊支域內(nèi)所有粒子相對(duì)應(yīng)的函數(shù)值進(jìn)行加權(quán)平均進(jìn)行近似.

圖1 二維空間SPH粒子近似示意圖 W為光滑函數(shù),支持域?yàn)棣蔴,S為計(jì)算區(qū)域Ω的表面

2.2 考慮黎曼解的SPH流體控制方程

液滴沖擊固壁面過(guò)程中,液滴和固壁面接觸區(qū)域粒子的材料性質(zhì)(如密度、壓強(qiáng)、速度等)梯度變化很大,會(huì)導(dǎo)致局部的非連續(xù)和間斷,利用傳統(tǒng)SPH的流體控制方程算法精度往往不高,在較大的We數(shù)沖擊作用下可能會(huì)導(dǎo)致計(jì)算結(jié)果的失真和計(jì)算中止.基于黎曼解法對(duì)于非連續(xù)問(wèn)題和間斷問(wèn)題求解的有效性,本文采用考慮黎曼解法的SPH流體控制方程：

2.3 密度和核梯度修正

在早期的文獻(xiàn)中,就近似格式而言,SPH通常被認(rèn)為是二階精度的方法[22].對(duì)(1)式右邊在x處進(jìn)行泰勒級(jí)數(shù)展開(kāi)

其中r為余項(xiàng).因?yàn)楣饣瘮?shù)滿足正則化條件和對(duì)稱性條件,(5)式可寫為

因此就核近似而言,SPH方法具有二階精度.

然而,(6)式未必對(duì)所有的情形都成立.例如,當(dāng)某個(gè)粒子臨近計(jì)算區(qū)域邊界時(shí),粒子的支持域與計(jì)算區(qū)域相交,支持域被邊界截?cái)?見(jiàn)圖2),正則化條件和對(duì)稱性條件都不能滿足,因此函數(shù)的SPH核近似不具有二階精度.

圖2 粒子的支持域與計(jì)算區(qū)域相交

利用SPH方法對(duì)偏微分方程進(jìn)行近似,其精度最終取決于離散形式的粒子近似.在傳統(tǒng)的SPH計(jì)算中,初始的粒子一般規(guī)則分布,質(zhì)量一致(粒子大小相等).隨著計(jì)算的進(jìn)行,粒子分布逐漸變得紊亂不規(guī)則,每個(gè)粒子的密度隨壓力應(yīng)力變化逐漸演變,因而每個(gè)粒子的大小也逐漸不一.因此傳統(tǒng)的SPH粒子近似格式很難保證一階甚至是0階連續(xù)性,因而不能精度再生線性函數(shù)甚至常數(shù).這也是導(dǎo)致傳統(tǒng)SPH方法粒子近似格式精度偏低的根本原因.

為提高傳統(tǒng)SPH方法的計(jì)算精度,參照文獻(xiàn)[23]對(duì)SPH方法中的密度和核梯度進(jìn)行了修正

對(duì)(1)式等號(hào)右邊項(xiàng)基于泰勒級(jí)數(shù)展開(kāi),可以得到

由于

可以得到一個(gè)新的核梯度表達(dá)式

粒子近似方程為

最終,核梯度修正模型為

2.4 粒子間相互作用的表面張力模型

在SPH方法中,考慮到所模擬液滴的尺度和支持域內(nèi)粒子數(shù)較少,因此,表面張力的作用尤其重要.對(duì)于表面張力的計(jì)算,研究者提出了多種求解的模型[24,25],大致可以分為兩類,即基于連續(xù)表面力(continuum surface force,CSF)的張力模型和基于原子/分子尺度的粒子間相互作用力(interparticle interaction force,IIF)模型.CSF模型通過(guò)求解界面曲率來(lái)計(jì)算表面張力,由于每一個(gè)時(shí)間步都需要求解界面曲率,因此計(jì)算的效率較低;且對(duì)于交界面的尖角以及邊界等粒子缺失和分布不均勻的地方曲率計(jì)算誤差較大,使得計(jì)算精度較差.IIF模型通過(guò)在粒子間施加相互作用力來(lái)模擬表面張力的影響,避免了計(jì)算界面曲率,對(duì)于交界面的尖角以及邊界處等粒子缺失和分布不均勻的地方張力計(jì)算穩(wěn)定高效.因此,本文構(gòu)造了一種新型的粒子間相互作用力模型來(lái)模擬表面張力的影響.

表面張力直接添加在動(dòng)量方程中

fij的表達(dá)式為

圖3 勢(shì)函數(shù)曲線圖

3 算例及分析

3.1 初始參數(shù)及計(jì)算模型

液滴沖擊固體表面的計(jì)算模型如圖4所示.與液滴沖擊固體表面問(wèn)題相關(guān)的物理量主要包括：液滴直徑D,液滴沖擊速度U,液體的密度ρ,黏性系數(shù)μ,以及表面張力系數(shù)σ,韋伯?dāng)?shù)We=ρU2D/σ,奧內(nèi)佐格數(shù)Oh=μ/(σρD)1/2.為便于分析,對(duì)以下參數(shù)進(jìn)行無(wú)量綱化處理：鋪展因子D?c=Dc/D,無(wú)量綱時(shí)間T?=UT/D.

為了準(zhǔn)確得到?jīng)_擊過(guò)程中液滴內(nèi)部及界面壓力的波動(dòng)情況,在數(shù)值模擬中,設(shè)置了三個(gè)監(jiān)測(cè)點(diǎn),監(jiān)測(cè)點(diǎn)1設(shè)置在液滴的圓心,監(jiān)測(cè)點(diǎn)2設(shè)置在液滴與固體表面的接觸點(diǎn),監(jiān)測(cè)點(diǎn)3設(shè)置在監(jiān)測(cè)點(diǎn)2正下方的固壁邊界處.

圖4 液滴沖擊固壁面模型圖

3.2 數(shù)值算例及結(jié)果分析

數(shù)值模擬中液滴的直徑D為4.2 mm,密度ρ為1000 kg/m3,黏性系數(shù)μ為0.001 N·s/m2.表面張力系數(shù)σ為0.0728 N/m,奧內(nèi)佐格數(shù)Oh為0.0018,韋伯?dāng)?shù)We為483.

圖5給出了液滴沖擊固壁面的液滴形態(tài)變化與壓力場(chǎng)演變過(guò)程.從圖5中我們可以清楚地觀察到液滴沖擊固壁面過(guò)程中液滴內(nèi)部壓力的整個(gè)波動(dòng)情況.隨著液滴開(kāi)始沖擊固壁面,接觸區(qū)域的壓力瞬間增大,這一點(diǎn)可以從監(jiān)測(cè)點(diǎn)的壓力曲線圖6中看出,并且初始韋伯?dāng)?shù)越大,壓力的峰值就越大.沖擊波一部分沿著固壁面向兩側(cè)傳遞,一部分在液滴內(nèi)部沿著沖擊方向反方向傳播,產(chǎn)生反方向的沖擊波.在接觸區(qū)域由于沖擊波作用,液滴沿固壁面向兩側(cè)逐步鋪展.同時(shí),由于相互作用粒子內(nèi)部能量逐漸耗散,液滴在慣性力,表面張力和重力的相互作用下運(yùn)動(dòng),直至達(dá)到平衡狀態(tài).

圖6給出了不同We數(shù)下的監(jiān)測(cè)點(diǎn)的壓力分布圖.從圖6可以看出,監(jiān)測(cè)點(diǎn)2即液滴與固壁面的接觸點(diǎn)的壓力在沖擊的瞬間增大,經(jīng)歷一次較大的壓力振蕩之后,壓力值急劇減小,隨后整個(gè)壓力場(chǎng)出現(xiàn)較小的波動(dòng).沖擊波傳播的整個(gè)過(guò)程中,在液滴內(nèi)部出現(xiàn)了壓力振蕩,這點(diǎn)從液滴圓心的檢測(cè)點(diǎn)曲線 p(1)可以看出,出現(xiàn)1次峰值,相比較曲線p(2)和 p(3),由于沖擊波反彈,壓力振蕩的幅度不大,約在0.4 ms的時(shí)候,沖擊波導(dǎo)致的壓力振蕩完全消失,之后液滴在慣性力的作用下運(yùn)動(dòng),直至趨于靜止.比較兩種不同韋伯?dāng)?shù)沖擊作用下的監(jiān)測(cè)點(diǎn)壓力分布可以得出：韋伯?dāng)?shù)越大,監(jiān)測(cè)點(diǎn)的壓力值就越大,能量轉(zhuǎn)化的幅度就越大.當(dāng)沖擊的速度較大或者韋伯?dāng)?shù)較高時(shí)(超過(guò)臨界韋伯?dāng)?shù)),粒子徑向運(yùn)動(dòng)的速度變大,在液滴與固壁面接觸區(qū)域兩側(cè)會(huì)產(chǎn)生斜向上的片狀射流,這部分流體所具有的能量足以克服表面張力的作用時(shí),在液滴與固壁面接觸區(qū)域兩側(cè)邊緣處會(huì)發(fā)生液滴的飛濺現(xiàn)象.

圖5 液滴形態(tài)變化與壓力場(chǎng)的演變過(guò)程(We=483)

圖6 不同We數(shù)下的監(jiān)測(cè)點(diǎn)的壓力分布圖(We=483,760)

圖7 鋪展因子隨時(shí)間的變化曲線

圖7 給出了液滴沖擊固壁面的鋪展因子D?c=Dc/D 隨無(wú)量綱時(shí)間t?=t(U/D)的變化曲線.實(shí)驗(yàn)觀察結(jié)果 f=2.8t?1/2由文獻(xiàn)[26]給出.從圖中可以得出：兩種不同韋伯?dāng)?shù)下的液滴沖擊固壁面鋪展階段,鋪展因子隨時(shí)間的變化曲線總體上比較接近實(shí)驗(yàn)觀察得出的結(jié)果;相同時(shí)刻,初始沖擊的韋伯?dāng)?shù)越大,鋪展因子就越大.采用改進(jìn)的SPH方法得到的結(jié)果與實(shí)驗(yàn)觀察得到的結(jié)果[26]基本一致.

4 結(jié)論

本文在傳統(tǒng)SPH方法的基礎(chǔ)上進(jìn)行了改進(jìn).為提高傳統(tǒng)SPH方法的精度和穩(wěn)定性,對(duì)密度和核梯度形式進(jìn)行了修正,保證粒子相互作用過(guò)程中總動(dòng)量始終守恒.采用了考慮黎曼解法的SPH流體控制方程,解決了接觸區(qū)域粒子材料性質(zhì)梯度變化所導(dǎo)致的間斷和不穩(wěn)定性.構(gòu)造了一種新型的IIF模型來(lái)模擬表面張力的影響,提高了交界面的尖角以及邊界處等粒子缺失和分布不均勻的地方曲率的計(jì)算精度.應(yīng)用改進(jìn)的SPH方法對(duì)液滴沖擊固壁面問(wèn)題進(jìn)行了數(shù)值模擬.計(jì)算結(jié)果表明：新型的IIF模型能夠較好地模擬表面張力的影響;改進(jìn)的SPH方法能夠精細(xì)地描述液滴與固壁面相互作用過(guò)程中液滴的內(nèi)部壓力場(chǎng)演變和自由面形態(tài)變化;得到了液滴的鋪展因子隨無(wú)量綱時(shí)間的變化曲線;相同時(shí)刻,韋伯?dāng)?shù)越大,鋪展因子就越大;數(shù)值模擬結(jié)果與實(shí)驗(yàn)觀察得到的結(jié)果基本一致.

[1]Tuan T 2012 Phys.Rev.Lett.108 036101

[2]Thoroddsen S T,Takehara K 2012 J.Fluid.Mech.706 560

[3]Zhang M K,Chen S,Shang Z 2012 Acta Phys.Sin.61 034701(in Chinese)[張明焜,陳碩,尚智2012物理學(xué)報(bào)61 034701]

[4]Liu M B,Chang J Z 2011 Int.J.Comput.Meth.8 637

[5]Worthington A M 1876 Proc.R.Soc.Lond.25 261

[6]Roisman I V,Opfer L,Tropea C,Raessi M,Mostaghimi J,Chandra S 2008 Colloids Surfaces A 322 183

[7]Qiang H F,Liu K,Chen F Z 2012 Acta Phys.Sin.61 204701(in Chinese)[強(qiáng)洪夫,劉開(kāi),陳福振2012物理學(xué)報(bào)61 204701]

[8]Bussmann M,Chandra S,Mostaghimi J 2000 Phys.Fluids 12 3121

[9]Eggers J,Fontelos M A,Josserand C,Zaleski S 2010 Phys.Fluids 22 301

[10]Ellis A S,Smith F T,White A H 2011 Q.J.Mech.Appl.Math.64 107

[11]Ma L Q,Chang J Z,Liu H T,Liu M B 2012 Acta Phys.Sin.61 054701(in Chinese)[馬理強(qiáng),常建忠,劉漢濤,劉謀斌2012物理學(xué)報(bào)61 054701]

[12]Yang B H,Wang H,Zhu X,Ding Y D,Zhou J 2012 CIESC J.10 3027(in Chinese)[楊寶海,王宏,朱恂,丁玉棟,周勁2012化工學(xué)報(bào)10 3027]

[13]Sikalo S,Wilhelm H D,Roisman I V,Jakirlic S,Tropea C 2005 Phys.Fluids 17 062103

[14]Liu M B,Liu G R,Zong Z,Lam K Y 2003 Comput.Fluids 32 305

[15]Liu M B,Liu G R,Zong Z 2008 Int.J.Comput.Meth.5 135

[16]Liu M B,Liu G R 2010 Arxiv.Comput.Methods Engrg.17 25

[17]Liu M B,Chang J Z 2010 Acta Phys.Sin.59 26

[18]Liu M B,Liu G R,Lam K Y 2003 Electron.Model.25 113

[19]Liu M B,Liu G R,Lam K Y 2003 J.Comput.Appl.Math.155 263

[20]Liu M B,Liu G R,Lam K Y,Zong Z 2003 Comput.Mech.30 106

[21]Liu M B,Shao J R 2012 Sci.China E 42 1

[22]Monaghan J J 1992 Annu.Rev.Astron.Astr.30 543

[23]Jiang T,Ouyang J,Zhao X K,Ren J L 2011 Acta Phys.Sin.60 054701(in Chinese)[蔣濤,歐陽(yáng)潔,趙曉凱,任金蓮2011物理學(xué)報(bào) 60 054701]

[24]Zhang S,Morita K,Fukuda K,Shirakawa N 2007 Int.J.Numer.Meth.Fl.55 225

[25]Liu M B,Liu G R 2005 Comput.Mech.35 332

[26]Rioboo R,Marengo M,Tropea C 2002 Exp.Fluids 33 112