最優(yōu)輸運(yùn)無網(wǎng)格方法及其在液滴表面張力效應(yīng)模擬中的應(yīng)用

周浩 李毅 劉海 陳鴻 任磊生

(中國空氣動(dòng)力研究與發(fā)展中心，超高速碰撞研究中心，綿陽 621000)

基于網(wǎng)格的數(shù)值方法(如有限元、有限體積、有限差分等)在大變形、不連續(xù)等問題中遇到挑戰(zhàn)，因此人們提出了多種無網(wǎng)格方法.最優(yōu)輸運(yùn)無網(wǎng)格方法是一種新型拉格朗日無網(wǎng)格方法，但是繼承了有限元方法在邊界表征、邊界處理等方面的優(yōu)勢，在表面張力模擬中具有較大潛力.基于拉格朗日方程，通過將表面張力勢能加入拉格朗日函數(shù)，得到的表面張力廣義力精確地作用在流體表面，而且表面張力系數(shù)是唯一的輸入?yún)?shù).給出了最優(yōu)輸運(yùn)無網(wǎng)格方法軸對(duì)稱離散格式.通過對(duì)二維/三維泊肅葉流、靜止和振動(dòng)的液滴、液滴變形等典型問題的仿真分析，驗(yàn)證了最優(yōu)輸運(yùn)無網(wǎng)格方法在表面張力問題模擬中的精度和收斂性.

1 引言

精確可靠的模擬表面張力對(duì)液滴變形與運(yùn)動(dòng)的影響在環(huán)境工程、微納米工程以及醫(yī)藥工程等領(lǐng)域有著十分重要的應(yīng)用.采用基于網(wǎng)格的數(shù)值方法模擬此類問題往往需要復(fù)雜的界面追蹤算法，典型的如流體體積函數(shù)法[1，2]和水平集方法[3]等.因此，越來越多的研究者采用拉格朗日粒子類方法模擬表面張力相關(guān)現(xiàn)象，利用粒子本身實(shí)現(xiàn)對(duì)邊界的自動(dòng)追蹤，極大地簡化了算法，其中最常用的數(shù)值方法是光滑粒子動(dòng)力學(xué)(smoothed particle hydrodynamics，SPH)[4，5]方法.在SPH 方法中，主要有兩種方法來考慮表面張力的影響:粒子間相互作用力(inter-particle interaction force，IIF)方法[6-10]和連續(xù)力方法(continuum surface force，CSF)[11].IIF 方法來源于如下思想，即表面張力來自于液體表面分子之間和內(nèi)部分子之間作用力的差異.受此啟發(fā)，在SPH 方法中，可以通過在宏觀粒子之間加入額外遠(yuǎn)程吸引力和近程排斥力來近似模擬表面張力.這種方法簡單靈活，粒子之間的額外作用力形式非常多.這種方法的缺點(diǎn)是表面張力系數(shù)不是輸入?yún)?shù)，而是需要先模擬某個(gè)標(biāo)準(zhǔn)算例，然后測量與這種宏觀粒子之間吸引力等價(jià)的表面張力系數(shù).此外，這種方法往往存在不收斂的問題.CSF方法通過給不同的流體賦以不同的色值(color)來區(qū)分不同的流體，利用不連續(xù)的色函數(shù)(color function)求解表面法向以及曲率來計(jì)表面張力，然后將只作用在表面的力轉(zhuǎn)換為連續(xù)的體積力.SPH方法中，界面法向向量及其曲率的計(jì)算誤差一般較大.無論哪種表面張力模擬方法，SPH 方法都存在如下固有缺點(diǎn):1) SPH 形函數(shù)不能嚴(yán)格滿足1 階精度(特別是粒子分布不均勻或者自由界面處)，導(dǎo)致收斂性不好，常常需要對(duì)形函數(shù)及其導(dǎo)數(shù)進(jìn)行修正[12，13];2)形函數(shù)不滿足Kronecker delta 性質(zhì)，因此位移邊界條件處理復(fù)雜，一般需要通過各種虛粒子處理邊界條件.因此，本文嘗試采用一種新型無網(wǎng)格方法模擬表面張力問題，即最優(yōu)輸運(yùn)無網(wǎng)格 (optimized transportation meshfree，OTM)方法[14].

