,

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

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

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

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

1 SPH數(shù)值解算

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

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

ρa(bǔ)=-p+ρg+μ2u。

(1)

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

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

(2)

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

(3)

(4)

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

(5)

其中，Ai指的是要累加的某種場(chǎng)量，本研究的場(chǎng)量是外力、密度、粘度；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ǔ)上進(jìn)行流體模擬，通過SPH算法將流體連續(xù)方程轉(zhuǎn)化為支持域內(nèi)粒子求和的離散方程，其壓力項(xiàng)和粘力項(xiàng)的求解式分別表示為：

(6)

(7)

對(duì)于密度項(xiàng)選取的光滑核函數(shù)為：

(8)

對(duì)于壓力項(xiàng)選取的光滑核函數(shù)為：

(9)

對(duì)于粘力項(xiàng)所用的光滑核函數(shù)選取

(10)

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

2 邊界處理

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

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

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

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

(11)

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

3 流體渲染

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

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

(12)

圖3 渲染流程Fig.3 Rendering process

(13)

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

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

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

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

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

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

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

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

Void BuildSBoundaryParticle(uint3DispatchThreadId:SV_DispatchThreadID)

{uint3 ijk = GetActiveIndex(DispatchThreadId);

float DistanceToNearestSurfaceParticle;

}

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

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

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

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

4.1 突水蔓延效果展示

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

1) 將最終構(gòu)建完善的粒子系統(tǒng)作為一個(gè)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)。

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

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

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

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

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

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

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

5 結(jié)論

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