基于光滑粒子動力學方法的液滴沖擊固壁面問題數值模擬*

蘇鐵熊 馬理強? 劉謀斌 常建忠

1 引言

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

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

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

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

2 光滑粒子動力學

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

SPH方法使用一組任意分布的粒子來代表所模擬的介質系統,基于粒子體系離散控制流體運動的偏微分方程,并通過適當的變換對所獲得的SPH運動方程進行求解,從而仿真流體系統的動力學特征.利用SPH方法對偏微分方程的近似包括兩步：核近似和粒子近似.對于傳統的SPH方法,在某點(或某個粒子)i上,對任意函數 f(x)的核近似(〈f(x)〉)可由下式定義

其中,f為三維坐標向量x的函數,Ω為包含x的積分區域,W是光滑函數,h定義了光滑函數W的影響區域,稱為光滑長度.

f(x)的離散形式的粒子近似是對相關粒子i支持域內所有粒子進行加權求和(如圖1)

式中N為粒子i的支持域內所有粒子總數,ρj和mj為粒子 j的密度和質量,Wij為粒子 j對粒子i產生影響的光滑函數.(2)式表明粒子i處的任一函數值可通過應用光滑函數對其緊支域內所有粒子相對應的函數值進行加權平均進行近似.

圖1 二維空間SPH粒子近似示意圖 W為光滑函數,支持域為κh,S為計算區域Ω的表面

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

液滴沖擊固壁面過程中,液滴和固壁面接觸區域粒子的材料性質(如密度、壓強、速度等)梯度變化很大,會導致局部的非連續和間斷,利用傳統SPH的流體控制方程算法精度往往不高,在較大的We數沖擊作用下可能會導致計算結果的失真和計算中止.基于黎曼解法對于非連續問題和間斷問題求解的有效性,本文采用考慮黎曼解法的SPH流體控制方程：

2.3 密度和核梯度修正

在早期的文獻中,就近似格式而言,SPH通常被認為是二階精度的方法[22].對(1)式右邊在x處進行泰勒級數展開

其中r為余項.因為光滑函數滿足正則化條件和對稱性條件,(5)式可寫為

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

然而,(6)式未必對所有的情形都成立.例如,當某個粒子臨近計算區域邊界時,粒子的支持域與計算區域相交,支持域被邊界截斷(見圖2),正則化條件和對稱性條件都不能滿足,因此函數的SPH核近似不具有二階精度.

圖2 粒子的支持域與計算區域相交

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

為提高傳統SPH方法的計算精度,參照文獻[23]對SPH方法中的密度和核梯度進行了修正

對(1)式等號右邊項基于泰勒級數展開,可以得到

由于

可以得到一個新的核梯度表達式

粒子近似方程為

最終,核梯度修正模型為

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

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

表面張力直接添加在動量方程中

fij的表達式為

圖3 勢函數曲線圖

3 算例及分析

3.1 初始參數及計算模型

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

為了準確得到沖擊過程中液滴內部及界面壓力的波動情況,在數值模擬中,設置了三個監測點,監測點1設置在液滴的圓心,監測點2設置在液滴與固體表面的接觸點,監測點3設置在監測點2正下方的固壁邊界處.

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

3.2 數值算例及結果分析

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

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

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

圖5 液滴形態變化與壓力場的演變過程(We=483)

圖6 不同We數下的監測點的壓力分布圖(We=483,760)

圖7 鋪展因子隨時間的變化曲線

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

4 結論

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

[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物理學報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)[強洪夫,劉開,陳福振2012物理學報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)[馬理強,常建忠,劉漢濤,劉謀斌2012物理學報61 054701]

[12]Yang B H,Wang H,Zhu X,Ding Y D,Zhou J 2012 CIESC J.10 3027(in Chinese)[楊寶海,王宏,朱恂,丁玉棟,周勁2012化工學報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)[蔣濤,歐陽潔,趙曉凱,任金蓮2011物理學報 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