OTM 方法采用粒子系統(tǒng)離散連續(xù)介質(zhì)，利用局部最大熵(local maximum entropy，LME)[15-17]形函數(shù)構(gòu)造連續(xù)場，并結(jié)合數(shù)學(xué)中的最優(yōu)輸運(yùn)理論計(jì)算系統(tǒng)的動(dòng)能積分.相比于其他無網(wǎng)格方法，OTM 方法在精度、收斂性、邊界處理以及數(shù)學(xué)理論基礎(chǔ)上均具有優(yōu)勢.OTM 方法具有無網(wǎng)格方法能夠處理大變形的優(yōu)點(diǎn)，同時(shí)又繼承了有限元方法的部分優(yōu)點(diǎn)，如形函數(shù)嚴(yán)格滿足一階精度、形函數(shù)在邊界處滿足弱Kronecker delta 性質(zhì)等.當(dāng)前，OTM 方法已經(jīng)成功應(yīng)用于金屬侵徹[18]、超高速碰撞[19-21]、裂紋擴(kuò)展[22]以及熱傳導(dǎo)等問題[23]，但是將其應(yīng)用于表面張力模擬還未見文獻(xiàn)報(bào)道.

本文從拉格朗日方程出發(fā)，從能量的角度考慮表面張力的影響，即將表面能加入到拉格朗日函數(shù)中，得到的表面張力廣義力精確地作用在流體表面，而且表面張力系數(shù)是唯一的輸入?yún)?shù).為了減少計(jì)算量，推導(dǎo)了軸對(duì)稱形式的OTM 離散格式，用于三維算例數(shù)值模擬，并給出了詳細(xì)的計(jì)算流程.通過模擬泊肅葉流并與解析解對(duì)比，驗(yàn)證了本文基本離散格式和軸對(duì)稱處理的正確性.通過模擬靜止液滴，并與Young-Laplace 公式對(duì)比，驗(yàn)證了表面張力處理的正確性.模擬了液滴的周期振蕩問題并與解析解對(duì)比，嚴(yán)格驗(yàn)證了液滴中速度分布的空間收斂性.最后，模擬了液滴在重力、表面張力以及界面共同作用下的形變問題，考察了液滴中的壓力分布并與理論解對(duì)比，進(jìn)一步驗(yàn)證了本文OTM 離散格式.

2 OTM 方法

2.1 連續(xù)介質(zhì)離散

OTM 方法中的連續(xù)介質(zhì)離散介于拉格朗日有限元與SPH 方法之間.OTM 方法將有限元網(wǎng)格中每個(gè)小網(wǎng)格轉(zhuǎn)換為1 個(gè)粒子，位于小網(wǎng)格形心.對(duì)于三角形，形心為三條中線的交點(diǎn).粒子存儲(chǔ)所有的物理量，與SPH 方法類似.同時(shí)，保留所有的網(wǎng)格節(jié)點(diǎn)，其中存儲(chǔ)了速度和位置，如圖1 所示.初始時(shí)刻，邊界上只有節(jié)點(diǎn)，且節(jié)點(diǎn)隨著流場運(yùn)動(dòng).邊界節(jié)點(diǎn)很好地表征了連續(xù)介質(zhì)的邊界位置，與拉格朗日有限元方法類似.此外，由于局部最大熵形函數(shù)滿足弱Kronecker delta 性質(zhì)，邊界條件處理簡單.拋棄有限元方法中單元和節(jié)點(diǎn)之間的固定連接關(guān)系，而是采用近鄰粒子搜索方式動(dòng)態(tài)確定，從而可以處理大變形，這與SPH 方法類似.

