林晨森 陳碩? 肖蘭蘭

1)(同濟(jì)大學(xué)航空航天與力學(xué)學(xué)院,上海 200092)

2)(上海工程技術(shù)大學(xué)機(jī)械與汽車(chē)工程學(xué)院,上海 201620)

耗散粒子動(dòng)力學(xué)(DPD)是一種針對(duì)介觀流體的高效的粒子模擬方法,經(jīng)過(guò)二十多年發(fā)展已經(jīng)在諸如聚合物、紅細(xì)胞、液滴浸潤(rùn)性等方面有了很多研究應(yīng)用.但是因?yàn)槠溥吔缣幚硎侄蔚牟煌晟?耗散粒子動(dòng)力學(xué)模擬仍局限于相對(duì)簡(jiǎn)單的幾何邊界問(wèn)題中.本文提出一種能自適應(yīng)各種復(fù)雜幾何邊界的處理方法,并能同時(shí)滿足三大邊界要求：流體粒子不穿透壁面、邊界處速度無(wú)滑移、邊界處密度和溫度波動(dòng)小.具體地,通過(guò)給每個(gè)壁面粒子賦予一個(gè)新的矢量屬性—局部壁面法向量,該屬性通過(guò)加權(quán)計(jì)算周?chē)诿媪Ｗ拥奈恢玫玫剑蝗缓笸ㄟ^(guò)定義周?chē)腆w占比概念,僅提取固體壁面的表層粒子參與模擬計(jì)算,減少了模擬中無(wú)效的粒子；最后在運(yùn)行中,實(shí)時(shí)計(jì)算每個(gè)流體粒子周?chē)腆w粒子占比,判斷是否進(jìn)入固體壁面內(nèi),如果進(jìn)入則修正速度和位置.我們將這種方法應(yīng)用于Poiseuille流動(dòng),驗(yàn)證了該方法符合各項(xiàng)要求,隨后還在復(fù)雜血管網(wǎng)絡(luò)和結(jié)構(gòu)化固體壁面上展示了該邊界處理方法的應(yīng)用.這種方法使得DPD模擬不再局限于簡(jiǎn)單函數(shù)描述的壁面曲線,而是可以直接從各種設(shè)計(jì)圖紙和實(shí)驗(yàn)掃描影像中提取壁面,極大地拓展了DPD的應(yīng)用范圍.

1 引 言

耗散粒子動(dòng)力學(xué)(DPD)由Hoogerbrugge和Koelman[1]在 1992年提出構(gòu)想,隨后由 Espa?ol和Warren[2]做出重要的系統(tǒng)論證并搭建了理論框架,經(jīng)過(guò)二十幾年的發(fā)展,至今已經(jīng)成為介觀尺度模擬流體動(dòng)力學(xué)的最主要的方法之一,被稱為“彌補(bǔ)宏觀尺度和介觀尺度之間空白的方法”[3].它的基本思想是一些離散的被稱為粒子的動(dòng)量載體在連續(xù)的空間和離散的時(shí)間上運(yùn)動(dòng),這些粒子代表的是一個(gè)小區(qū)域內(nèi)的大量分子的集體行為,這樣就可以在粗粒化粒子尺度上進(jìn)行研究,從而忽略更小尺度的結(jié)構(gòu)和運(yùn)動(dòng)狀況.經(jīng)過(guò)人們的不斷探索和嘗試,DPD在以下應(yīng)用場(chǎng)景中展現(xiàn)出獨(dú)特的模擬優(yōu)勢(shì): 1)聚合物,采用珠簧鏈模型,可以建立和Flory-Huggins理論之間的對(duì)應(yīng)關(guān)系,研究微相分離體系和自組裝現(xiàn)象,可用于藥物輸運(yùn)、介孔材料制備等研究; 2)紅細(xì)胞,DPD紅細(xì)胞模型可以準(zhǔn)確對(duì)應(yīng)真實(shí)紅細(xì)胞的力學(xué)特性,已經(jīng)再現(xiàn)了單個(gè)紅細(xì)胞穿過(guò)狹縫、多個(gè)紅細(xì)胞聚集成串、紅細(xì)胞在血管狹窄處堵塞流動(dòng)等現(xiàn)象[4,5]; 3)液滴浸潤(rùn),對(duì)DPD勢(shì)能稍加修改就可以模擬氣液共存的系統(tǒng)[6],對(duì)液滴在各種化學(xué)和幾何異質(zhì)表面上運(yùn)動(dòng)行為的研究為預(yù)防表面結(jié)冰和微流控芯片中的液滴控制帶來(lái)新的啟示[7,8].此外,隨著近年來(lái)計(jì)算機(jī)圖形處理器(GPU)作為通用計(jì)算加速器的蓬勃發(fā)展,DPD已經(jīng)實(shí)現(xiàn)了在GPU上的加速運(yùn)行,達(dá)到了比CPU快幾十倍的運(yùn)算速度[9],更快的計(jì)算時(shí)間意味著可以計(jì)算更大的模擬體系,這使得該方法的應(yīng)用邊界進(jìn)一步擴(kuò)大,成為愈發(fā)有力的模擬研究工具.

