王江坤，何坤金*，曹紅飛，王金強，張 燕

（1.河海大學物聯(lián)網(wǎng)工程學院，江蘇常州213022； 2.河海大學企業(yè)管理學院，江蘇常州213022）

（?通信作者電子郵箱kjinhe@163.com）

0 引言

隨著虛擬現(xiàn)實技術的發(fā)展，人們的消費方式產(chǎn)生了相應的變化，多數(shù)消費者不愿通過線下門店直接了解產(chǎn)品，而是希望通過虛擬場景還原產(chǎn)品的工作過程，從而到達對產(chǎn)品的多方面了解。選購太陽能熱水器時同樣如此，用戶希望通過自己設定環(huán)境溫度、加熱溫度、注水量等參數(shù)模擬熱水器的工作過程，以仿真得到的數(shù)據(jù)直觀地了解熱水器的信息。熱水器仿真的核心即為模擬內部的水流以及水流的加熱過程。通過計算機仿真模擬水流并對水流加熱過程進行分析，能夠有效把握水流加熱狀態(tài)的合理性，從而模擬水流的熱運動，真實還原熱水器的工作原理，使用戶通過交互控制少量系統(tǒng)定義的參數(shù)、修改加熱條件，實現(xiàn)可交互的水流加熱仿真。

近年來，計算機流體動力學（Computational Fluid Dynamics，CFD）領域發(fā)展迅速，Myuller等［1］將光滑粒子流體動力學（Smoothed Particle Hydrodynamics，SPH）方法［2-3］引入計算機圖形學領域模擬流體，奠定了SPH方法流體模擬的基礎，但模擬的流體有較為明顯的壓縮性；為改善流體的壓縮性，Becker等［4］提出了弱可壓縮的 SPH（Weakly Compressible SPH，WCSPH）方法，模擬出更為真實的流體效果，該方法被用于模擬海嘯［5］、潰壩［6］、液滴碰撞［7］等場景；Grenier等［8］和Tartakovsky等［9］為模擬多種性質流體共存現(xiàn)象，在粒子間添加不同形式的排斥力，分別建立了氣-液兩相流與液-液兩相流模型；Xiong等［10］利用兩相流模型處理水加熱過程水蒸氣和液體兩相粒子，模擬流體表面的沸騰現(xiàn)象；李世杰等［11］改進兩相流模型并應用在農(nóng)業(yè)生產(chǎn)水射流上，模擬水射流對土壤沖蝕，具有一定的工程指導意義；常元章等［12］通過對傳統(tǒng)納維斯托克斯方程添加彈性應力，統(tǒng)一牛頓流體與粘彈性流體，通過調節(jié)參數(shù)可以模擬不同粘彈性、不同性質的流體。

在處理邊界碰撞問題時，Schechter等［13］在邊界動態(tài)生成一組虛粒子對邊界粒子施加相互作用力，防止粒子穿透邊界，但動態(tài)生成的虛粒子增加了算法的時間復雜性，降低了系統(tǒng)實時性；為優(yōu)化邊界虛粒子的生成，劉旭等［14］使用泊松采樣算法對容器邊界進行采樣，深入分析三維自由表面模型，快速定位虛粒子生成位置；袁志勇等［15］通過引入平均曲率與構建水平集網(wǎng)格用于邊界修正，更加準確地處理邊界問題，具有較好的魯棒性；Hu等［16］使用邊界虛粒子方法操作一個容器中的流體倒入另一個容器，對化學實驗等高風險操作具有重要的實用意義。

在處理流體加熱問題時，可以采用基于神經(jīng)網(wǎng)絡［17］的溫度預測模型，對大量實驗數(shù)據(jù)分析進行流體溫度預測，雖然具有較高的精度，但加熱過程流體無流動感且真實性較差；Cleary等［18］首次將熱傳導方程應用于SPH方法，建立了一般的流體熱傳導數(shù)值模型，該方法已被用于模擬鑄造［19-21］和熔化金屬液滴固化表面［22］，具有較高的流動性；Farrokhpanah等［23］在SPH相變過程中添加潛熱方程，模擬具有相變的瞬態(tài)熱傳導過程，提高了熱傳導計算精度。目前的研究成果促進了SPH流體熱傳導方面的研究，但上述方法主要研究了金屬溶液固化/熔化過程流體溫度與粘度的關系，忽略了流體內部熱運動，僅適用于模擬物質的相變過程。對于流體加熱（非相變）過程而言，流體的運動主要由加熱過程中不均勻受熱產(chǎn)生的熱運動為主，并且伴隨著流體熱運動，流體內部存在對流給熱過程，又會反過來影響流體熱運動。

