基于GIS光滑粒子流體動力學改進方法仿真①

王建華,冉煜琨

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

1 引言

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

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

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

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

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

2 提出的分析方法

2.1 SPH 方法的基本方程

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

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

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

2.2 流體分析的控制方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中,n為墻壁的內向法線向量.

3 分析模型的開發

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

圖1 粒子模型的開發流程

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

圖2 降低STL 數據的建筑物底部

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

圖3 地面粒子

4 分析與評價

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

4.1 分析條件

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

圖4 分析模型的示意圖

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

表1 分析條件

4.2 結果評價

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

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

圖5 海嘯爆發分析示意圖

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

圖6 海嘯流入150 s 后被淹沒的區域

4.3 照片級的可視化

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

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

圖7 可視化方法的程序

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

4.4 與網格法的比較分析

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

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

5 結論和展望

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

未來,本文將著力開發能夠對特定區域進行縮放的多步驟分析技術,對不同區域采用不同的粒子距離.因此,不同方法的結合和縮放分析將是本文后續研究方向.