在近年來(lái)的DPD研究中,邊界條件的設(shè)置一直是研究者們關(guān)心的熱點(diǎn),各種相應(yīng)的處理方法被提出以應(yīng)對(duì)特定的問(wèn)題.發(fā)展至今,DPD模擬中對(duì)邊界的處理大致分為兩類(lèi).第一類(lèi)是采用周期性邊界條件,即不存在傳統(tǒng)意義的固體邊界,一般應(yīng)用于不受限流體.在此基礎(chǔ)上進(jìn)行改進(jìn),還進(jìn)化出模擬雙Poiseuille流用以測(cè)量流體黏性的邊界方法[10]和模擬庫(kù)特流的Lees-Edwards邊界方法[11]等.第二類(lèi)方法是構(gòu)造固體的邊界,也即將問(wèn)題變成受限流體問(wèn)題.以上第一類(lèi)周期性邊界條件雖然能較好地處理某些特定的問(wèn)題,但應(yīng)用范圍很窄,大多數(shù)流動(dòng)問(wèn)題是在管流內(nèi)或者受到某些限制的,例如在微通道內(nèi)的大分子流動(dòng)[12]、多孔介質(zhì)中的多相流[13]、結(jié)構(gòu)化表面上的液滴行為[14]等,都需要用到有形的固體邊界.固體邊界條件還可細(xì)分為彈性壁面和非彈性壁面,彈性壁面指可以發(fā)生彈性變形的壁面,例如血管壁等,因?yàn)榧夹g(shù)原因,DPD目前還不能很好地處理這種移動(dòng)的彈性的固-液界面,也鮮見(jiàn)相關(guān)彈性壁面的文獻(xiàn).相比之下非彈性壁面更為簡(jiǎn)單,在模型中較易實(shí)現(xiàn),是采用最多的邊界處理方法.在以往的DPD研究中,非彈性壁面往往還被進(jìn)一步簡(jiǎn)化成沒(méi)有結(jié)構(gòu)的平板或簡(jiǎn)單的幾何形狀,究其原因主要是缺乏對(duì)復(fù)雜幾何邊界的處理方法.這使得DPD方法的應(yīng)用受限,在模擬進(jìn)一步的復(fù)雜問(wèn)題,例如紅細(xì)胞在不規(guī)則血管網(wǎng)絡(luò)中的流動(dòng)現(xiàn)象、液滴在微結(jié)構(gòu)表面的運(yùn)動(dòng)時(shí)無(wú)法準(zhǔn)確地處理這樣的邊界.一種對(duì)復(fù)雜幾何邊界自適應(yīng)的邊界處理方法將提升DPD對(duì)微流體問(wèn)題的研究能力.

2 耗散粒子動(dòng)力學(xué)的邊界條件

2.1 耗散粒子動(dòng)力學(xué)模型

DPD模型中包含一系列的質(zhì)點(diǎn)粒子,這些粒子在無(wú)網(wǎng)格的空間上運(yùn)動(dòng),彼此之間根據(jù)三種力互相碰撞: 根據(jù)勢(shì)能推導(dǎo)出來(lái)的保守力降低粒子之間切向速度的耗散力粒子連線方向上的隨機(jī)力后兩種力可以看作是配對(duì)的布朗阻尼器,和郎之萬(wàn)動(dòng)力學(xué)或者布朗動(dòng)力學(xué)不一樣,DPD的這對(duì)力是動(dòng)量守恒的.對(duì)于粗粒化的DPD系統(tǒng),布朗阻尼器是在粒子之間體現(xiàn)出黏性力的同時(shí)還要體現(xiàn)熱噪音的最簡(jiǎn)單模型.因?yàn)閯?dòng)量是守恒的,所以當(dāng)尺度達(dá)到一定量級(jí)后模型中流體的行為就完全符合宏觀的流體力學(xué)了.