圖1 利用有限元網(wǎng)格將連續(xù)介質(zhì)離散為粒子和節(jié)點(diǎn)(a)有限元網(wǎng)格;(b)用粒子與節(jié)點(diǎn)離散連續(xù)介質(zhì)Fig.1.Continuum discretized by particles (red symbols)and nodes (white symbols) in the OTM method by means of finite element mesh:(a) Finite element mesh;(b) continuum discretized by particles and nodes.

2.2 LME 形函數(shù)

最優(yōu)輸運(yùn)無網(wǎng)格方法與更新拉格朗日有限元方法最大的1 個(gè)區(qū)別是采用LME 形函數(shù)，因此不需要網(wǎng)格，屬于無網(wǎng)格方法.給定1 個(gè)節(jié)點(diǎn)集xa(a1，2，···，n)，每個(gè)節(jié)點(diǎn)的形函數(shù)記為Na(x) .利用形函數(shù)可以從離散數(shù)值插值出1 個(gè)連續(xù)場:

其中，uh(x)是插值出來的連續(xù)場，ua是物理量u在節(jié)點(diǎn)a處的值.在有限元和SPH 等方法中，形函數(shù)及其導(dǎo)數(shù)都具有簡單解析表達(dá)式，但是LME形函數(shù)沒有顯式的解析表達(dá)式，而是需要數(shù)值求解.

其中，由于x是固定值所以將其省略.為了減少計(jì)算量，形函數(shù)一般都是局域的，即當(dāng) |x-xa| 大于某一臨界值時(shí)，形函數(shù)Na(x) 為零.因此，可以定義如下形函數(shù)影響域平均大小:

形函數(shù)影響域越大，計(jì)算量越大.LME 形函數(shù)的基本思路是讓形函數(shù)熵最大，同時(shí)影響域又盡量小，因此，構(gòu)造了如下優(yōu)化問題:

其中，β是1 個(gè)可調(diào)參數(shù)，用于調(diào)節(jié)兩個(gè)目標(biāo)的權(quán)重.最后1 個(gè)約束條件表示形函數(shù)嚴(yán)格滿足一階精度.上述優(yōu)化問題解的存在性、唯一性以及求解方法文獻(xiàn)[15-17]中有詳細(xì)論述，本文只給出最后結(jié)果.LME 形函數(shù)可以寫成如下形式:

其中Z函數(shù)的定義為

λ*(x)是以下無約束優(yōu)化問題的解:

定義:

那么，可以采用如下迭代方法求解λ*(x) :

初始時(shí)刻，λ(0)一般取零.流場發(fā)生變化后，λ(0)一般取上1 個(gè)時(shí)刻的λ*.當(dāng)β等于零時(shí)，得到的最大熵形函數(shù)是全局的，不適合數(shù)值計(jì)算.當(dāng)β趨于無限大時(shí)，即不考慮熵最大，只要求形函數(shù)影響域最小，那么LME 形函數(shù)可以連續(xù)過渡到有限元線性形函數(shù).LME 形函數(shù)具有以下特點(diǎn):1)滿足一階精度;2)不依賴網(wǎng)格;3)邊界處滿足弱Kronecker delta 性質(zhì);4)形函數(shù)非負(fù).

得到λ*(x) 與形函數(shù)后，形函數(shù)導(dǎo)數(shù)為

2.3 基本離散方程

首先考慮無黏、可壓縮以及等溫流體.當(dāng)初始狀態(tài)已知，節(jié)點(diǎn)的位置描述了流場的形變的流動(dòng)，實(shí)際上完全描述了流場，因此，節(jié)點(diǎn)的位置可以當(dāng)成系統(tǒng)的廣義坐標(biāo).本文下標(biāo)p表示粒子，用下標(biāo)a，b，c表示節(jié)點(diǎn).由于無黏，這是1 個(gè)保守系統(tǒng)，拉格朗日方程為

其中，E是系統(tǒng)動(dòng)能，V 是系統(tǒng)勢能，fa是廣義節(jié)點(diǎn)力，va是節(jié)點(diǎn)速度，xa是節(jié)點(diǎn)位置，N是節(jié)點(diǎn)總數(shù).系統(tǒng)動(dòng)能和勢能為所有粒子動(dòng)能和內(nèi)能之和:

其中mp表示粒子質(zhì)量，vp表示粒子速度，ep表示粒子單位質(zhì)量內(nèi)能，ρp表示粒子密度.

粒子內(nèi)能和壓力的關(guān)系通過物態(tài)方程(equation of state，EOS)描述:

與有限元一樣，粒子位置和速度可以通過節(jié)點(diǎn)位置和速度插值得到:

其中 δxp和 δxa分別為粒子和節(jié)點(diǎn)的虛位移，Na(xp)表示節(jié)點(diǎn)a處形函數(shù)在粒子p處的值，?Na(xp) 表示形函數(shù)梯度，求和是對(duì)粒子的所有近鄰節(jié)點(diǎn).粒子密度變化和位置變化的關(guān)系由連續(xù)性條件給出:

將(14e)式代入(15)式，得到

利用求導(dǎo)鏈?zhǔn)揭?guī)則，并且利用(13)式和(16)式，得到

利用(12a)式、(14c)式和(14d)式得到

在(18)式中交換求和順序，拉格朗日方程變成:

其中，M為相容質(zhì)量矩陣

(19)式與更新拉格朗日有限元格式有兩點(diǎn)區(qū)別:1)有限元格式中的單元積分對(duì)應(yīng)了OTM 中的粒子積分;2)有限元中使用基于網(wǎng)格的形函數(shù)，而OTM 方法中使用無網(wǎng)格的局部最大熵形函數(shù).注意到局部最大熵形函數(shù)可以連續(xù)過渡到有限元線性形函數(shù)，因此，OTM 方法可以看成是一種特殊的單點(diǎn)積分拉格朗日有限元方法.

在有限元中方法，為避免矩陣求逆，常常將相容質(zhì)量矩陣按行求和，得到集中質(zhì)量矩陣，大大減少了計(jì)算量和存儲(chǔ)量，OTM 方法中同樣如此.利用形函數(shù)的單位分解性質(zhì)，可以將粒子的質(zhì)量分配給其相鄰的節(jié)點(diǎn)，得到如下節(jié)點(diǎn)質(zhì)量定義:

注意到粒子質(zhì)量和局部最大熵形函數(shù)都是嚴(yán)格大于零的，因此，節(jié)點(diǎn)質(zhì)量也是嚴(yán)格大于零的.這樣連續(xù)介質(zhì)既可以看成是粒子系統(tǒng)，也可以看成是節(jié)點(diǎn)系統(tǒng).很明顯，這兩種系統(tǒng)的總質(zhì)量和總動(dòng)量都是嚴(yán)格相等的.

系統(tǒng)的總勢能通過粒子系統(tǒng)計(jì)算.系統(tǒng)的總動(dòng)能既可以通過粒子系統(tǒng)計(jì)算，也可以通過節(jié)點(diǎn)系統(tǒng)近似計(jì)算:

將(23)式代入(11)式，得到如下離散方程:

從(20)式和(21)式可知，節(jié)點(diǎn)質(zhì)量剛好就是相容質(zhì)量矩陣按行求和.利用節(jié)點(diǎn)廣義力和節(jié)點(diǎn)質(zhì)量可以更新節(jié)點(diǎn)的速度和位置:

用(25b)式更新節(jié)點(diǎn)位置后，粒子位置根據(jù)(14a)式更新.注意到節(jié)點(diǎn)位置變化代表了流場的形變，因此兩個(gè)相鄰時(shí)刻的連續(xù)映射φ:xk →xk+1可以通過更新的節(jié)點(diǎn)位置插值得到

上述映射的梯度即為變形梯度張量:

變形梯度張量的行列式在粒子位置處的值表征了粒子的密度變化

(28)式可以用于更新粒子密度.粒子密度更新后，根據(jù)物態(tài)方程更新粒子壓力.

2.3.1 流體黏性的處理

以上推導(dǎo)沒有考慮黏性，因而是保守系統(tǒng).處理非保守力的基本想法將其看成外力，然后利用虛功原理推導(dǎo)對(duì)應(yīng)的廣義力.假設(shè)fp是作用在粒子p上的外力，那么此力的虛功為

