一種基于物質點法的熔冰仿真方法

摘要：為了解決拉格朗日法模擬熔冰時產生的穿模現象，提出一種基于物質點法的熔冰仿真方法。首先，將歐拉空間中的每個拉格朗日粒子視為固-液雙相耦合的集合體，實現固相到液相的平穩過渡；其次，在物質點法背景網格上計算熱能，采用預處理共軛梯度法求解相變過程的溫度線性系統；最后，對潛熱現象進行處理，引入虛擬水模型，通過限制虛擬水的移動來實現對冰塊外部水流現象的仿真。實驗結果表明，該方法不僅能夠利用物質點法處理自碰撞的優勢解決穿模現象，而且能模擬出真實的潛熱過程的水流現象。

關鍵詞：熔冰現象；物質點法；虛擬水；水流處理

中圖分類號：TP391.9文獻標識碼： ADOI：10.3969/j.issn.1007-791X.2024.04.0050引言

近年來，基于物理的流體、固體以及各種相之間的交互耦合模擬已經成為計算機圖形學和虛擬現實領域的研究熱點。趙靜等［1］已經對流固交互仿真做了全面的介紹，熔冰現象作為流固交互現象的特例，表現為水從固相轉變為液相。由于相變過程中表現出復雜的力學特性（如不可壓縮性、非線性本構模型解、不同相之間的平穩轉換等問題），如何高效、真實地實現相變的物理仿真，成為計算機圖形學研究領域的難點和熱點。

2002年Carlson等［2］對MAC算法進行改進，將固體和接近固體的材料處理為非常高粘度的流體，首次完成了熔融流動的仿真；2005年Keiser等［3］將固體力學與N-S方程合并，采用統一的處理器制作了流體、固體及相變的動畫，為相變的仿真提供了新的方向；2010年Iwasaki等［4］提出了一種基于粒子的模擬方法對熔水粒子與冰粒子之間的熱傳導、交互耦合進行了仿真，通過計算水-冰、水-水粒子間的界面力，對冰與水滴之間的流動進行了處理，實現了快速熔冰模擬，但是并沒有對潛熱過程進行處理，在物理世界中，物體從固相到液相的轉移中，會吸收能量用來打開分子鍵，但是物體的溫度卻不會變化，也就是潛熱過程。如果不能對潛熱階段進行良好處理，可能會導致整個系統能量不守恒；2014年Lii等［5］在 Iwasaki的基礎上，對潛熱過程進行了處理，他們使用分區的方法，對于固體和流體分別迭代計算，并利用其中一個確定另一個的邊界條件，這種方法雖然簡單，但是需要很小的時間步長才能獲得穩定的效果，效率較低；2017年Yang等［6］將Helmholtz自由能應用在模擬相變中，使用一個額外的變量來分別表示不同階段，提供了相之間的連續界面，從能量的角度實現了固-液相變。物質點法是一種伽遼金方法，使用統一的方式處理固體和流體，使用拉格朗日粒子攜帶物質的信息，并用歐拉網格作為背景來計算粒子所受力。物質點法本身支持自動的基于網格的碰撞處理和雙向耦合，已經成功地應用于各種物質特性的模擬，包括顆粒狀物質［7］、毛發［8］，相變也不例外。2014年Stomakhin等［9］在物質點法背景網格上求解熱方程，開發了非通用的Chorin樣式投影，實現了基于物質點法的相變模擬，為采用物質點法解算器模擬熔冰現象奠定了理論基礎；2019年Ding等［10］開發了一個熱力模型，利用混合理論解決了水、氣體和面團之間的相互作用，對日常生活中的烘焙面包、煎餅等相變現象進行了模擬；2021年Su等［11］將POM-POM模型引入計算機圖形學，設計了統一的粘性本構模型，引入二階精度廣義單步單解格式，同時將粘彈性液體模擬器與非傅立葉熱擴散解算器耦合，成功模擬出相變現象（如巧克力融化），進一步證明了物質點法在模擬相變現象中的靈活性與可擴展性；2021年Chen等［12］使用能量自動捕捉表面梯度，并且采用粒子重采樣的方法對表面進行捕捉，成功模擬出蠟燭燃燒等現象，但是Chen等人的方法專注于對平穩相變的模擬，并沒有對冰塊熔化過程的表面水流進行處理。