在DPD中,粒子的運(yùn)動(dòng)符合牛頓第二定律,保守力、耗散力和隨機(jī)力的計(jì)算公式如下:

其中 rij=|ri-rj| 是粒子i和粒子j之間的距離;eij是連接粒子i和粒子j的單位向量; vij是粒子i和粒子j之間的相對(duì)速度; θij是高斯白噪音,同時(shí)具有對(duì)稱性,即θij=θji,以此保證整體的動(dòng)量守恒; γ和σ分別是耗散力參數(shù)和隨機(jī)力參數(shù),它們之間滿足耗散漲落定理:

其中kB是玻爾茲曼常數(shù),T是平衡溫度.簡(jiǎn)單地根據(jù)距離的衰變函數(shù)來(lái)定義保守力的權(quán)重函數(shù),

耗散力和隨機(jī)力的權(quán)重函數(shù)采用以下形式:

其中在經(jīng)典DPD中s=2,s也可以取其他值用來(lái)調(diào)節(jié)流體的黏度.

2.2 邊界條件的要求

根據(jù)長(zhǎng)期的實(shí)踐研究,人們總結(jié)出優(yōu)秀的固體邊界條件至少要滿足:

1)流體粒子不能穿透壁面;

2)流體在固體壁面處無(wú)滑移;

3)固體壁面不能影響周?chē)黧w的物理性質(zhì),如固體壁面附近的溫度、密度要和流體內(nèi)部一致,不能波動(dòng)太大.

一種簡(jiǎn)單并且被廣泛應(yīng)用的方法是在壁面的位置排布凍結(jié)的粒子,通過(guò)凍結(jié)粒子和流體粒子的作用來(lái)反映固體對(duì)流體的力.這種方法和MD中處理固體壁面的方法很像,但在DPD中會(huì)面臨新的問(wèn)題: 因?yàn)榱Ｗ又g的作用勢(shì)相對(duì)較軟流體粒子很容易穿透固體壁面粒子從而進(jìn)入固體內(nèi)部.為了解決這一明顯的錯(cuò)誤,研究者們提出了兩種方法:1)提高固-液排斥作用力; 2)設(shè)置虛擬反彈邊界.提高固-液排斥力一般通過(guò)提高固體粒子密度或者增大壁面粒子與流體粒子間的排斥力系數(shù)aij來(lái)實(shí)現(xiàn),操作起來(lái)非常簡(jiǎn)單,對(duì)于防止穿透的效果也非常好,但人為設(shè)置的強(qiáng)烈壁面排斥力雖然防止了流體的穿透,也嚴(yán)重影響了固體壁面附近流體的物理性質(zhì),例如溫度和密度,進(jìn)一步影響邊界附近的流動(dòng)行為使模擬結(jié)果偏離正確結(jié)果.設(shè)置虛擬反彈條件則不改變固-液間的作用強(qiáng)度,只對(duì)模擬過(guò)程中穿透進(jìn)入固體內(nèi)部的流體粒子進(jìn)行位置和速度的修正,主流的處理方法有以下三種: 鏡像反射、原路折回反射和麥克斯韋反射.通過(guò)反復(fù)比較和實(shí)踐,這三種方法均存在明顯的缺點(diǎn)[15],且當(dāng)壁面的形狀較復(fù)雜時(shí),難以界定粒子是否穿透進(jìn)入固體壁面內(nèi)部.