由虛位移的獨(dú)立性，得到相應(yīng)的廣義力為

例如，作用在1個(gè)粒子上的重力fpmpg.根據(jù)(30)式得到廣義力考慮到重力實(shí)際上是保守力，因此也可以將重力勢能加入拉格朗日函數(shù)中，得到重力對(duì)應(yīng)的廣義力.簡單計(jì)算可以證明兩者是一致的.

材料的很多行為(如彈塑性，黏彈性等)都可以統(tǒng)一用1 個(gè)應(yīng)力張量描述，而應(yīng)力張量的影響也可以看成外力，如下式中的動(dòng)量守恒方程:

應(yīng)力張量散度的物理意義為單位體積受到的力，因此，應(yīng)力張量虛功為

其中用到了高斯積分公式，并且假設(shè)了邊界上虛位移或者應(yīng)力張量為零，這在固壁邊界以及自由表面邊界上都成立.同理，得到廣義力為

如果只考慮牛頓流體的黏性，那么黏性應(yīng)力張量為

其中μ為黏性系數(shù).

同理，(31)式也能通過勢能和拉格朗日方程得到.從(33)式可以看到，1 個(gè)粒子在無限短時(shí)間內(nèi)受到的力只和應(yīng)力張量本身有關(guān)，而和應(yīng)力應(yīng)變關(guān)系無關(guān).因此，可以臨時(shí)假設(shè)1 個(gè)線性應(yīng)力應(yīng)變關(guān)系，從而得到簡單的彈性勢能.對(duì)彈性勢能求變分，就得到了相應(yīng)的廣義力，詳細(xì)證明見附錄A.兩種方法的結(jié)果是一致的.

2.3.2 表面張力處理

考慮二維情況，液體邊界被邊界節(jié)點(diǎn)表征，如圖2 所示.

圖2 二維條件下的邊界節(jié)點(diǎn)(白點(diǎn))Fig.2.Boundary nodes (white symbols) in two-dimensional coordinate.

表面勢能正比于表面積:

其中rab|xa-xb| ，rac|xa-xc| 為節(jié)點(diǎn)之間的間距.γ是表面張力系數(shù).對(duì)(35)式求變分，得到廣義力

可以看到，(36)式非常簡潔，可以解釋為節(jié)點(diǎn)a同時(shí)受到兩個(gè)相鄰節(jié)點(diǎn)b和c的吸引力，吸引力大小剛好等于表面張力系數(shù).(36)式中沒有可調(diào)參數(shù)，表面張力系數(shù)是唯一的輸入?yún)?shù).此外，表面張力對(duì)應(yīng)的廣義力只作用在表面節(jié)點(diǎn)上，即精確的作用在流體表面.

2.3.3 軸對(duì)稱處理

記(r，θ，h)為柱坐標(biāo)系，并且假設(shè)環(huán)向沒有位移，而且所有物理量都和環(huán)向坐標(biāo)無關(guān)，即 δuθ0，?/?θ0.定義:

那么，(11)—(14)式依然成立.速度矢量散度在柱坐標(biāo)下的表達(dá)式為

因此，粒子密度變化和粒子位置變化關(guān)系式(15)需要修改為

于是

其中nr是徑向單位向量.于是，軸對(duì)稱條件下的廣義力為

(18)—(26)式仍然成立.注意(27)式中的變形梯度張量的行列式表征的是粒子在(r，h)平面上的面積的變化:

因此，可以先更新粒子在(r，h)平面上的面積Sk+1，然后再更新密度:

很明顯，外力對(duì)應(yīng)的廣義力(30)式仍然成立.黏性力對(duì)應(yīng)的廣義力也仍然成立，詳細(xì)證明見附錄B.

三維軸對(duì)稱條件下，表面勢能為