第4期趙靜等一種基于物質點法的熔冰仿真方法

燕山大學學報2024

一系列擴展于拉格朗日粒子法的求解器，盡管已經取得了較可觀的效果，但是這些方法需要解決來自于非線性固體本構模型的復雜數值計算問題，通常使用線性近似求解，并且需要外部的迭代以收斂到非線性的精確解，這需要耗費巨大的算力，效率很低。基于上述情況，物質點法成為了一個十分具有吸引力的解決方案，與傳統的拉格朗日法相比，基于物質點法的求解器不需要額外處理粒子間的自碰撞問題，并且可以使用物質點法網格直接計算線性系統，保證計算效率。因此，選擇物質點法作為基礎求解器，同時對潛熱過程進行處理。

1熔冰過程溫度線性系統求解

1.1基于物質點法的固液模型

物質點法提供了一種統一的粒子模擬框架，可以很容易地實現流體材料與固體材料的耦合，保證固相到液相的平穩過渡，因此，采用物質點法作為基礎的處理器進行固液模型的構建。物質點法作為一種混合的歐拉-拉格朗日方法，在運動過程中歐拉網格不變，使用拉格朗日粒子表示材料并進行移動，以保留運動過程中的細節。同時，物質點法的混合特性允許使用規則的笛卡爾網格自動處理碰撞現象，而不需要額外的計算，可以有效地處理熔冰時冰塊與熔水之間的自碰撞問題。

使用X∈Ω0表示物體坐標下的材料點，x∈Ωt表示世界坐標。是一種映射關系，滿足x=（X，t）。物質運動由質量、動量守恒控制。控制方程如下：

ρDρDt=0ρDvDt=·σ+ρg，（1）

式中，ρ表示密度，DρDt表示流體密度的物質導數，v表示速度，σ表示柯西應力，g表示重力加速度，P表示皮奧拉基爾霍夫應力。

1.2熱傳導模擬

以構建的物質點法固液模型為基礎，將歐拉空間中的每個拉格朗日粒子視為固-液雙相耦合的集合體，為粒子設置溫度屬性。為了更真實模擬出冰塊受熱熔化的現象，主要考慮兩個熱能來源：1）拉格朗日粒子間的熱傳導；2）來自空氣的熱能。示意圖如圖1所示。

熱傳導的控制方程為

Tt=k2T，（2）

式中，T表示溫度，t表示時間，k表示熱擴散系數，2表示拉普拉斯算子。

在本文提出的系統中，在歐拉網格上求解熱傳導產生的線性系統，因此對于拉格朗日粒子之間的熱傳導不需要做額外的處理。對于來自空氣溫度的熱能，將周圍空氣的溫度設定為一個恒定的環境溫度，不斷地向冰體發射熱光子，單位時間內熱源產生的熱能由熱光子攜帶：

E=εδT′4HΔtN，（3）

式中，E表示熱能，ε表示發射率因子，δ表示Stefan-Boltzmann常數，T′表示熱源的溫度，H表示熱源的面積，Δt表示光子輻射的時間間隔，N表示每個時間步發射的光子數。

1.3溫度線性系統

熱傳導的過程遵循熱力學方程，在求解熱力學方程的過程中會產生一個復雜的線性系統。物質點法允許基于網格的隱式積分求解，其條件獨立于拉格朗日粒子的數量，因此使用物質點法背景網格對熱力學方程產生的線性系統進行求解。通過添加一層物質點法網格，將拉格朗日粒子的溫度插值到增設的物質點法網格上，在網格上使用一種改進矩陣條件數的預處理共軛梯度法對溫度泊松方程進行迭代求解，求解結束后將溫度信息返回到拉格朗日粒子。

將式（2）進行基于空間的離散化：