為了解決壁面附近液滴粒子的均勻性問(wèn)題,Xu和Wang[16]賦予壁面粒子相對(duì)流體粒子的虛擬速度,用于計(jì)算對(duì)流體粒子的耗散力,從而增加了流體粒子的耗散阻力實(shí)現(xiàn)了壁面的無(wú)滑移邊界條件,并得到平滑的壁面附近溫度密度分布,但仍需要依靠人為指定的虛擬邊界來(lái)界定哪些流體粒子已穿透需反彈,并不能自動(dòng)適應(yīng)復(fù)雜壁面結(jié)構(gòu).為了解決DPD對(duì)復(fù)雜幾何壁面的適應(yīng)性問(wèn)題,Mehboudi等[17]用三角形單元描述邊界形狀,這樣可以容易地從各種微機(jī)電系統(tǒng)設(shè)計(jì)圖中獲得邊界形狀,并經(jīng)測(cè)試獲得了很好的邊界效果,但因?yàn)樾枰诹Ｗ酉到y(tǒng)中引入網(wǎng)格處理的邊界方法,所以在程序簡(jiǎn)易型和計(jì)算量上仍存在推廣的壁壘.劉謀斌和常建忠[18]將整體計(jì)算域用規(guī)則背景網(wǎng)格覆蓋,分為流體區(qū)域網(wǎng)格和固體障礙區(qū)域網(wǎng)格,固體障礙區(qū)域網(wǎng)格根據(jù)周?chē)W(wǎng)格的種類(lèi)來(lái)計(jì)算法向量,實(shí)現(xiàn)了對(duì)諸如多孔介質(zhì)等問(wèn)題的邊界處理,但仍存在密度波動(dòng),并且用規(guī)則網(wǎng)格背景對(duì)復(fù)雜幾何的描述仍需要近似.Li等[19]提出了一種模擬運(yùn)行時(shí)動(dòng)態(tài)檢測(cè)粒子是否穿透壁面并動(dòng)態(tài)計(jì)算周?chē)诿娣ㄏ蛄康姆椒?這種方法在復(fù)雜微流控芯片內(nèi)通道內(nèi)展現(xiàn)出很好的適應(yīng)性,并對(duì)旋轉(zhuǎn)雙筒等動(dòng)態(tài)壁面的情況也能適應(yīng).對(duì)于固定和移動(dòng)壁面,這種方法不加區(qū)分,運(yùn)行時(shí)均需每一步都重新計(jì)算壁面法向量.

3 新型的適用復(fù)雜形狀的邊界條件

我們提出一種新型的可以應(yīng)用在任何復(fù)雜幾何形狀上的邊界條件LWNM(local wall normal method).具體地,新定義一個(gè)粒子的矢量屬性——固體壁面局部法向量(LWN),并將其記錄在每個(gè)固體粒子上,例如壁面粒子i的局部壁面法向量記為lwni.如果壁面是可移動(dòng)的,在每個(gè)時(shí)間步計(jì)算固液之間的作用力前,需要先更新每個(gè)固體粒子的法向量; 如果固體壁面是靜止的,則無(wú)需更新,只需在模擬初始化時(shí)計(jì)算一次即可,額外的計(jì)算開(kāi)銷(xiāo)接近于0.

圖1是計(jì)算壁面粒子i的法向量lwni的示意圖,其中S代表固體壁面區(qū)域,rcw是計(jì)算壁面局部法向量時(shí)的截?cái)喟霃?可以不同于粒子之間作用力的截?cái)喟霃絩c; 固體粒子j是固體粒子i的rcw范圍內(nèi)的其他固體粒子,rij是粒子j到粒子i的距離,rij是粒子j到粒子i的向量.

圖1 計(jì)算計(jì)算壁面粒子i的法向量lwniFig.1.Computing the local wall normal vector lwni of wall particle i.

計(jì)算固體粒子i的局部壁面法向量的方向向量lwnti:

將lwnti的長(zhǎng)度進(jìn)行標(biāo)準(zhǔn)化后得到壁面粒子i的壁面局部法向量lwni.:

為了判斷流體粒子是否已經(jīng)穿透壁面進(jìn)入固體區(qū)域,為每個(gè)流體粒子新增一個(gè)額外屬性φ,定義φ為流體粒子的固體體積分?jǐn)?shù),通過(guò)其周?chē)墓腆w壁面粒子的位置計(jì)算:

其中ρw是固體壁面的平均粒子密度,是一個(gè)三維的權(quán)重函數(shù),對(duì)rc范圍內(nèi)的空間進(jìn)行權(quán)重的處理,這個(gè)權(quán)重函數(shù)也被應(yīng)用在多體耗散粒子動(dòng)力學(xué)(MDPD)的保守力計(jì)算中.如圖2,當(dāng)粒子完全浸入到固體壁面中時(shí),粒子截?cái)喟霃絻?nèi)被固體粒子充滿,根據(jù)(10)式計(jì)算得φ應(yīng)為1,實(shí)際中因?yàn)楸诿媪Ｗ拥拿芏炔▌?dòng),所得φ在1左右波動(dòng)(壁面粒子密度越高,波動(dòng)范圍越小,密度越小,波動(dòng)范圍越大); 當(dāng)粒子遠(yuǎn)離固體壁面時(shí),粒子截?cái)喟霃絻?nèi)沒(méi)有固體粒子,φ為0; 當(dāng)流體粒子位于理想固-液的分界面上時(shí),半徑rc球內(nèi)的一半被固體占據(jù),另一半被流體占據(jù),根據(jù)(10)式得φ應(yīng)為0.5,如果φ ＞ 0.5即對(duì)應(yīng)流體粒子更多地進(jìn)入固體,有更多的固體粒子在截?cái)喟霃絻?nèi),如果φ ＜ 0.5即對(duì)應(yīng)流體粒子更遠(yuǎn)離固體,有更少的固體粒子在截?cái)喟霃絻?nèi).因此以φ=0.5為分界線,判斷流體粒子是否穿透固體壁面.