式中，前兩項(xiàng)為節(jié)點(diǎn)b和節(jié)點(diǎn)c對(duì)節(jié)點(diǎn)a的吸引力，后兩項(xiàng)為對(duì)稱軸對(duì)邊界節(jié)點(diǎn)的吸引力.可見，笛卡爾坐標(biāo)系下的OTM 計(jì)算格式只需要經(jīng)過微小改動(dòng)，就能計(jì)算三維軸對(duì)稱問題.

2.4 OTM 方法詳細(xì)步驟

笛卡爾坐標(biāo)系下，采用OTM 方法模擬表面張力的詳細(xì)步驟總結(jié)如下.

1)初始化:給定節(jié)點(diǎn)的位置和速度，給定粒子的位置和其他物理量.

2)近鄰粒子搜索:找到所有粒子-節(jié)點(diǎn)對(duì).

3)對(duì)于所有粒子-節(jié)點(diǎn)對(duì)，計(jì)算形函數(shù)Na，p和形函數(shù)導(dǎo)數(shù)?Na，p.

4)計(jì)算節(jié)點(diǎn)質(zhì)量、粒子速度梯度、黏性應(yīng)力張量和廣義節(jié)點(diǎn)力:

5)更新節(jié)點(diǎn)速度和位置:

6)根據(jù)邊界條件調(diào)整節(jié)點(diǎn)位置，然后更新粒子位置:

7)計(jì)算變形梯度張量并更新粒子密度和壓力:

8)完成1 個(gè)時(shí)間步長，回到步驟2).

3 典型算例

所有的初始化借助于有限元網(wǎng)格.網(wǎng)格節(jié)點(diǎn)變成節(jié)點(diǎn)，網(wǎng)格形心放置粒子，粒子體積即為網(wǎng)格大小.物態(tài)方程如下:

其中，ρ0表示材料初始密度，ρ表示當(dāng)前密度，c0表示聲速，p表示壓力.

3.1 泊肅葉流

兩塊無限長的二維平板中間充滿了靜止液體，液體在平行于平板方向的體積力作用下開始運(yùn)動(dòng)，稱之為泊肅葉(Poiseuille)流.如果是無限長的三維圓管，那么稱之為Hagen-Poiseuille 流.對(duì)于二維泊肅葉，參數(shù)和文獻(xiàn)[24]一樣，解析解見文獻(xiàn)[24].將平板距離變成圓管半徑，即為Hagen-Poiseuille流，解析解見文獻(xiàn)[25].OTM 模擬結(jié)果見圖3.

圖3 速度分布的OTM 模擬與解析解對(duì)比 (a) 二維泊肅葉流;(b) 三維軸對(duì)稱泊肅葉流Fig.3.Comparison of OTM (symbols) and analytic solutions (solid curves) for velocity profile:(a) Two-dimensional Poiseuilleflow;(b) axisymmetric Hagen-Poiseuille flow.

OTM 模擬結(jié)果與理論解吻合很好.為了模擬無限長管道，用到了周期邊界條件.

3.2 靜止液滴

二維液滴的半徑R=0.2 m，密度ρ=1 kg/m3，表面張力系數(shù)γ=1 Pa·m，黏性系數(shù)η=0.5 Pa·s，聲速c0=50 m/s.靜止?fàn)顟B(tài)下液滴的理論壓力值可以根據(jù)Young-Laplace 公式得到ρ=γ/R=5 Pa.如果是三維液滴，那么其理論壓力值為ρ=2γ/R=10 Pa.OTM 模擬的平均壓力分別為5.004 Pa 和 9.967 Pa.壓力分布如圖4 所示.

圖4 靜止液滴的壓力場 (a) 二維液滴;(b) 軸對(duì)稱三維液滴Fig.4.Pressure field:(a) Two-dimensional rod;(b) axisymmetric three dimensional drop.

由于弱可壓假設(shè)，壓力場中有誤差.但是液滴中的平均壓力非常接近理論值(誤差0.5%以內(nèi))，證明了本方法的精度.

3.3 液滴的周期振蕩

