適用復(fù)雜幾何壁面的耗散粒子動力學(xué)邊界條件*

林晨森 陳碩? 肖蘭蘭

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

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

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

1 引 言

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

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

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

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

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

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

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

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

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

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

2.2 邊界條件的要求

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

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

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

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

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

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

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

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

圖1是計(jì)算壁面粒子i的法向量lwni的示意圖,其中S代表固體壁面區(qū)域,rcw是計(jì)算壁面局部法向量時的截?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的長度進(jìn)行標(biāo)準(zhǔn)化后得到壁面粒子i的壁面局部法向量lwni.:

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

其中ρw是固體壁面的平均粒子密度,是一個三維的權(quán)重函數(shù),對rc范圍內(nèi)的空間進(jìn)行權(quán)重的處理,這個權(quán)重函數(shù)也被應(yīng)用在多體耗散粒子動力學(xué)(MDPD)的保守力計(jì)算中.如圖2,當(dāng)粒子完全浸入到固體壁面中時,粒子截?cái)喟霃絻?nèi)被固體粒子充滿,根據(jù)(10)式計(jì)算得φ應(yīng)為1,實(shí)際中因?yàn)楸诿媪Ｗ拥拿芏炔▌?所得φ在1左右波動(壁面粒子密度越高,波動范圍越小,密度越小,波動范圍越大); 當(dāng)粒子遠(yuǎn)離固體壁面時,粒子截?cái)喟霃絻?nèi)沒有固體粒子,φ為0; 當(dāng)流體粒子位于理想固-液的分界面上時,半徑rc球內(nèi)的一半被固體占據(jù),另一半被流體占據(jù),根據(jù)(10)式得φ應(yīng)為0.5,如果φ ＞ 0.5即對應(yīng)流體粒子更多地進(jìn)入固體,有更多的固體粒子在截?cái)喟霃絻?nèi),如果φ ＜ 0.5即對應(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ù)估-修正的策略來防止流體粒子穿透壁面.在每一個時間步,預(yù)估粒子下一時刻的位置并計(jì)算φ,當(dāng)φ ＞ 0.5,即表明粒子下一個時間步會穿透壁面,需要對流體粒子速度進(jìn)行修正,修正后的新速度為

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

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

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

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

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

圖4 LWNM邊界條件的實(shí)施過程圖 (a)固體粒子構(gòu)造壁面; (b)提取表面層固體粒子; (c)計(jì)算LWN; (d)模擬時效果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ì)算,無論是LWN的計(jì)算還是之后φ的計(jì)算,都只需要截?cái)喟霃絻?nèi)鄰居粒子的信息,而這些信息是標(biāo)準(zhǔn)DPD模型中計(jì)算力都已經(jīng)計(jì)算過的,無論CPU計(jì)算還是GPU計(jì)算都可以在現(xiàn)有的模型上經(jīng)過少量的改動實(shí)現(xiàn)LWNM方法.LWNM的局限性在于如果壁面是運(yùn)動的,LWN需要不斷重新計(jì)算,而且表面固體層LWN的計(jì)算需要更多固體內(nèi)部的粒子的參與,計(jì)算效率將有所降低.

4 算例驗(yàn)證

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

其中d是兩板之間距離的一半,ν是運(yùn)動黏度,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個流體粒子和16200個固體粒子.當(dāng)流體充分發(fā)展后,在z方向上按0.2的厚度設(shè)置統(tǒng)計(jì)層,對流體粒子的x方向上的分速度進(jìn)行統(tǒng)計(jì),統(tǒng)計(jì)結(jié)果繪制于圖5.

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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