圖2 流體粒子固體體積分?jǐn)?shù)(φ)的四種情況(φ=0遠(yuǎn)離壁面,0 ＜ φ ＜ 0.5 未穿透壁面,0.5 ＜ φ ＜ 1.0 已穿透壁面,φ=1.0完全浸入壁面)Fig.2.Four scenarios for solidfraction (φ): φ=0 particle away from wall; 0 ＜ φ ＜ 0.5 particle near the wall; 0.5 ＜ φ＜ 1.0 particle penetrating the wall; φ=1.0 particle submerged in wall.

采取預(yù)估-修正的策略來(lái)防止流體粒子穿透壁面.在每一個(gè)時(shí)間步,預(yù)估粒子下一時(shí)刻的位置并計(jì)算φ,當(dāng)φ ＞ 0.5,即表明粒子下一個(gè)時(shí)間步會(huì)穿透壁面,需要對(duì)流體粒子速度進(jìn)行修正,修正后的新速度為

其中U和a是壁面局部的速度和加速度.如果模擬問(wèn)題中壁面靜止,(11)式可以簡(jiǎn)化成

其中en是流體粒子對(duì)應(yīng)的固體壁面局部法向量.如圖3所示,在流體粒子的截?cái)喟霃絻?nèi),有若干帶局部法向量的壁面粒子,en取其加權(quán)平均后歸一化的向量,權(quán)重函數(shù)用的是MDPD保守力中采用的三維權(quán)重函數(shù):

圖3 選用一定范圍內(nèi)固體粒子的lwni計(jì)算修正流體粒子時(shí)的法向量enFig.3.Computing en from nearby lwni for correcting penetrating fluid particle.

在對(duì)即將穿透固體壁面的流體粒子進(jìn)行速度修正時(shí),(11)式的本質(zhì)是原路返回反射方法,原路返回反射方法對(duì)壁面附近的密度和溫度波動(dòng)影響較小,但是并不能很好地滿足無(wú)滑移邊界條件.為了改進(jìn)這種方法,研究者們對(duì)此提出了很多修改方法[19—23].本文采用文獻(xiàn)[19]提出的辦法對(duì)固體和流體粒子之間的耗散力系數(shù)進(jìn)行修改.假設(shè)流體粒子距離固體壁面的距離為h,則修改后的耗散力系數(shù)為

圖4是一個(gè)不規(guī)則流體通道LWNM邊界條件實(shí)施的流程圖.操作步驟如下: 1)提取固體區(qū)域的邊界層粒子,如圖4(b),為了簡(jiǎn)化計(jì)算,只有固體邊界層粒子會(huì)被賦予局部壁面法向量,固體內(nèi)部區(qū)域的粒子在LWN創(chuàng)建完成后可以刪除,只留下殼層固體粒子,這樣可以進(jìn)一步減少系統(tǒng)總粒子數(shù),加快計(jì)算速度; 2)通過(guò)(8)和(9)式計(jì)算殼層固體粒子的LWN,每個(gè)固體粒子都有獨(dú)立的局部壁面法向量,如圖4(c),紅色箭頭即為局部壁面法向量; 3)在模擬過(guò)程中,通過(guò)(10)式預(yù)判流體粒子下一步是否會(huì)進(jìn)入固體壁面內(nèi),如果是則通過(guò)(11)式和(15)式修正粒子速度,模擬時(shí)的快照如圖4(d).

圖4 LWNM邊界條件的實(shí)施過(guò)程圖 (a)固體粒子構(gòu)造壁面; (b)提取表面層固體粒子; (c)計(jì)算LWN; (d)模擬時(shí)效果Fig.4.Workflow of LWNM: (a)Constructing wall with frozen particles; (b)identifying surface wall particles; (c)computing LWN; (d)snapshot during simulation.