為此，本文提出一種基于SPH的流體熱運動仿真方法，著重研究SPH方法求解熱傳導問題時流體內部的熱傳遞與熱運動過程。首先，基于SPH方法將水流粒子化，建立水粒子間約束關系，以及粒子與邊界的碰撞檢測方法；其次，對粒子設置溫度屬性，在給定溫度場下加熱粒子，并模擬加熱過程中粒子的熱運動；最后，提供自定義系統(tǒng)參數(shù)接口，便于用戶調節(jié)加熱條件，實現(xiàn)可交互的控制水流的加熱過程。

1 加熱仿真系統(tǒng)框架

基于SPH方法的交互式水流加熱仿真系統(tǒng)主要包括SPH水流物理模型的構建、水粒子與容器的邊界碰撞檢測、水流的熱運動模擬以及參數(shù)化控制等，其流程如圖1所示。

根據(jù)上述設計流程，水流的加熱仿真步驟如下：

步驟1 水流物理模型的建立?；赟PH方法將連續(xù)水流離散為一系列粒子，建立粒子間的約束關系，分析粒子受力并更新每一時刻粒子的運動狀態(tài)。

步驟2 粒子與容器邊界的碰撞處理。對容器的網(wǎng)格進行采樣，獲取邊界信息，并將離散后的水粒子注入三維容器中，判斷粒子是否與邊界網(wǎng)格碰撞，若發(fā)生碰撞則進行碰撞處理。

步驟3 粒子熱運動模擬。實時計算容器內溫度場，加熱水粒子，對不同溫度的粒子添加外力，以模擬水流熱運動。

步驟4 用戶可交互加熱。設置系統(tǒng)加熱參數(shù)接口，使用戶可以通過修改少量系統(tǒng)參數(shù)改變加熱仿真條件。

其中，步驟1和步驟3是本文研究的主要內容。首先，通過以上流程實例化SPH水粒子，建立粒子間約束關系以及粒子與容器的碰撞檢測方法，使粒子運動限定在容器內；然后，實時計算容器內的溫度場加熱水粒子，根據(jù)粒子溫差施加外力模擬水流熱運動；最后，通過編輯修改加熱仿真參數(shù)，實現(xiàn)用戶對水流加熱仿真的可交互控制。

圖1 本文方法流程Fig.1 Flowchart of theproposed method

2 流體仿真模型

為了在計算機上實現(xiàn)水流在容器內的模擬，首先需要基于SPH方法對拉格朗日流體控制方程［24］求解，建立水流的物理模型；其次通過邊界碰撞檢測方法，將水流保存至容器內。

2.1 SPH水流物理模型

利用SPH將連續(xù)的流體離散為具有一定體積有限數(shù)量的粒子，并采用核函數(shù)逼近的方式將場函數(shù)的積分形式表示成離散粒子求和的形式，SPH核函數(shù)逼近的表達式為：

式中：Wh為核函數(shù)［24］，h為光滑長度，決定了Wh的影響域。

離散后的粒子i的某物理量A可以由光滑半徑內的所有粒子對應的物理量A加權近似得到，得到式（1）的離散形式為：

式中：下標“j”用于描述粒子“i”核半徑內的所有粒子，r ij=r ir j，r i、r j表示粒子i、j的位置，mj、ρj表示周圍粒子的質量和密度。

在模擬水流時，需滿足水流的動態(tài)性，水粒子的狀態(tài)（位置、速度等）會隨著時間的推移變化。由于SPH方法是一種純拉格朗日方法，粒子的運動遵循文獻［24］的拉格朗日形式的流體控制方程：

