,

(山東科技大學(xué) 電子通信與物理學(xué)院，山東 青島266590)

礦井突水災(zāi)害是煤礦的主要災(zāi)害之一[1]，礦井突水造成淹井和傷亡事故，給國家、社會和家庭帶來嚴重的直接或間接的經(jīng)濟損失和人身損害。因此，礦井突水蔓延研究對礦井治水安全具有重要意義。張麗娟[2]研究了基于OSG的礦井突水應(yīng)急虛擬仿真系統(tǒng)關(guān)鍵技術(shù)，改進MC和光滑流體動力學(xué)(smoothed particle hydrodynamics，SPH)算法，基于表面紋理和粒子系統(tǒng)完成了水災(zāi)的真實感表達。劉欣倩[3]研究了基于Unity 3D粒子系統(tǒng)進行礦井突水事故的虛擬仿真，展示了突水淹沒巷道的情景。李長春[4]進行了巷道建模和水淹模擬,宏觀展現(xiàn)水淹巷道情景，但缺乏水災(zāi)蔓延動態(tài)研究。李翠平等[5]結(jié)合水災(zāi)流動的水力特征，建立了能夠真實模擬水流蔓延三維動態(tài)仿真模型，但只是宏觀體現(xiàn)。汪金花等[6]、王鵬[7]進行了水災(zāi)仿真實驗，主要研究的是水災(zāi)數(shù)學(xué)模型和最優(yōu)避難路徑。趙怡晴等[8]構(gòu)建了礦井突水仿真的一體化模型，宏觀展示了水災(zāi)淹沒范圍。

對于虛擬環(huán)境中礦井突水災(zāi)害逃生培訓(xùn)來說，巷道內(nèi)突水蔓延的動態(tài)模擬能夠為安全訓(xùn)練提供較好的輔助作用，因此對水體本身進行具有物理規(guī)律的真實感模擬是突水蔓延模擬的需求?；谖锢淼牧黧w模擬方法分為基于網(wǎng)格的歐拉法和基于粒子的拉格朗日法。用歐拉網(wǎng)格法模擬流體的計算效率較高，Stam[9]提出基于網(wǎng)格法的流體模擬，提高了流體的穩(wěn)定性和計算效率。邵緒強等[10]提出一種基于shallow water方程的物理模擬方法，加速大規(guī)模流體的物理模擬計算。1983年，Reeves[11]首次提出了粒子系統(tǒng)的方法并廣泛應(yīng)用到虛擬仿真中。Hu等[12]基于SPH算法在游戲引擎使用粒子系統(tǒng)進行了流體的模擬?；诹Ｗ拥睦窭嗜辗ㄌ貏e是SPH粒子法能夠避免歐拉方法網(wǎng)格扭曲的缺陷而且能夠處理流體自由邊界問題，適合求解高速碰撞等問題[13-14]。

本研究采用SPH算法求解N-S方程獲得流體的物理模型并應(yīng)用到粒子系統(tǒng)，采用固定邊界粒子法進行邊界處理，結(jié)合三維八叉樹方法進行表面重構(gòu)實現(xiàn)水體的渲染，在Unreal Engine 4中實現(xiàn)了水流體的運動仿真效果，并在三維巷道地形中實現(xiàn)了突水水平流動、上向升漲和下向蔓延的仿真。使得沉浸性和交互性的礦井水災(zāi)逃生培訓(xùn)和災(zāi)情預(yù)測成為可能，在礦井突水災(zāi)害逃生演習(xí)中起到重要作用。

1 SPH數(shù)值解算

光滑粒子流體動力學(xué)，是一種基于物理的純拉格朗日粒子算法[15]。在SPH方法中，將流體視為質(zhì)點系，系統(tǒng)的狀態(tài)用包含各種物理量(密度、壓力、速度等)的質(zhì)點來描述，通過求解質(zhì)點組的動力學(xué)方程和跟蹤每個質(zhì)點的運動軌道，求得整個系統(tǒng)的力學(xué)行為。SPH的公式構(gòu)造不受粒子分布的隨意性影響，可以很自然地處理一些具有極大變形的模擬情況[16]，適用于礦井水災(zāi)蔓延模擬的研究。