LWNM不僅適合CPU計(jì)算,也同樣適合GPU計(jì)算,無(wú)論是LWN的計(jì)算還是之后φ的計(jì)算,都只需要截?cái)喟霃絻?nèi)鄰居粒子的信息,而這些信息是標(biāo)準(zhǔn)DPD模型中計(jì)算力都已經(jīng)計(jì)算過(guò)的,無(wú)論CPU計(jì)算還是GPU計(jì)算都可以在現(xiàn)有的模型上經(jīng)過(guò)少量的改動(dòng)實(shí)現(xiàn)LWNM方法.LWNM的局限性在于如果壁面是運(yùn)動(dòng)的,LWN需要不斷重新計(jì)算,而且表面固體層LWN的計(jì)算需要更多固體內(nèi)部的粒子的參與,計(jì)算效率將有所降低.

4 算例驗(yàn)證

我們?cè)赑oiseuille流動(dòng)中驗(yàn)證此邊界條件的效果.在一個(gè)三維模擬腔中,x和y方向是周期性邊界條件,在z方向布置上下兩塊平板,平板和流體間為無(wú)滑移邊界條件,對(duì)流體施加沿著x方向的體積力.根據(jù)納維-斯托克斯方程可以得出這個(gè)情況的精確解[24]:

其中d是兩板之間距離的一半,ν是運(yùn)動(dòng)黏度,F是單位質(zhì)量上的體積力.模擬中采用的參數(shù)設(shè)置如下: ρ=4,kBT=1.0,a=18.75,γ=4.5,σ=3.0,rc=rcw=1.0,F=0.02.模擬盒子的大小是30.0 × 30.0 × 34.0,其中上下壁面壁厚為 2.0,總共包含121500個(gè)流體粒子和16200個(gè)固體粒子.當(dāng)流體充分發(fā)展后,在z方向上按0.2的厚度設(shè)置統(tǒng)計(jì)層,對(duì)流體粒子的x方向上的分速度進(jìn)行統(tǒng)計(jì),統(tǒng)計(jì)結(jié)果繪制于圖5.

圖5 LWNM邊界方法統(tǒng)計(jì)結(jié)果和理論解的對(duì)比(Vx)Fig.5.Comparison of LWNM and theoretical results (Vx).

從圖5中可見(jiàn)模擬結(jié)果和(16)式給出的理論解吻合得相當(dāng)好,證明了LWNM能提供真正的無(wú)滑移邊界條件.圖6顯示了在流體區(qū)域的溫度和密度分布,可以看到壁面附近的區(qū)域無(wú)論溫度還是密度的波動(dòng)都非常小.

驗(yàn)證了LWNM邊界方法優(yōu)秀的速度無(wú)滑移控制和溫度密度控制后,需要具體展示LWNM最大的優(yōu)勢(shì),即對(duì)處理復(fù)雜幾何邊界條件時(shí)的適應(yīng)能力.以下是LWNM在兩個(gè)DPD應(yīng)用上的使用情況.