Ti，jt=kΔx2（-4Ti，j+Ti+1，j+Ti-1，j+Ti，j+1+Ti，j-1），（4）

式中，i，j表示序號為i，j的網格節點，Δx表示網格間距，k表示熱擴散率。

為了允許更大的時間步長，采用隱式時間積分器對溫度進行求解，將式（2）進行基于隱式時間積分的離散化：

Tn+1-TnΔt=kΔx2DTn+1，（5）

式中，n表示當前時間步，n+1表示下一時間步，D表示系數矩陣，T表示網格溫度矩陣。

2熔冰過程中潛熱現象繪制

2.1虛擬水

冰與水的密度不同，因此一定數量的冰粒子并不能轉換為等數量的水粒子。此外，如果冰粒子到達熔點后，忽略潛熱過程，直接將冰粒子轉換為水粒子，可能導致能量不守恒。為此，引入虛擬水覆蓋在冰的表面，介于冰粒子和水粒子之間的一種形態，作為冰粒子與水粒子交互耦合的橋梁（如圖2所示）。將冰粒子到水粒子的相變拆分為兩個部分：冰粒子到虛擬水的轉換、虛擬水到水粒子的轉換。

冰粒子持續吸收熱能，用以打開粒子中水分子鍵，設定一個冰粒子由大量的水分子組成，當冰粒子的溫度達到熔點之后繼續吸收熱能，分子鍵打開的那部分水分子以虛擬水的形態包裹在冰粒子的周圍，直至冰粒子中所有的分子鍵被打開。

一個冰粒子完全轉換為水粒子需要的潛熱能量為

Etotal=miceEw，（6）

式中，mice表示一個冰粒子的質量，Ew表示單位質量的晶體變成同溫度的液態物質所需吸收的熱能。

當冰粒子的溫度低于熔點時，吸收熱量用以升溫；溫度高于熔點之后，吸收的熱量增加了冰粒子的潛熱。因此，到達熔點之后，冰粒子轉化為虛擬水的占比為

ΔEratio=ΔEEtotal，（7）

式中，ΔEratio表示當前時間步冰粒子轉化為虛擬水的占比，ΔE表示當前時間步吸收的用于打開分子鍵的熱能。如果持續吸收的能量超過Etotal，則一個冰粒子完全轉換為水粒子，因此ΔEratio的取值為［0，1］。

虛擬水的體積代表了圍繞著冰粒子的薄薄一層水，其計算公式為

Δvv，water=γviceΔEratio，（8）

式中，Δvv，water表示虛擬水在當前時間步的體積變化量，γ表示同等質量的冰粒子與水粒子的體積比，一般取值為0.92，vice表示冰粒子體積。

當虛擬水的體積大于一個水粒子的體積時，動態地生成水粒子。

2.2表面水流

由于內部冰粒子并不直接與空氣接觸，其吸收的熱能主要來源于外部冰粒子的熱傳導，因此內部冰粒子不含任何的虛擬水，虛擬水主要聚集于吸收熱能較多的外部冰粒子。引入虛擬水模型之后，通過虛擬水的移動來實現表面水流現象。從兩個維度對虛擬水的運動進行限制：

1） 垂直方向：由于重力，虛擬水會移動到較低區域。

2） 水平方向：如果一個粒子的虛擬水比另一個粒子多，那么該粒子的虛擬水會移動到同一水平的另一個粒子，相當于對周圍鄰居的水量進行一個平均的操作。

虛擬水從一個冰粒子轉移到另一個冰粒子的過程中，體積總變化量為零。

3邊界碰撞

由于基于粒子的碰撞會使形變梯度產生不一致性，導致模擬中出現偽影現象。因此在物質點法解算器中需要將邊界條件應用在歐拉網格上。物質點法中常用的邊界處理有：粘性邊界條件、滑移邊界條件。這兩類邊界條件皆適用于本文的系統。

使用粘性邊界條件對邊界進行限制時，將邊界之外的網格速度賦值為零，其表達式為

vt+1i=BCsticky（v^t+1）=0，（9）