其中：g表示粒子外力加速度；p表示粒子受到的壓力；u、ρ分別表示粒子的速度、密度；μ為流體的粘性系數(shù)。

依據(jù)式（3）的流體控制方程與式（2）插值核函數(shù)建立粒子間的約束關系，分析粒子受力，將水流粒子受力分解為壓力、粘滯力和外力：

通過式（4）得到水粒子i所受到的合力，結合牛頓第二運動定律推導出粒子i的加速度如式（5）所示。粒子在該加速度下不斷地更新位置和速度，實時模擬水流的運動。根據(jù)上述計算進行編程，其程序流程如圖2所示。

式（4）壓力項中的壓強計算公式?jīng)Q定了粒子的壓縮性，為最大限度實現(xiàn)流體的不可壓縮性。本文采用基于Tait方程的壓強計算方法計算水流粒子壓強，定義如下：

其中：ρ0為水流的初始密度；cs為水中聲速；τ為常數(shù)，本文取值為6。使用式（6）計算的壓強不會使水流密度產(chǎn)生較大的波動，保證了水流的弱可壓縮性。

圖2 SPH程序流程Fig.2 Flowchart of SPH program

2.2 粒子與容器邊界的碰撞處理

水流的加熱過程一般是在有限的局部空間中，這就需要判斷流體粒子在運動過程中是否與容器邊界發(fā)生碰撞，進而進行碰撞處理，將粒子運動限定在容器內。本文的碰撞檢測與處理方法與文獻［25］類似但稍有不同與改進，采用基于半邊數(shù)據(jù)結構網(wǎng)格的邊界處理方法，可以對復雜曲面邊界快速檢測，具體方法如下：

1）系統(tǒng)運行時對容器模型網(wǎng)格預處理，以半邊數(shù)據(jù)結構存儲網(wǎng)格以及網(wǎng)格的法向方向，并將采樣信息保存在網(wǎng)格中。

2）每個粒子從自身位置朝著運動方向射線檢測得到粒子與邊界的預碰撞點，由預碰撞點信息得到對應半邊網(wǎng)格中的碰撞點法向量。

3）將粒子位置沿著法線方向投影，計算粒子與邊界的距離，并判斷該值是否小于邊界碰撞閾值。小于則表示正在發(fā)生碰撞，需要修正粒子在下一時刻的位置和速度。

4）粒子與邊界碰撞后發(fā)生反彈運動（將粒子和邊界均看作理想剛體）。位置上將粒子沿著半邊網(wǎng)格法線方向投影到對稱表面；速度上分別修改碰撞后法線方向與切線方向的速度，合并后得粒子碰撞后下一時刻的速度：

式（7）中：ζ為粒子碰撞邊界后法向方向的彈性系數(shù)，ψ為粒子碰撞邊界后切向方向的摩擦系數(shù)。

利用上述原理構建加熱仿真中水流的物理模型，其相關參數(shù)如表1所示，并以太陽能熱水器的簡易保熱墻為容器，注入水粒子，效果如圖3所示。

表1 水流模型參數(shù)Tab.1 Parameters of water flow model

圖3 充滿粒子的保熱墻Fig.3 Thermal wall filled with particles

3 交互式水流加熱

以太陽能熱水器加熱工作為例分析水流的加熱過程：首先，保熱墻的集熱管受到太陽輻射吸收熱量進行電熱轉換，集熱管內的電?。ㄎ挥诩療峁軆缺谏蟼龋┓懦龃罅康臒?。然后，將電弧所放熱量按照一定的熱傳導方式傳導至水流，使保熱墻內各個位置水流產(chǎn)生溫差從而形成水流的熱運動。針對太陽能熱水器自身特點［26］，本文對熱水器加熱過程進行部分簡化，采用第一類邊界條件的單壁穩(wěn)態(tài)熱傳導模型，得到一個滿足實時性條件的溫度場計算方法。在給定的溫度場下傳導熱量計算水粒子的溫度變化，對不同溫度的粒子施加外力，進而模擬水流的熱運動。最后，設置可修改的仿真參數(shù)接口，使用戶自主的調節(jié)加熱條件，實現(xiàn)對水流的可交互加熱。