1)具有復(fù)雜結(jié)構(gòu)的超疏水表面.自然界中的很多材料常具有規(guī)律性的微結(jié)構(gòu),比如人們熟知的荷葉,表面有很多微小的突起,這些微小的突起改變了表面的親疏水性質(zhì),使材料可以達(dá)到超疏水或者超親水的浸潤(rùn)屬性.這些生物材料的表面特性給了人們啟發(fā),使得在工程應(yīng)用中也越來(lái)越多地使用帶微結(jié)構(gòu)的表面來(lái)達(dá)到特定的目的,例如除冰.當(dāng)冰在固體表面凝結(jié)冰層逐漸增厚,會(huì)導(dǎo)致很多災(zāi)難,比如高壓輸電線的變重、管道的爆裂、路面的變滑,最致命的還是飛機(jī)機(jī)翼的結(jié)冰,會(huì)直接使飛機(jī)失去升力,還曾經(jīng)導(dǎo)致了美國(guó)3407航班的空難.當(dāng)冰層已經(jīng)累積起來(lái)后除冰會(huì)非常困難,所以一般在冰層開(kāi)始凝結(jié)時(shí)就進(jìn)行干預(yù)防止凝結(jié),常用方法有加熱法和化學(xué)法,但是加熱法太浪費(fèi)能源,化學(xué)法又可能腐蝕表面.近年來(lái),具有微結(jié)構(gòu)的表面給除冰帶來(lái)了新的研究思路,通過(guò)加工表面的微結(jié)構(gòu)使表面的疏水性增強(qiáng),當(dāng)液態(tài)水沾到表面時(shí)不容易附著,在微小外力作用下很容易發(fā)生滾動(dòng)并脫落,由此防止了在固體表面上結(jié)冰.如何通過(guò)調(diào)整材料表面的微結(jié)構(gòu)而不是改變表面的化學(xué)性質(zhì)來(lái)調(diào)節(jié)親疏水性是近幾年來(lái)研究的熱點(diǎn),研究者們?cè)O(shè)計(jì)出了各種各樣的微結(jié)構(gòu),并測(cè)試它們的效果.Liu和Kim[25]還設(shè)計(jì)出了一種帶二級(jí)結(jié)構(gòu)的微結(jié)構(gòu)表面,可以使完全浸潤(rùn)材料如硅,在經(jīng)過(guò)表面微結(jié)構(gòu)的設(shè)計(jì)后達(dá)到超疏水的表面性能,實(shí)驗(yàn)證明可以使任何流體在滴落后彈開(kāi)而不附著.更進(jìn)一步還有學(xué)者通過(guò)表面微結(jié)構(gòu)的密度梯度,來(lái)控制液滴反彈的高度和方向.僅改變表面的微結(jié)構(gòu),而不用改變材料的化學(xué)性質(zhì)就能帶來(lái)如此巨大的表面性能的改變,這給材料工程帶來(lái)了一個(gè)新的思路和研究熱點(diǎn).DPD等粒子方法也常被用來(lái)模擬液滴和表面碰撞的過(guò)程,由于微結(jié)構(gòu)表面已經(jīng)不屬于簡(jiǎn)單形狀的表面,傳統(tǒng)邊界條件已經(jīng)不再適用,所以在以往研究中往往通過(guò)構(gòu)造高密度的壁面防止流體粒子的穿透,這種方法的弊端在前文中已經(jīng)闡述.采用LWNM邊界方法,可以在滿足各項(xiàng)邊界要求的條件下進(jìn)行模擬.圖7(a)展示了文獻(xiàn)中微結(jié)構(gòu)表面的實(shí)驗(yàn)照片,圖7(b)展示了應(yīng)用本邊界條件,對(duì)固體壁面粒子賦予局部法向量后的可視化結(jié)果,圖中的箭頭代表每個(gè)壁面粒子的法向量.

圖6 LWNM邊界方法統(tǒng)計(jì)結(jié)果和理論解的對(duì)比(溫度和密度)Fig.6.Comparison of LWNM and theoretical results (temperature and density).

圖7 (a)微結(jié)構(gòu)表面對(duì)液滴親疏水性的影響[26]; (b)粒子模擬中微結(jié)構(gòu)表面的LWNFig.7.(a)Microstructures on surface affect the hydrophilicity; (b)LWNs (grey arrows)of surface with microstructures.

2)復(fù)雜通道.很多應(yīng)用中需要用到幾何形狀復(fù)雜的管道,例如分叉和匯集血管中的血液流動(dòng)、微流控芯片中的復(fù)雜管路等.模擬介觀尺度的血液和其中的紅細(xì)胞是DPD的重要應(yīng)用之一,紅細(xì)胞模擬的長(zhǎng)度尺寸通常在幾到幾百個(gè)微米之間,非常適合介觀方法DPD,此外DPD是一種粒子方法,在模擬復(fù)雜結(jié)構(gòu)的流體時(shí)非常靈活,又滿足守恒定律,可以從系統(tǒng)的狀態(tài)追溯到復(fù)雜流體的相關(guān)本構(gòu)方程.這些都使DPD非常適合用來(lái)進(jìn)行紅細(xì)胞相關(guān)的模擬研究.紅細(xì)胞模擬的研究按照時(shí)間發(fā)展主要分為三個(gè)階段: 第一階段模擬單個(gè)紅細(xì)胞,主要解決紅細(xì)胞的建模問(wèn)題,包括紅細(xì)胞的變形和舒展,并解決單個(gè)紅細(xì)胞在簡(jiǎn)單流動(dòng)中的動(dòng)力學(xué)行為[27]; 第二階段模擬兩個(gè)紅細(xì)胞,主要關(guān)注兩個(gè)紅細(xì)胞之間的相互作用,包括聚集行為和解聚集行為[28]; 第三階段也正是當(dāng)前最被關(guān)注的階段,模擬大量的紅細(xì)胞在真實(shí)血管中的流動(dòng)行為,比如管流或剪切流中的運(yùn)動(dòng)行為等[29].在前兩個(gè)階段,紅細(xì)胞一般都是處于無(wú)限大不受限的流體中,或者在簡(jiǎn)單的圓管中,但是在第三個(gè)階段就必須考慮復(fù)雜的血管形狀的影響,比如血管的分叉、狹窄、以及堵塞等.對(duì)復(fù)雜的真實(shí)的血管進(jìn)行建模,并施加合適的邊界條件,是模擬最終能夠幫助醫(yī)學(xué)研究并有更廣闊應(yīng)用的前提之一.