在上1 個(gè)算例的基礎(chǔ)上，將黏性減小到η5×10-3Pa·s，并且附加1 個(gè)散度為0 的初始速度場vxx，vy-y給所有節(jié)點(diǎn).二維液滴的理論振蕩周期為T[26].為了驗(yàn)證OTM方法的收斂性，采用不同的粒子總數(shù)模擬這個(gè)問題，并且檢查兩條相互垂直線x0和y0 上的物理量分布，時(shí)刻tT/2，如圖5 所示.

圖5 t=T/2 時(shí)刻的壓力和速度分布 (a) x=0 上的壓力分布;(b) x=0 上的速度分布;(c) y=0 上的壓力分布;(d) y=0 上的速度分布Fig.5.Pressure and velocity profile at t=T/2:(a) Pressure profile at x=0;(b) velocity profile at x=0;(c) pressure profile at y=0;(d) velocity profile at y=0.

可以看到，壓力分布的收斂不明顯，但是速度分布的收斂很明顯.改變表面張力系數(shù)的大小，得到的周期與理論值的對(duì)比如圖6 所示.周期是根據(jù)液滴右上部分質(zhì)心軌跡測量得到，與文獻(xiàn)[26]一樣.

圖6 二維液滴振蕩 (a) 振蕩周期理論解和數(shù)值解得對(duì)比;(b) 表面張力系數(shù)為 γ=1 時(shí)液滴右上部分質(zhì)心的軌跡Fig.6.Two-dimensional rod oscillation:(a) Comparison of period between the numerical (symbols) and analytical (solid curve) results;(b) center of mass position history when γ=1 .

對(duì)于軸對(duì)稱的三維液滴，散度為零的初始速度場取為vrr，vh-2h，其他物理量不變.振蕩周期的理論解為T[27].模擬結(jié)果與理論的對(duì)比如圖7 所示.模擬的周期根據(jù)對(duì)稱軸端點(diǎn)的位置軌跡測量.可以看到，兩種情況下OTM計(jì)算結(jié)果和理論非常吻合.

圖7 軸對(duì)稱三維液滴振蕩周期與理論的對(duì)比Fig.7.Three-dimensional droplet oscillation，comparing of period between the numerical (symbols) and analytical (solid curve) results.

3.4 液滴形變

將液滴半徑擴(kuò)大到 1 m，表面張力系數(shù)擴(kuò)大到10 Pa·m，加上大小為10 m/s2的重力.y=—1 m 處放置1 個(gè)光滑的固定平板.液滴在重力、表面張力以及光滑平板的共同作用下產(chǎn)生變形，如圖8 所示.

在圖8(d)中，壓力等值線與重力垂直，和理論相符.液滴內(nèi)部的壓力場需要滿足兩個(gè)條件，一是壓力沿Y軸的線性分布規(guī)律 δpρgδY，另一個(gè)是自由表面附近的壓力滿足Young-Laplace 方程pγ/R.這兩個(gè)條件實(shí)際上決定了液滴的形狀，即液滴自由表面部分的微分方程為

圖8 壓力場 (a) t=0 s;(b) t=0.02 s;(c) t=0.4 s;(d) t=4 sFig.8.Pressure field at several typical times:(a) t=0 s;(b) t=0.02 s;(c) t=0.4 s;(d) t=4 s.

其中H和p0分別為液滴最高點(diǎn)的Y坐標(biāo)和壓力.上述理論曲率半徑在實(shí)際計(jì)算時(shí)誤差太大，因此，本文采用如下方法計(jì)算曲率半徑:對(duì)于每個(gè)邊界節(jié)點(diǎn)，找到與其相鄰的4 個(gè)表面節(jié)點(diǎn)，然后用1 個(gè)圓來擬合這5 個(gè)節(jié)點(diǎn)，即自由變量為圓的中心坐標(biāo)和半徑.圓的半徑即為此節(jié)點(diǎn)的曲率半徑.

液體中自由邊界的右邊部分如圖9(a)所示.采用曲率半徑計(jì)算的邊界壓力和根據(jù)EOS 計(jì)算的壓力比較如圖9(b)所示.