式中，t+1表示下一時間步，v^t+1表示碰撞尚未結束的網格速度，在更新粒子位置之前，可以多次對粒子進行碰撞處理，BCsticky為粘性邊界碰撞處理函數。

使用滑移邊界條件進行限制時，將垂直于邊界的速度賦值為零，只保留平行于邊界的速度，其表達式為

vt+1i=BCslip（v^t+1）=v^t+1-n（nTv^t+1），（10）

式中，n表示曲面法線，BCslip為滑移邊界碰撞處理函數。

4實驗結果與分析

所有實驗皆在Windows 10操作系統下基于C++、Python與Taichi編程語言進行搭建，并通過Houdini對實驗效果進行渲染，其中，C++編程語言的代碼開發環境為Visual Studio 2019，Python 3.7與Taichi 0.9.2，實驗中的流體與固體模型使用斯坦福開源的三維代碼、3Dmax與Maxon Cinema 4D進行建模。硬件環境為Intel Core i5-10500 CPU @ 3.10 GHz，主顯卡為Nvidia GeForce RTX2060，顯存容量為6 GB。實驗的具體參數見表1。表1實驗場景數據表

Tab.1Table of experimental scenario data

圖序號粒子數k求解器時間步長/ms楊氏模量/MPa速率/（s/frame）3a50MPM1.152×10-30.0283b50SPH1.15—0.0044a50MPM2.46×1020.064b50MPM2.46×1020.05560MPM0.82×10-31.026a780MPM3.375×10-34.986b780MPM3.375×10-33.68圖3是冰塊熔化的對比圖，實驗中采用“S”字母形狀的冰塊進行效果展示，紅色框圖為細節展示。圖3（b）為傳統粒子法實現的冰塊熔化效果，出現了明顯的穿模現象（如圖3（b）紅色框圖所示）并且不會形成表面水流。圖3（a）為本文方法實現的熔冰效果，由實驗結果可知，通過使用基于物質點法的仿真方法，可以解決粒子穿模現象，同時引入虛擬水可以模擬出真實的表面水流現象（如圖3（a）紅色框圖所示），證明了本文方法的可行性。

圖4分別展示了文獻［4］和本文方法對穿模現象的處理。為了更好地展示冰粒子與水粒子的位置關系，采用紫色粒子表示水粒子，右圖為左圖紅色框區的放大展示。由圖4（a）可以看出，文獻［4］采用的方法存在明顯的粒子穿模問題；而本文采用基于物質點法的求解器，可以很好地解決該問題，模擬出更真實的熔冰過程。

圖5為線性系統求解的對比圖，周圍空氣的溫度作為一個恒定的環境溫度不斷地向冰體發射熱光子，粒子初始顏色為藍色，當受熱到達熔點之后，粒子顏色變為紅色。由圖5（b）中第50幀與第55幀的對比可以看出，采用未改進矩陣條件數的方法對線性系統進行求解時，一些粒子到達熔點之后又重新回到了未達到熔點時的狀態，出現了不合理的視覺效果。而采用改進矩陣條件數的思想預先對線性系統進行處理，避免了這一問題，同時減少了線性系統的迭代次數。

圖6分別展示了文獻［9］和本文方法對表面水流現象的處理。文獻［9］主要聚焦于固體與流體之間的平穩轉變，并沒有對熔冰過程的表面水流進行處理。由圖6（a）可以看出，文獻［9］的方法并不會形成水流效果，而本文通過限制虛擬水的移動對冰粒子表面的薄薄水層進行了模擬，獲得了更合理的視覺效果。

5結論

提出一種基于物質點法的求解器對熔冰現象進行了仿真，利用物質點法處理自碰撞的優勢解決了粒子法求解器中常出現的穿模現象，使用背景網格對溫度線性系統進行了求解，減小了系統的冗余度，提高了仿真的速度。同時引入虛擬水模型對熔冰的表面水流進行了模擬，從實驗效果可以看出，添加表面水流的仿真效果更具真實性。