基于粒子的流體模擬中，每個粒子的運動都遵循牛頓第二定律：F=ma。在SPH方法中，流體粒子的質(zhì)量取決于粒子的密度，所以一般用密度來代替質(zhì)量。作用在粒子上的力由外力、壓力和粘性力組成。因此水作為一種不可壓縮粘性流體，其模擬可用以下N-S方程的簡化形式描述：

ρa=-p+ρg+μ2u。

(1)

其中：ρ為液體的密度，p為壓力，u為速度，為梯度，μ為粘性系數(shù)，2為拉普拉斯算子。忽略其他相關(guān)較小的力，外部力一般指重力，外力F=ρg。μ2u表示粘力項。

假設(shè)流體中一個位置為ri的點，此處的密度為ρ(ri)、壓力為p(ri)、速度為u(ri)，那么根據(jù)公式(1)，可以計算出此處粒子的加速度

(2)

采用跳蛙(leap-frog)算法對加速度進行時間積分后，得到粒子在下一時刻的速度與位置[17]。在半個積分時間步得到速度，并利用這一速度計算新的位置，其位置和速度表達式為：

(3)

(4)

SPH方法中，整個系統(tǒng)是由具有獨立質(zhì)量、占有獨立空間的有限個粒子表示的?？梢杂梢韵铝Ｗ咏品ǖ玫?。假設(shè)流體中某點r(此處不一定有粒子)，在光滑核半徑h范圍內(nèi)有數(shù)個粒子，則該處場量A的計算可近似為：

(5)

其中，Ai指的是要累加的某種場量，本研究的場量是外力、密度、粘度；r指的是該粒子的當(dāng)前位置，ri為i處的粒子位置；h指的是光滑核半徑；函數(shù)W是光滑核函數(shù)。

圖1 光滑函數(shù)影響閾內(nèi)的粒子示意圖Fig.1 Schematic diagram of the particle in the threshold of the influence of smooth function

本研究在SPH方法的基礎(chǔ)上進行流體模擬，通過SPH算法將流體連續(xù)方程轉(zhuǎn)化為支持域內(nèi)粒子求和的離散方程，其壓力項和粘力項的求解式分別表示為：

(6)

(7)

對于密度項選取的光滑核函數(shù)為：

(8)

對于壓力項選取的光滑核函數(shù)為：

(9)

對于粘力項所用的光滑核函數(shù)選取

(10)

由式(2)，式(5)～(10)即可求出粒子的加速度，進而根據(jù)式(3)、(4)獲取粒子下一時刻的速度與位置。

2 邊界處理

基于SPH的粒子法是純拉格朗日粒子方法，粒子與邊界的作用力必須轉(zhuǎn)化為粒子與粒子之間的相互作用才可以進行計算。本研究采用一種精確并節(jié)省時間的固定邊界粒子法[20]，將粒子與邊界的作用轉(zhuǎn)化為粒子與粒子之間的作用，得以應(yīng)用SPH方法計算邊界問題。

如圖2所示，固定邊界粒子法在邊界上布置固定粒子，用邊界粒子與流體粒子反應(yīng)代替固體邊界與流體粒子發(fā)生反應(yīng)，通過內(nèi)部粒子的壓力，估計邊界粒子的壓力。

圖2 固定邊界粒子法Fig.2 Fixed boundary particle method

選取以光滑長度h為半徑的區(qū)域所有內(nèi)部流體粒子，然后計算內(nèi)部流體粒子與邊界粒子的間距rij，再通過式(11)近似求得此處的壓力

(11)

式中，rij表示邊界粒子i支持域內(nèi)的粒子j到邊界粒子i的距離，Pi表示粒子i的壓力，h表示光滑長度。通過公式可以看出，距離邊界點較遠的粒子，h-rij較小，其對邊界點的壓力貢獻值較小；反之，距離邊界點越近，h-rij較大，對Pi的影響越大。求得Pi后，在求解i粒子的受力時，Pi會對i產(chǎn)生影響，通過這種方式保證了粒子不會超出固壁邊界，符合正常的物理規(guī)律，模擬出的粒子邊界效果更為真實。計算過程中，邊界粒子和內(nèi)部粒子一樣，搜索臨近粒子，計算壓力和密度，不同的是邊界處的粒子速度始終為零，即粒子不發(fā)生運動。