3.1 溫度場計算

由于集熱管電弧為水流的外熱源，在計算溫度場過程中，可將集熱管內的熱傳導過程視為不隨時間而變化穩(wěn)態(tài)熱傳遞狀態(tài)。通過對連續(xù)時間離散化，計算每個時間點上集熱管內的溫度場，從而得到單位時間步長上流體粒子溫度變化。當集熱管內流體處于穩(wěn)態(tài)導熱狀態(tài)且無內熱源時，其導熱方程為：

其中：T為溫度；ρ為流體密度；c為流體熱容；x、y、z為空間三維坐標軸方向；λ為熱傳導系數(shù)。

在加熱過程中，將集熱管電弧對管內水流的加熱簡化為第一類邊界條件下單壁的穩(wěn)態(tài)熱傳導問題，并且不考慮水體溫度對其溫度場的影響，則管內溫度場為：

其中：δ為集熱管內直徑；t1、t0分別為集熱管電弧發(fā)熱溫度和對側最低溫度。

3.2 粒子熱運動模擬

根據(jù)水流加熱的真實狀態(tài)可知，從電弧放熱傳導至水流的熱量遠大于水流內的熱量損失。因此，可以在實時的加熱仿真中忽略掉熱量損失的情況，并考慮加熱過程中水流的吸熱效率恒定，以簡化熱傳導模型的復雜度，增加仿真系統(tǒng)的實時性。將熱傳導問題轉化為線性方程求解，單位時間內r位置處的水流粒子吸熱量如式（10）所示，在該溫度場下，二維流體粒子加熱效果如圖4所示，流體粒子有規(guī)律地均勻加熱。

其中：α為傳熱速率系數(shù)，c為水的比熱容。

圖4 不同時刻粒子的加熱溫度Fig.4 Heatingtemperaturesof particlesat different times

由于水流受熱不均產(chǎn)生溫度差，即為水粒子間的內能的差異，而內能差異影響水粒子的運動狀態(tài)?；赟PH方法分析粒子熱運動，對不同內能的粒子施加外力，該力由光滑核半徑內鄰域粒子內能差值在核函數(shù)上的累加決定，表達式為：

其中：σ為熱力系數(shù)，可調節(jié)內能對外力的影響；Qi、Qj分別為粒子i、j的熱能。

伴隨著水流的運動，水流內部存在對流換熱的情況，即內部的粒子在運動過程中相接觸，粒子間表面存在熱量傳遞現(xiàn)象。不考慮對流換熱過程的熱量損失，熱流側放熱量與冷流側吸熱量相同，冷熱流側熱量變化量為：

式中：β為對流給熱系數(shù)；S為粒子接觸面積；tc、th分別為冷、熱流側粒子的溫度。

對流換熱后的粒子更新內能與溫度，并重新計算粒子熱運動。將上述循環(huán)添加至SPH程序流程中，形成水流熱運動主循環(huán)，流程如圖5所示。在該循環(huán)下二維流體粒子加熱效果如圖6所示?？梢钥闯?，伴隨著粒子熱運動與對流給熱過程，水粒子整體上呈現(xiàn)出高溫水粒子上運動，低溫水粒子向下運動，冷熱粒子循環(huán)的效果。且高溫粒子在熱運動的過程中將熱量傳導至鄰近的低溫粒子，使粒子群的溫度分布具有一定的隨機性，并且加速了粒子群整體的溫度變化。

圖5 粒子熱運動主循環(huán)Fig.5 Main circleof particlethermal motion

圖6 熱運動下粒子的加熱效果Fig.6 Heatingeffect of particlesunder thermal motion

最后，通過定義可編輯的參數(shù)接口，如入水量、注水速度、目標溫度等，以用戶輸入的參數(shù)設定水粒子初始屬性，如粒子數(shù)量、初始溫度等，使用戶自主修改加熱條件，以人機交互的方式實現(xiàn)多種條件下水流的加熱過程。

