王建華,冉煜琨

(成都理工大學工程技術(shù)學院,樂山 614000)

1 引言

我國有著遼闊的海岸線和大量島嶼,雖然我國很多沿海城市沒有處在板塊交界處,發(fā)生海嘯的可能性不像鄰國日本那樣多.但近幾年,較多的海底大地震給我國臺灣島和福建沿海城市造成海嘯的可能性變大.針對可能發(fā)生海嘯的安全隱患問題,利用一定的分析評估方法是防災(zāi)減災(zāi)的重要手段之一.

目前在實踐層面上,廣泛使用的海嘯分析是基于淺水長波理論[1]的數(shù)值分析,該理論中假設(shè)水粒子的垂直加速度非常小,取垂直方向上流速的均值.這種分析方法對二維現(xiàn)象進行了模擬,可用于理解從地震中心至沿海城市廣闊區(qū)域上的海嘯傳播.然而,該方法難以確定3D 屬性在其中占主導(dǎo)地位的局部水流的自由表面形狀[2].一般情況下,海嘯預(yù)測的內(nèi)容是海嘯到達地面后,隨時間發(fā)展,海水覆蓋的區(qū)域變化情況.此外,由于水波的三維性會造成的影響,建筑物附近的水流預(yù)測也是預(yù)測內(nèi)容之一,但前者通常為預(yù)測的主要內(nèi)容.

GIS 數(shù)據(jù)以其在時空數(shù)據(jù)的存儲表達分析等方面的優(yōu)勢,使得在多種3D 仿真和可視化領(lǐng)域得到廣泛關(guān)注,如潰壩洪水演進[3]、洪水淹沒分析[4]等.傳統(tǒng)GIS方法大多采用數(shù)字高程模型(Digital Elevation Model,DEM)[5],利用有源淹沒或者無源淹沒進行分析,其缺點是忽略了物理學和動力學的因素,時效性和可靠性較弱.

在3D 流動分析中,當水波發(fā)生顯著變化時,網(wǎng)格的數(shù)值分析方法[6]也經(jīng)常被使用,如有限差分法和有限元法[7]等.在眾多理論分析方法中,光滑粒子流體動力學(Smoothed Particle Hydrodynamics,SPH)[8]逐漸得到認可.SPH是一種無網(wǎng)格的拉格朗日粒子法,其基本思想是將連續(xù)的流體用相互作用的粒子組來描述,各粒子上承載各種物理量,包括質(zhì)量、速度、加速度等.通過求解粒子組的動力學方程和跟蹤每個粒子的運動軌跡,求得整個系統(tǒng)的力學特征[8].SPH在海嘯3D 仿真[9]、洪水災(zāi)難評估[3,10]等方面得到了應(yīng)用,驗證了SPH具有一定的可行性和準確性.也有研究者將SPH 擴展到不可壓縮的光滑粒子流體動力學(Incompressible Smoothed Particle Hydrodynamics,ISPH)[11].