3 流體渲染

流體效果的渲染基于八叉樹的自適應(yīng)流體的表面重建方法[21]，跟蹤流體表面的粒子，構(gòu)建自適應(yīng)距離場。對于任意流體的粒子i,計算一個再歸一化矩陣Bi[22]:

Bi=[∑jW(ri)?(rj-ri)Vj]-1。

(12)

圖3 渲染流程Fig.3 Rendering process

(13)

然后根據(jù)八叉樹結(jié)點對不同粒子間距采樣，構(gòu)建自適應(yīng)距離場，重建出僅由一層粒子構(gòu)成的薄膜結(jié)構(gòu)的網(wǎng)格表面用于流體表面繪制。再結(jié)合Unreal Engine 4的shader著色程序進行最終的渲染。

Unreal Engine 4中，紋理(Texture)是基于GPU上需要運用特定的數(shù)據(jù)結(jié)構(gòu)，紋理的坐標、紋理的查詢分別相當(dāng)于數(shù)組中元素的索引和元素的讀取。

首先，把距離場的FDistanceFieldVolumeTexture數(shù)據(jù)對應(yīng)提交到GPU中的3D紋理DistanceFieldTexture中；隨后，通過調(diào)用函數(shù)使流體的距離場數(shù)據(jù)與全局距離場的網(wǎng)格相對應(yīng)，此時需要變換矩陣使GPU距離場的數(shù)據(jù)與世界空間的數(shù)據(jù)相互轉(zhuǎn)化；最后針對每一個視口網(wǎng)格生成一個全局距離場大小相同、密度不同的3DTexture，并在此網(wǎng)格進行更新。

以上步驟完成后，需要將所有數(shù)據(jù)渲染到當(dāng)前相機的平面，呈現(xiàn)在屏幕。根據(jù)攝像機坐標系下的坐標值，可以求得在當(dāng)前攝像機下最近點的三維坐標，將屏幕作為投影平面，相機的中點作為起點，向模型距離場發(fā)出的采樣射線與場景中的流體表面相交，計算出此交點的屬性，就獲取流體距離場中每個網(wǎng)格的屬性數(shù)據(jù)，紋理采樣然后進行繪制。繪制之前啟用深度緩沖(Depth Buffer)對比采樣點的距離生成正確的深度感知效果。最后調(diào)用一個DrawRectangle把距離場信息貼到對應(yīng)的View全屏上。

流體表面繪制自定義函數(shù)如下：

void WaterShaderOnSurface(float4 surfPos)

{

float3 surfNorm = CalcSurfaceGrad(surfPos.xyz);

float3 WorldPosition = GetWorldPosition(surfPos.xyz);

float FluidDist =(surfPos.w - lastPos.w) * WorldRayDirLen;

float3 WorldNormal = normalize(surfNorm);

WaterShader.OnEnterSurface(lastPos.xyz+surfPos.xyz)*0.5f, FluidDist);

}

4 水災(zāi)仿真實驗

在進行礦井突水蔓延仿真實驗時選擇的驗證平臺為：硬件平臺為Intel Core i7-7700 CPU @3.60 GHz處理器、16 GB內(nèi)存、NVIDIA Titan XP顯卡；軟件平臺為Windows 10操作系統(tǒng)，仿真平臺為Unreal Engine4和Visual Studio 2015。

在Unreal Engine 4源代碼中定義解算相關(guān)函數(shù)，通過解算獲得粒子的速度、位置等參數(shù)；然后在Unreal Engine 4編輯器中創(chuàng)建粒子Particle藍圖類，修改粒子的數(shù)量、周期和顏色等參數(shù)。

在巷道中，流體的邊界主要是與巷道墻壁和設(shè)備的交界面，把巷道壁和巷道內(nèi)設(shè)備作為剛體處理，從Unreal Engine 4距離場數(shù)據(jù)的3D紋理中獲取邊界信息后，在邊界上布置一層固定粒子，利用粒子與邊界粒子間的作用計算粒子在邊界處與墻體碰撞后的流動狀態(tài)。