4 實驗結果與分析

本文仿真實驗在VS2017和Unity2018平臺上進行測試。以太陽能熱水器的加熱仿真為例，在該平臺對虛擬場景進行構建與渲染，并將可控參數(shù)在UI（User Interface）中顯示，效果如圖7所示。

圖7 場景渲染及UI設計Fig.7 Scenerenderingand UIdesign

圖7中的半透明模型為熱水器保熱墻的簡易模型，使用該模型近似替代原模型進行仿真，可以降低模型復雜度，增加系統(tǒng)的實時性。并且該仿真系統(tǒng)所用粒子數(shù)最多為1 500，幀率保持在30 frame/s，滿足了系統(tǒng)高效性要求。

根據(jù)上述的建模方法對熱水器中水流加熱過程進行模擬。利用光滑粒子流體動力學方法以及拉格朗日流體控制方程建立水流的物理模型，實例化水粒子并通過碰撞檢測方法將水粒子保存在簡易保熱墻中。當注水達到指定容量后用戶可以自主選擇開始加熱。加熱后電弧放熱形成溫度場，粒子實時吸熱更新溫度，吸熱的效率由式（10）的傳熱速率系數(shù)α決定。為體現(xiàn)太陽能熱水器與陽光照射的關系，本文通過陽光直射角度的變化修改α的值。在加熱過程中，由于水粒子的位置不同，粒子的吸熱量不同，使得不同位置的水粒子產(chǎn)生內能差，從而改變水粒子的運動狀態(tài)，形成水流熱運動。

水流加熱的速率不僅與傳熱速率系數(shù)α相關，還受式（12）粒子間熱對流給熱系數(shù)β影響。表2給出了初始注水溫度為20℃時，β分別取值0.3和0.5時，不同水量加熱到指定溫度的耗時表?？梢钥闯?，總體上加熱耗時與注水量以及目標溫度的大小成正相關，增加β可以加速水粒子間的熱量傳遞，使容器內水流的整體溫度較快的提升。

表2 不同條件下加熱耗時 單位：sTab.2 Heatingtimesunder different conditions unit：s

為了直觀地觀察水加熱過程的溫度變化，使用淺色表示注入的冷水，深色表示加熱后的熱水，以深淺插值變化表示水溫的變化，效果如圖8所示。

圖8 水加熱過程的顏色變化Fig.8 Changeof color duringwater heating

在熱水器的加熱仿真中，用戶可以通過修改注水溫度、注水量、目標溫度等可控參數(shù)，從而改變水粒子的實例化數(shù)量和相應屬性，實現(xiàn)交互式控制水流的加熱過程。

在流體加熱模擬方法上，文獻［17］所述的神經(jīng)網(wǎng)絡方法雖然精度較高、算法計算量較小但加熱過程無流體運動感，且需要對大量樣本數(shù)據(jù)構建神經(jīng)網(wǎng)絡模型。本文在保持系統(tǒng)高效性與交互性的前提下，通過利用光滑粒子流體動力學方法可對水流加熱進行有效模擬，具有較強的流動感，并且可以通過溫度場計算方法與相關的參數(shù)值調節(jié)水流加熱速率，加熱時間接近合理化。

5 結語

在現(xiàn)有太陽能熱水器加熱原理的基礎上，本文提出了一種基于光滑粒子流體動力學的流體熱運動仿真方法。該方法可以真實地模擬熱水器保熱墻內的注水、與水的流動過程，以及加熱后水流內部的熱運動現(xiàn)象。在保證系統(tǒng)實時性與高效性條件下，能夠以人機交互的方式實現(xiàn)多種條件下熱水器中水流的加熱仿真，以直觀的仿真動畫與系統(tǒng)參數(shù)使用戶在虛擬場景中直接地了解產(chǎn)品信息。本文的研究成果對于多種虛擬場景中交互式水流加熱仿真具有實際的參考意義。在現(xiàn)有的研究基礎上，如何針對水加熱較為真實的模擬沸騰現(xiàn)象以及水汽化后氣-液兩相流的仿真建模將是未來的重點。