本文在SPH 方法的基礎(chǔ)上,對邊界問題和緩解不可壓縮問題進行描述和解決,建立能夠?qū)[爆發(fā)進行有效分析的工具.提出了一系列的預(yù)處理程序和地理分析模型.該模型能夠通過從一個地理信息系統(tǒng)(GIS)中,得到的3D 位置信息(SHP)和數(shù)字高程模型(DEM)來展示地形、海拔和建筑物的外部形態(tài).

2 提出的分析方法

2.1 SPH 方法的基本方程

SPH 方法的基本概念是將位于空間點x和時間t的函數(shù) φ (x,t)作為一個積分表達式進行逼近,其定義如下:

式中,W為核函數(shù),即一種加權(quán)函數(shù).在SPH 方法中,一般將光滑長度在h內(nèi),具有非零正值,且滿足統(tǒng)一條件的緊支撐函數(shù)作為一個近似.基于離散化,使用分布在空間中的粒子 xi對式(1)進行逼近,并在長度h內(nèi)使用其附近的粒子值獲得加權(quán)總和:

式中,下標i和j表示粒子數(shù)量,ρj和mj分別表示與粒子數(shù)量j相關(guān)的平均密度和代表性質(zhì)量.rij(=|rij|)表示粒子距離,rij(=xi?xj)表示粒子相對位置向量.符號<·>為SPH 方法鄰近粒子數(shù)值的近似值.

2.2 流體分析的控制方程

為求解不可壓縮的流體問題,應(yīng)該在滿足質(zhì)量守恒和動量守恒定律的前提下,得出兩個自變量,即流速u和壓力p.控制方程的拉格朗日形式定義如下:

式中,v為動力粘性系數(shù),g為重力加速度.

假設(shè)水的密度是恒定的,基于這一假設(shè),可以將質(zhì)量守恒定律(式(3))改寫為:

將式(4)和式(5)應(yīng)用到粒子i,得到:

2.3 SPH 方法中緩解不可壓縮條件

ISPH 方法[11],使用投影法對運動方程的速度場和壓力場進行分離,以達到對壓力場的隱式計算,和速度場的顯式計算.另外,使用投影法,通過定義一個暫定狀態(tài),可以對速度場和壓力場進行分離.

下面將投影法應(yīng)用到SPH 方法中,從n至n+1 對變量進行更新.首先,對式(4)中的時間導(dǎo)數(shù)項進行前向差分近似,并將中間速度u 定義在一個中間狀態(tài),以對速度進行分離:

在分離后的加速度分量之間,計算中間速度u 如下:

這里假定式(9)右邊的第一項和第二項分別對應(yīng)于式(4)的壓力梯度項和其他項.然后,對從中間狀態(tài)至下一個時間步的速度進行更新:

如上所述,通過實施兩個單獨過程,以對速度的狀態(tài)進行更新.而壓力是通過求解壓力Poisson 公式獲得:

為了緩解不可壓縮條件,作如下定義:

在SPH 方法中,從粒子的分布中對密度進行數(shù)值計算.在分析過程過程中很難始終滿足恒定密度條件,因為僅在鄰近粒子為固定數(shù)量且粒子保持完全的均勻分布時,密度才是恒定的.因此,有必要采用一個折中方案,即:密度長期沒有變化,但允許密度在瞬間存在一定程度的誤差.

式(13)給出的壓力Poisson 方程完全服從無散度條件下的模型(松弛參數(shù)為零).此外,當瞬間密度與初始密度吻合(或密度小到可以忽略不計),則壓力Poisson方程和模型可被視為相同,因為壓力Poisson 方程源項的第2 項可以忽略不計.根據(jù)該公式,通過密度差項逐步消去分析過程中產(chǎn)生的累積誤差,使得即使長期計算,密度也幾乎保持恒定,從而得到了較好的體積守恒方案.

2.4 利用虛擬標記進行邊界處理

在粒子法的邊界處理中,一般會假設(shè)一個分析模型,該模型的墻壁邊界基本符合物理邊界.然而,對于曲線或梯度邊界,很難得到粒子在邊界內(nèi)部的均勻分布.因此,在粒子模型的實際開發(fā)中,粒子通常被放置在格子(將CAD 數(shù)據(jù)劃分為網(wǎng)格結(jié)構(gòu))的中心點或交點處.通過這一簡單預(yù)處理,實現(xiàn)了粒子在一個區(qū)域的均勻分布;然而,粒子在邊界的分布呈階梯狀,這樣的分布不同于實際物理邊界中的平滑分布.由此導(dǎo)致邊界附近的流動不自然.

針對該問題,本文處理的方法是:將虛擬標記與墻壁粒子對稱地放置在墻壁邊界上.雖然虛擬標記與SPH 方法的計算沒有直接關(guān)系,但其測量點被用于對墻壁粒子施加適當?shù)倪吔鐥l件.這里,為滿足邊界表面上的滑移條件,假設(shè)墻壁粒子的速度為vw,則該速度與虛擬標記流動速度vv的對稱性映像為:

式中,M為二階張量,其利用一個內(nèi)向法線向量n=(n1,n2,n3)t和Kronecker δ執(zhí)行鏡像運算:

為滿足非滑移條件,計算與虛擬標記流速 vv滿足點對稱關(guān)系的墻壁邊界流速 v′w,并將其應(yīng)用到墻壁粒子.其中,v′w使用點對稱張量R計算如下:

接下來,確定虛擬標記的壓力,使之滿足墻壁表面不均勻壓力的Neumann 條件[12]即可:

式中,n為墻壁的內(nèi)向法線向量.

3 分析模型的開發(fā)

在對海嘯爆發(fā)進行建模時需要考慮到作為主要障礙物的建筑物.眾所周知,3D 分析的準確度很大程度上取決于分析模型的可靠性.為了盡可能地將模型海拔高度和建筑物形態(tài)與實際的地理位置相匹配,本文建立了一個利用GIS 數(shù)據(jù)進行粒子分析模型的程序.根據(jù)海拔上的DEM 數(shù)據(jù)和建筑物上的SHP 數(shù)據(jù),開發(fā)出粒子分析模型.DEM 數(shù)據(jù)是通過航測得到的高程數(shù)據(jù),并以等間隔的網(wǎng)格形式保存.DEM 數(shù)據(jù)剔除建筑物高度數(shù)據(jù),給出了地表高程.由此,DEM 數(shù)據(jù)一般用作高程數(shù)據(jù),而SHP 數(shù)據(jù)(GIS 地理空間數(shù)據(jù)文件的一種格式)一般用于定義建筑物的幾何形狀.在SHP數(shù)據(jù)中,每個建筑物的輪廓被定義為2D 數(shù)據(jù),其中還包括了建筑物的高度信息.當使用通用3D CAD 數(shù)據(jù)定義地表和建筑物的輪廓時,可以使用粒子分析的通用預(yù)處理程序(例如Meshman 服務(wù))將該3D CAD 數(shù)據(jù)轉(zhuǎn)換為粒子模型數(shù)據(jù).圖1給出了本文粒子模型的開發(fā)流程.