生成固定邊界粒子的部分函數(shù)如下：

Void BuildSBoundaryParticle(uint3DispatchThreadId:SV_DispatchThreadID)

{uint3 ijk = GetActiveIndex(DispatchThreadId);

float DistanceToNearestSurfaceParticle;

}

計算粒子下一刻的速度函數(shù)：

Velocity.xyz += DiffuseParticles.DeltaSeconds * acc.xyz /2* 100

得到粒子更新的速度，相應(yīng)位置確定為：

Position.xyz += DiffuseParticles.DeltaSeconds * Velocity.xyz。

4.1 突水蔓延效果展示

在Unreal Engine 4中用畫刷工具繪制實驗所需的礦井巷道三維模型，進行突水蔓延模擬的巷道模型部分結(jié)構(gòu)如下：

1) 將最終構(gòu)建完善的粒子系統(tǒng)作為一個Actor加入到巷道模型中的突水位置；

2) 主要參數(shù)值設(shè)置為：

Spawn.Rate.Constanat=40 000.0;

圖4 巷道部分結(jié)構(gòu)圖

Spawn.Rate Scale.Constanat=2.0;

LifeTime.Constant=1 000.0;

Start Velocity =(100.0,0.0,0.0)。

實驗中，最先發(fā)射出的粒子在落到巷道表面時重力與邊界粒子的支持力相互抵消，只有后面的粒子對其施加一個壓力產(chǎn)生加速度向前，其他粒子由壓力、重力和粘性力共同產(chǎn)生加速度，由SPH水模型的求解流程獲得粒子下一刻的速度與位置，表現(xiàn)出水流向前蔓延的效果。

突水下向蔓延過程中，重力、粘力、壓力的合力提供給水流沿巷道地表面的加速度，表現(xiàn)出水流沿斜面向下的流動狀態(tài)。突水上向蔓延過程中，水流持續(xù)蔓延至巷道低洼處，每個粒子占據(jù)一定的空間體積，底部粒子不斷地往上層累積，水面不斷抬高。

當(dāng)突水蔓延過程發(fā)生在巷道左側(cè)分叉口，由于左側(cè)墻壁不再有粒子對流體產(chǎn)生壓力，左側(cè)流體表面粒子只受到內(nèi)部相鄰流體粒子向左的壓力，上部的流體粒子受自身重力和相鄰內(nèi)部流體粒子的壓力，部分粒子會產(chǎn)生向左的加速度，因此產(chǎn)生流體在分叉口分流的現(xiàn)象。

最終部分巷道內(nèi)突水蔓延仿真效果如圖5所示，其中右側(cè)黑色+標記處為突水點。通過仿真的實驗結(jié)果可以看到，在突水發(fā)生后，水流首先涌入低洼處平坦地勢的巷道進行下向蔓延，在水流灌滿低洼處巷道后水位開始上漲，隨后水流會沿地勢由低到高灌滿整個巷道，并在分叉處分流。

圖5 巷道突水蔓延效果Fig.5 Simulation results of water inrush in roadway

本研究的方法在Unreal Engine 4虛擬環(huán)境中實現(xiàn)了突水蔓延的狀態(tài)，巷道內(nèi)部蔓延效果圖如圖6所示，在具有真實感的突水仿真的同時，實現(xiàn)蔓延的動態(tài)變化趨勢。效果直觀可見，相對于其他路徑算法研究的突水蔓延，具有客觀的實用性。

圖6 巷道內(nèi)部蔓延效果圖Fig.6 Effect of internal spread of roadway

5 結(jié)論

運用SPH算法求解粒子的壓力項、密度項和粘力項得到粒子的加速度最終求得粒子的速度與位置，最終在復(fù)雜地形的巷道中實現(xiàn)了具有真實物理規(guī)律的突水蔓延仿真，包括水流上向蔓延、水流下向蔓延和水流在交叉口處的蔓延狀態(tài)。通過實驗，所研究的突水蔓延可以應(yīng)用于礦井水災(zāi)逃生模擬中，但本研究的突水蔓延仿真局限于較小的巷道范圍，下一步需要對大型巷道蔓延進行研究。