圖8展示了應(yīng)用LWNM邊界方法進(jìn)行復(fù)雜血管生成的過(guò)程.首先要進(jìn)行血管建模,這個(gè)模型可以是由CAD軟件繪制,或是由臨床掃描圖像轉(zhuǎn)化而來(lái),如圖8(a)和圖8(b); 然后對(duì)該模型進(jìn)行網(wǎng)格布點(diǎn),所有的網(wǎng)格節(jié)點(diǎn)都將對(duì)應(yīng)粒子模型中的壁面粒子,對(duì)模型提取表面層的粒子,此步可以通過(guò)對(duì)固體粒子計(jì)算固體體積分?jǐn)?shù)來(lái)得到,如圖8(c);最后對(duì)固體粒子計(jì)算局部法向量,形成可供DPD模擬的固體壁面模型,如圖8(d).在此邊界條件的支持下,血液細(xì)胞的模擬將有很多新的問(wèn)題可以模擬研究.

圖8 復(fù)雜血管的局部法向量的生成Fig.8.LWN generation in complex blood vessel.

微流控芯片中的流體管道通常具有很多轉(zhuǎn)向和分叉匯集,我們采用LWNM邊界處理方法做了一個(gè)復(fù)雜形狀的管道.先通過(guò)貝塞爾曲線繪制TJU形狀的路徑,然后由圓形截面沿該路徑形成復(fù)雜管道,接著用固體壁面粒子填充管道外的空間,再采用LWNM邊界條件僅保留邊界層固體粒子并計(jì)算局部壁面法向量,液滴采用MDPD模型,模型參數(shù)為 All=—40,Bll=25,rd=0.75,固-液之間的參數(shù)Asl=—12,壁面對(duì)液滴的接觸角約為130°.每個(gè)液滴由1607個(gè)DPD粒子組成.模擬完成后對(duì)液滴的材質(zhì)設(shè)置成水,管道壁面的材質(zhì)設(shè)置成毛玻璃,再在頂部添加方形面光源,在管道下方添加背景底板,最后渲染得到圖9.

圖9 應(yīng)用LWNM實(shí)現(xiàn)液滴在TJU(同濟(jì)大學(xué))形狀的復(fù)雜管道中運(yùn)動(dòng)Fig.9.Droplets in TJU (Tongji University)shape tube with LWNM.

模擬結(jié)果表明,對(duì)這種非常不規(guī)則的管道形狀,LWNM邊界處理方法仍表現(xiàn)出了很強(qiáng)的適應(yīng)能力,對(duì)微流控芯片中復(fù)雜管道用粒子方法進(jìn)行模擬研究提供了有力的輔助.

5 結(jié) 論

本文提出一種適用于復(fù)雜幾何壁面的邊界條件LWNM.該方法通過(guò)對(duì)壁面粒子的周?chē)Ｗ拥奈恢眠M(jìn)行加權(quán)計(jì)算,得到壁面局部的法向量,從而在流體粒子穿透壁面時(shí)提供可靠的位置速度修正方向; 通過(guò)定義流體粒子的周?chē)腆w體積分?jǐn)?shù)來(lái)判斷流體粒子是否已經(jīng)穿透壁面,如果流體粒子的φ大于閾值則認(rèn)為流體已經(jīng)穿透固體,并根據(jù)壁面的局部法向量施加垂直壁面向外的作用力; 同時(shí)通過(guò)對(duì)固液粒子之間黏滯力的修改實(shí)現(xiàn)無(wú)滑移邊界條件.通過(guò)Poiseuille流算例證明了LWNM邊界條件可以達(dá)到速度無(wú)滑移條件,并能出色地控制壁面附近流體的密度和溫度波動(dòng); 通過(guò)帶微結(jié)構(gòu)的表面和復(fù)雜管道的例子,展示了LWNM對(duì)多應(yīng)用的廣闊前景.