圖9 邊界位置和壓力 (a) 邊界位置;(b) 邊界壓力與Y 軸坐標(biāo)的關(guān)系，壓力根據(jù)表面曲率和EOS 計(jì)算Fig.9.Boundary position and pressure profile:(a) Boundary position;(b) boundary pressure versus Y coordinate，computed by Young-Laplace equation (circles) and EOS (points).

兩種方法得到壓力值大致相等，說明模擬的自由表面形狀和理論相符.靠近固壁處誤差較大，經(jīng)檢查是因?yàn)楣瘫诟浇?jié)點(diǎn)的分布更加無序，導(dǎo)致搜索算法對(duì)初始參數(shù)敏感.壓力在Y軸方向的線性分布和理論相符，擬合直線的斜率也和理論dp/dY-ρg-10Pa/m 相符.

4 結(jié)論

本文基于拉格朗日方程，通過將表面張力勢能加入拉格朗日函數(shù)，得到了OTM 方法框架下的表面張力效應(yīng)表達(dá)式，并且給出了軸對(duì)稱的處理方式.模擬了泊肅葉流，靜止和振蕩的液滴，液滴在重力、表面張力以及光滑平板共同作用下的形變等典型算例，通過與解析解對(duì)比，驗(yàn)證了OTM 方法在模擬表面張力現(xiàn)象中具有如下優(yōu)勢:1)表面邊界節(jié)點(diǎn)精確的表征了邊界的位置，保證了表面張力勢能的準(zhǔn)確計(jì)算，而且表面張力對(duì)應(yīng)的廣義力精確的作用在流體表面;2)表面張力對(duì)應(yīng)的廣義力形式簡潔，具有直觀的物理意義，而且表面張力系數(shù)是唯一的輸入?yún)?shù)，不需要任何可調(diào)參數(shù);3)可以保證速度分布的空間收斂性.

本文在連續(xù)介質(zhì)離散時(shí)，節(jié)點(diǎn)和節(jié)點(diǎn)、節(jié)點(diǎn)和粒子之間沒有任何固定聯(lián)系，而是通過近鄰粒子搜索動(dòng)態(tài)確定.但是，在計(jì)算表面張力勢能時(shí)，為了方便計(jì)算表面積，假設(shè)了表面節(jié)點(diǎn)之間的連接關(guān)系不變.因此，所有表面節(jié)點(diǎn)可以看成是1 個(gè)表面有限元網(wǎng)格.表面有限元網(wǎng)格與內(nèi)部節(jié)點(diǎn)在拉格朗日方程框架下相互作用和協(xié)調(diào)運(yùn)動(dòng).下一步，將考慮表面拓?fù)浣Y(jié)構(gòu)發(fā)生變化的情況，如多個(gè)液滴的融合過程.

附錄A 黏性應(yīng)力張量對(duì)應(yīng)廣義力的勢能方法推導(dǎo)

下面給出(33)式的勢能推導(dǎo)方法.將應(yīng)力和應(yīng)變張量寫成向量形式:

應(yīng)變與位移的關(guān)系為

其中L是微分矩陣

u是位移.廣義胡克定律為

其中c是1 個(gè)對(duì)稱矩陣:

以下4 個(gè)等式成立:

將(16)式寫成矩陣形式:

其中約定了形函數(shù)導(dǎo)數(shù)?Na，p為列向量.對(duì)于固定的應(yīng)力張量σp，由于c可以取得足夠大，使得應(yīng)變?yōu)樾∽冃?，于是彈性勢能?/p>

對(duì)彈性勢能求變分，得到

將向量形式改寫為常用的張量形式，得到

(A12)式第二項(xiàng)中有應(yīng)變，因此是小量，可以忽略

(A13)式和 (33)式一致.

附錄B 軸對(duì)稱條件下黏性力對(duì)應(yīng)的廣義力

兩個(gè)張量的雙點(diǎn)乘是1 個(gè)標(biāo)量，因此，可以在柱坐標(biāo)下求解

考慮到軸對(duì)稱，得到

利用定義(37a)式和(37c)式，虛功為

(B3)式與(32)式相同.