目前只考慮了小規模的冰塊熔化模擬，由于物質點法的局限性，使得處理大規模的熔化仿真變得困難，頻繁的粒子與網格之間的信息交換，降低了系統的計算速度。不過，隨著機器學習算法與神經網絡的流行與發展，將當前模型與神經網格、機器學習相結合可能是個未來解決方案。

參考文獻

［1］趙靜，許立群，安泳鋼，等.基于物理的流固交互仿真方法綜述［J］.燕山大學學報，2022，46（5）：385-397.

ZHAO J，XU L Q，AN Y G，et al.Review of physics based fluid solid interaction simulation methods［J］.Journal of Yanshan University，2022，46（5）：385-397.

［2］CARLSON M，MUCHA P J，VAN HORN III R B，et al.Melting and flowing［C］// Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer Animation，San Antonio，USA，2002：167-174.

［3］KEISER R，ADAMS B，GASSER D，et al.A unified lagrangian approach to solid-fluid animation［C］// Proceedings Eurographics/IEEE VGTC Symposium Point-Based Graphics，Stony Brook，USA，2005：125-148.

［4］IWASAKI K，UCHIDA H，DOBASHI Y，et al.Fast particle-based visual simulation of ice melting［J］.Computer Graphics Forum，2010，29（7）：2215-2223.

［5］LII S Y，WONG S K.Ice melting simulation with water flow handling［J］.The Visual Computer，2014，30（5）：531-538.

［6］YANG T，CHANG J，LIN M C，et al.A unified particle system framework for multi-phase，multi-material visual simulations［J］.ACM Transactions on Graphics，2017，36（6）：224.

［7］KLR G，GAST T，PRADHANA A，et al.Drucker-prager elastoplasticity for sand animation［J］.ACM Transactions on Graphics，2016，35（4）：103.

［8］HAN X，GAST T F，GUO Q，et al.A hybrid material point method for frictional contact with diverse materials［J］.Proceedings of the ACM on Computer Graphics and Interactive Techniques，2019，2（2）：17.

［9］STOMAKHIN A，SCHROEDER C，JIANG C，et al.Augmented MPM for phase-change and varied materials［J］.ACM Transactions on Graphics，2014，33（4）：138.

［10］DING M，HAN X，WANG S，et al.A thermomechanical material point method for baking and cooking［J］.ACM Transactions on Graphics，2019，38（6）：192.

［11］SU H，XUE T，HAN C，et al.A unified second-order accurate in time MPM formulation for simulating viscoelastic liquids with phase change［J］.ACM Transactions on Graphics，2021，40（4）：119.

［12］CHEN J，KALA V，MARQUEZ-RAZON A，et al.A momentum-conserving implicit material point method for surface tension with contact angles and spatial gradients［J］.ACM Transactions on Graphics，2021，40（4）：111.

A simulation method for melting ice based on material point method

ZHAO Jing1，2，SUN Mengmeng1，2，WANG Feng1，2，TANG Yong1，2

（1. School of Information Science and Engineering，Yanshan University，Qinhuangdao，Hebei 066004，China;

2. The Key Laboratory for Computer Virtual Technology and System Integration of Hebei Province，

Yanshan University，Qinhuangdao，Hebei 066004，China）

Abstract： In order to solve the problem of ice penetration caused by Lagrangian method，a method of ice penetration simulation based on material point method was proposed.Firstly，each Lagrangian particle in Euler space is regarded as an aggregate of solid-liquid two-way coupling to achieve a smooth transition from solid to liquid phase.Secondly，the thermal energy is calculated on the material point method background grid and the temperature linear system of the phase transition process is solved by pretreated conjugate gradient method.Finally，the latent heat phenomenon is treated，and the virtual water model is introduced to realize the simulation of the external flow phenomenon of ice by limiting the movement of the virtual water.The experimental result shows that this method can not only use material point method method to deal with the self collision advantage to solve the interpenetration，but it also simulates the real latent heat process of water flow phenomenon.

Keywords：" ice melting;material point method;virtual water;water flow treatment