圖1 粒子模型的開發(fā)流程

下面將解釋使用上述兩類GIS 數(shù)據(jù),將地表輪廓定義為STL 數(shù)據(jù)(類似于三角面片的一種3D CAD 格式),以及STL 數(shù)據(jù)轉(zhuǎn)換為粒子數(shù)據(jù)的一系列程序.遵循的基本流程是:從DEM 數(shù)據(jù)中估計出每個建筑物位置的高程數(shù)據(jù),并使用SHP 數(shù)據(jù)基于建筑物的平面輪廓對建筑物的高度進行掃描,從而將高程數(shù)據(jù)轉(zhuǎn)換為3D CAD 數(shù)據(jù).由于地表通常包含斜坡;因此,在一個建筑物的底部和地表之間可能會產(chǎn)生間隙.如果在粒子模型中不進行任何修正,則粒子可能會進入該間隙.因此,當制備建筑物的CAD 數(shù)據(jù)時,將建筑物的底部降低一個恒定的量(5 m),如圖2所示.

圖2 降低STL 數(shù)據(jù)的建筑物底部

本文使用的粒子分析預(yù)處理程序能夠為CAD 數(shù)據(jù)中定義的多個表面,獨立地生成粒子.并在粒子的位置與另一個粒子的位置發(fā)生重疊時,按照優(yōu)先次序移除重復(fù)粒子,開發(fā)出具有均勻粒子分布的分析模型.由此,將定義地表的STL 數(shù)據(jù)和定義建筑物的數(shù)據(jù),分別存儲在單獨的文件中.接著,如圖3所示對建筑物和地面粒子進行制備,通過將地面粒子設(shè)為高優(yōu)先級以使得建筑物粒子和地面粒子在與地表相對應(yīng)的區(qū)域中發(fā)生重疊時,自動移除建筑物粒子.

圖3 地面粒子

4 分析與評價

下面將基于前文解釋的方法,開發(fā)一個分析模型,并使用SPH 方法對一個沿海城市海嘯爆發(fā)進行分析.

4.1 分析條件

直徑為2 米粒子的分析模型如圖4所示.在該分析中,假設(shè)海嘯的波形會急劇上升.通過放置在目標區(qū)域下的水流對河流進行模擬.假設(shè)河流的平均海拔低于海平面4 m,該河底部上為高度8 m的水體.

圖4 分析模型的示意圖

基于長波理論,該河的流速公式為c=波速為9 m/s.本文采用了滑移邊界條件.對粒子直徑分別為1 m、2 m和4 m的3個案例進行了仿真,并比較得出的結(jié)果.時間增量被固定為0.005 s,時間步的總數(shù)為30 000.數(shù)值分析的其他條件如表1所示.

表1 分析條件

4.2 結(jié)果評價

在海嘯流入開始后60 s (時間步為12 000)時,3個案例的結(jié)果如圖5(a)所示.海嘯粒子沒有滲透入建筑物,這說明了已成功開發(fā)出分析模型,在建筑物和地面之間沒有間隙.在對3個案例的分析結(jié)果進行比較時,本文發(fā)現(xiàn)被淹沒的陸地面積隨空間分辨率的增大而增加.在對粒子直徑1 m和2 m,以及2 m和4 m 案例被淹沒面積之間的差異進行比較時,前者間的差異要小于后者間的差異,這意味著數(shù)值分析的結(jié)果隨空間分辨率的增加而收斂.圖5(b)給出了圖5(a)中框內(nèi)區(qū)域的放大視圖,該區(qū)域中差別較為明顯.圖5(b)中定義的點A、B、C 處測量到水位升高情況.在不同分辨率模型之間的比較中,海嘯到達時間和水位升高歷史表現(xiàn)出相似趨勢.

