適用復雜幾何壁面的耗散粒子動力學邊界條件*

林晨森 陳碩? 肖蘭蘭

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

2)(上海工程技術大學機械與汽車工程學院,上海 201620)

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

1 引 言

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

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

2 耗散粒子動力學的邊界條件

2.1 耗散粒子動力學模型

DPD模型中包含一系列的質點粒子,這些粒子在無網格的空間上運動,彼此之間根據三種力互相碰撞: 根據勢能推導出來的保守力降低粒子之間切向速度的耗散力粒子連線方向上的隨機力后兩種力可以看作是配對的布朗阻尼器,和郎之萬動力學或者布朗動力學不一樣,DPD的這對力是動量守恒的.對于粗粒化的DPD系統,布朗阻尼器是在粒子之間體現出黏性力的同時還要體現熱噪音的最簡單模型.因為動量是守恒的,所以當尺度達到一定量級后模型中流體的行為就完全符合宏觀的流體力學了.

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

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

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

耗散力和隨機力的權重函數采用以下形式:

其中在經典DPD中s=2,s也可以取其他值用來調節流體的黏度.

2.2 邊界條件的要求

根據長期的實踐研究,人們總結出優秀的固體邊界條件至少要滿足:

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

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

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

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

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

3 新型的適用復雜形狀的邊界條件

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

圖1是計算壁面粒子i的法向量lwni的示意圖,其中S代表固體壁面區域,rcw是計算壁面局部法向量時的截斷半徑,可以不同于粒子之間作用力的截斷半徑rc; 固體粒子j是固體粒子i的rcw范圍內的其他固體粒子,rij是粒子j到粒子i的距離,rij是粒子j到粒子i的向量.

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

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

將lwnti的長度進行標準化后得到壁面粒子i的壁面局部法向量lwni.:

為了判斷流體粒子是否已經穿透壁面進入固體區域,為每個流體粒子新增一個額外屬性φ,定義φ為流體粒子的固體體積分數,通過其周圍的固體壁面粒子的位置計算:

其中ρw是固體壁面的平均粒子密度,是一個三維的權重函數,對rc范圍內的空間進行權重的處理,這個權重函數也被應用在多體耗散粒子動力學(MDPD)的保守力計算中.如圖2,當粒子完全浸入到固體壁面中時,粒子截斷半徑內被固體粒子充滿,根據(10)式計算得φ應為1,實際中因為壁面粒子的密度波動,所得φ在1左右波動(壁面粒子密度越高,波動范圍越小,密度越小,波動范圍越大); 當粒子遠離固體壁面時,粒子截斷半徑內沒有固體粒子,φ為0; 當流體粒子位于理想固-液的分界面上時,半徑rc球內的一半被固體占據,另一半被流體占據,根據(10)式得φ應為0.5,如果φ ＞ 0.5即對應流體粒子更多地進入固體,有更多的固體粒子在截斷半徑內,如果φ ＜ 0.5即對應流體粒子更遠離固體,有更少的固體粒子在截斷半徑內.因此以φ=0.5為分界線,判斷流體粒子是否穿透固體壁面.

圖2 流體粒子固體體積分數(φ)的四種情況(φ=0遠離壁面,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.

采取預估-修正的策略來防止流體粒子穿透壁面.在每一個時間步,預估粒子下一時刻的位置并計算φ,當φ ＞ 0.5,即表明粒子下一個時間步會穿透壁面,需要對流體粒子速度進行修正,修正后的新速度為

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

其中en是流體粒子對應的固體壁面局部法向量.如圖3所示,在流體粒子的截斷半徑內,有若干帶局部法向量的壁面粒子,en取其加權平均后歸一化的向量,權重函數用的是MDPD保守力中采用的三維權重函數:

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

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

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

圖4 LWNM邊界條件的實施過程圖 (a)固體粒子構造壁面; (b)提取表面層固體粒子; (c)計算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計算,也同樣適合GPU計算,無論是LWN的計算還是之后φ的計算,都只需要截斷半徑內鄰居粒子的信息,而這些信息是標準DPD模型中計算力都已經計算過的,無論CPU計算還是GPU計算都可以在現有的模型上經過少量的改動實現LWNM方法.LWNM的局限性在于如果壁面是運動的,LWN需要不斷重新計算,而且表面固體層LWN的計算需要更多固體內部的粒子的參與,計算效率將有所降低.

4 算例驗證

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

其中d是兩板之間距離的一半,ν是運動黏度,F是單位質量上的體積力.模擬中采用的參數設置如下: ρ=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個固體粒子.當流體充分發展后,在z方向上按0.2的厚度設置統計層,對流體粒子的x方向上的分速度進行統計,統計結果繪制于圖5.

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

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

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

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

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

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

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

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

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

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

圖9 應用LWNM實現液滴在TJU(同濟大學)形狀的復雜管道中運動Fig.9.Droplets in TJU (Tongji University)shape tube with LWNM.

模擬結果表明,對這種非常不規則的管道形狀,LWNM邊界處理方法仍表現出了很強的適應能力,對微流控芯片中復雜管道用粒子方法進行模擬研究提供了有力的輔助.

5 結 論

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