由圖5,在使用直徑為1 m的粒子與使用直徑為2 m的粒子案例中觀察到了相同的趨勢;然而,使用直徑為2 m 粒子的案例中觀察到的趨勢,則未出現(xiàn)在粒子直徑為4 m的案例中.即使使用Meshman 從相同的STL 數(shù)據(jù)中生成粒子,但在不同分辨率下,地形表示和建筑物形狀均會發(fā)生變化.特別是當使用低分辨率時,其未能展示建筑物之間的小路,從而導(dǎo)致海嘯爆發(fā)中出現(xiàn)阻塞物.舉例來說,單車道的平均道路寬度為3.5 m,而使用直徑為4 m的粒子時,分析模型不會展示出此類道路.此外,低分辨率將導(dǎo)致地面梯度的再現(xiàn)性較低,從而使分析準確度較低.

圖5 海嘯爆發(fā)分析示意圖

圖5中沒有展現(xiàn)的是,在使用直徑1 m的粒子時,整個東側(cè)的海拔大約比西部海拔低1 m.然而,在使用直徑為2 m 或4 m的粒子時,均未觀察到這一趨勢.因此,使用1 m 粒子時,可能會造成高程重現(xiàn)性的下降.綜上,應(yīng)該采用粒子直徑為2 m的空間分辨率對沿海城市地區(qū)的海嘯爆發(fā)進行建模,該分辨率能夠重現(xiàn)主要道路,從而實現(xiàn)較高的準確度.圖6給出了海嘯開始流入后150 s (時間步30 000 處)時的被淹沒區(qū)域.

圖6 海嘯流入150 s 后被淹沒的區(qū)域

4.3 照片級的可視化

在傳統(tǒng)的可視化方法中,粒子一般被表示為原始形狀.使用這一可視化方法,可以對計數(shù)的類型進行切換 但照片缺乏真實感.

本文使用建筑領(lǐng)域可視化軟件Lumion 開發(fā)了3D 動畫.Lumion 提供了使用視差效果進行3D 動畫創(chuàng)建的功能,支持創(chuàng)造帶有真實感的3D 動畫,還可以導(dǎo)入通用CAD 數(shù)據(jù),可加入與景物、日照條件和陰影相關(guān)的數(shù)據(jù).因此,只有通過分析模型中虛擬城市的CAD數(shù)據(jù)和海嘯的CAD 數(shù)據(jù),才能實現(xiàn)3D 可視化.此外,在對海嘯的表面進行定義后,將建筑物的真實圖像加入建筑數(shù)據(jù),并將水流紋理加入海嘯數(shù)據(jù)中,以產(chǎn)生真實感,進行可視化的一系列程序如圖7所示.可視化的示例如圖8所示.

圖7 可視化方法的程序

圖8 3D 可視化照片的示例

4.4 與網(wǎng)格法的比較分析

網(wǎng)格法在諸多分析中經(jīng)常使用.在3D 流動分析中,當水波的波形變化顯著變化時,使用網(wǎng)格的數(shù)值分析方法(例如限差分法和有限元法等)會變得非常復(fù)雜.在發(fā)生較大形變的過程中,自由表面形狀可能會發(fā)生顯著變化,會發(fā)展成不連續(xù)的表面,網(wǎng)格也可能會崩潰.與之相比,沒有使用網(wǎng)格的粒子法在處理此類場景時具有一定優(yōu)勢,其可以使用一個簡單的算法對復(fù)雜的自由表面變化進行分析.因此,粒子法能夠有效表達當爆發(fā)到陸地上以及在建筑物周圍流動,水波發(fā)生復(fù)雜變化的3D 現(xiàn)象.另外,考慮復(fù)雜目標(包括建筑物)的分析模型時,粒子法的開發(fā)成本更低,因為在粒子法僅需要考慮粒子在目標區(qū)域中的放置問題.

在海嘯爆發(fā)的3D 分析中,無論使用的是網(wǎng)格方法還是粒子方法,都應(yīng)該事先確定空間分辨率(粒子法背景下為平均粒子距離)和分析準確度之間的關(guān)系,特別是在3D 分析中,由于極高的分析成本,應(yīng)該盡可能降低分辨率.

5 結(jié)論和展望

本文利用3D 粒子法(SPH 方法)開發(fā)出沿海城市地區(qū)的海嘯爆發(fā)分析模型.基于粒子法的空間分辨率和被淹沒區(qū)域之間的關(guān)系,得出了進行海嘯爆發(fā)分析所需的最小粒子距離.然后創(chuàng)建了具有真實感的3D 動畫快照.3D 仿真驗證了所提方法的有效性.通過比較在相同條件下進行的分析結(jié)果,確定了收斂解所需的最低空間分辨率.

未來,本文將著力開發(fā)能夠?qū)μ囟▍^(qū)域進行縮放的多步驟分析技術(shù),對不同區(qū)域采用不同的粒子距離.因此,不同方法的結(jié)合和縮放分析將是本文后續(xù)研究方向.