基于B樣條物質點法的潰壩流模擬研究

徐云卿, 周曉敏,2, 趙世一, 徐聖飛, 孫 政,2

(1. 江西理工大學 土木與測繪工程學院, 江西 贛州 341000;2. 江西理工大學 江西省環境巖土與工程災害控制重點實驗室, 江西 贛州 341000)

0 引 言

水庫大壩是土木工程中常見的工程結構,在防洪、能源生產、蓄水等方面都發揮著重要作用．但大壩在具有社會經濟效益的同時,也存在著潰壩的潛在風險,一旦發生潰壩,將會對人民生命財產安全造成極大的威脅．因此,開展水壩的潰決過程及其演進規律研究,對于潰壩洪水災害的預測和防治具有重要意義．

潰壩流是一種典型的自由表面流動問題,研究內容主要包括波前到達時間、波前流速、波前位置和沿程各處的水位等[1]．目前,研究方法主要有理論分析[2-4]、物理試驗[5-6]、數值模擬[7-8]．隨著數學理論和計算機技術的高速發展,數值模擬廣泛應用于潰壩水流的研究．數值模擬研究相比較于物理試驗具有更高的靈活性,不局限于場地設置,計算周期相對較短,能夠大幅度降低試驗成本等優點．Han等[8]采用Preissmann法的一維模型耦合采用有限差分法的潰壩二維模型求解潰壩擴散波動情況;胡四一和譚維炎[9]以TVD格式預測了高壩瞬潰的潰壩波演進過程;王大國等[10]采用CBOS有限元法建立了水波模型,精確模擬和分析了下游河床無水時潰壩模型的自由水面運動特征;張建偉等[11]采用光滑粒子流體動力學方法(SPH)建立了包含下游水位的潰壩模型,對于潰壩流沖擊過程中波前到達時間、波前流速、波前位置能夠做到精確模擬．

但潰壩流涉及到水波演進大變形問題,有網格類的有限差分法(FDM)[8]、有限體積法(FVM)[7]、有限元法(FEM)[10]等方法在模擬時會出現網格畸變并引起誤差．盡管無網格化的SPH[11]避免了網格畸變,但由于追蹤的是網格邊界上質量、動量和能量的通量流動,需求解非線性對流項,增加了求解難度,且不易于追蹤各質點的運動時間歷程,使得計算效率低．因此,需要開發一種可以減少誤差和提高計算效率的數值模擬方法．

物質點法(MPM)[12]作為一種相對新型的粒子型算法,其具有Lagrange 粒子和 Euler背景網格雙重描述,既包含了Lagrange粒子型無網格算法的優勢,又可避免網格畸變和求解非線性對流項,可方便地追蹤自由表面和粒子信息的時間歷程; 同時利用Euler型背景網格求解動量方程和梯度,具有計算效率和處理基本邊界條件[13-14]的優勢．物質點法已廣泛應用各個工程領域問題[15-19]中．然而,傳統的物質點法采用線性插值形函數來實現物質點和背景網格節點之間的信息映射．由于網格邊界處的形函數梯度不連續,當物質點穿越背景網格邊界時,將產生網格穿越誤差,引起數值振蕩．Gan等[20]采用高階B樣條基函數[21]替換線性插值基函數,節點自由度空間代替背景網格節點,減少了網格穿越誤差,提高了物質點法的計算精度和收斂性．目前B樣條物質點法(BSMPM)已經用于模擬自由表面流體問題,如具有自由表面[22-23]的復雜流動問題,流體-結構相互作用問題[24]等．

1 計 算 方 法

1.1 弱可壓縮B樣條物質點法

1.1.1 B樣條基函數

B樣條基函數的定義有很多種,盡管所得的表達式不同,但其本質是一致的．本文選用計算穩定、易于理解的Cox-de Boor遞歸公式[25]來計算B樣條基函數．在Cox-de Boor遞歸公式中,零階基函數Si,0(q=0)如下所示:

(1)

對于q≥1,基函數由Cox-de Boor遞歸公式定義為

(2)

式中,ξi是第i個節點,i=1,2,…,n+q+1,n是基函數的總數目,q是基函數的階數．由式(1)可以看出,如果基函數的階數為零(q=0),則基函數呈階梯分布,即在第i節點區間[ξi,ξi+1)上基函數Si,0(ξ)等于1,而其他節點區間等于0．

采用Si,0(ξ)和Si+1,0(ξ),通過式(2)可以得到Si,1(ξ);當所有Si,1(ξ)計算完成,采取同樣的方法可以得到Si,2(ξ),持續這個計算過程直到所有需要的Si,q(ξ)計算完畢,其遞歸過程如圖1所示．

圖1 Cox-de Boor遞歸示意圖Fig. 1 Schematic diagram of the Cox-de Boor recursion

(a) 參數網格空間 (b) 節點自由度空間(a) The parametric grid space(b) The node degree of freedom space圖2 網格空間示意圖Fig. 2 Illustration of grid spaces

假設0/0=0,Si,q(ξ)的導數也可以用遞歸公式定義,即

(3)

本文研究的是二維潰壩流模擬實驗,基函數是雙變量基函數,假定節點向量N={ξ1,ξ2,…,ξn+q,ξn+q+1}和H={η1,η2,…,ηm+p,ηm+p+1},ξ和η是參數的方向,n和m分別為沿圖2(a)中ξ,η方向的基函數數目．利用張量積結構雙變量基函數定義為

S(i,j),(q,p)(ξ,η)=Si,q(ξ)·Sj,p(η),

(4)

式中,i=1,2,…,n;j=1,2,…,m．三變量B樣條基函數也可以類似于雙變量基函數推導得到．圖2(a)參數網格中的空心圓和實心圓分別表示節點向量和物質點．

圖2(b)給出的是節點自由度空間,圖中(i,j)表示節點自由度空間中對應的(i,j)自由度;物質點(ξ,η)在節點自由度(i,j)上對應的B樣條基函數記作S(i,j),(q,p)(ξ,η),其中(q,p)分別表示在節點矢量N和H方向上B樣條基函數的階數．

1.1.2 B樣條物質點法

B樣條物質點法是基于物質點法改進而來的,B樣條插值基函數階數為一階時與物質點法中的傳統線性插值形函數是一樣的,即B樣條插值基函數階數為一階時的B樣條物質點法就是物質點法[26],兩者求解步驟基本一致．兩者的不同之處在于:在物質點法中,物質點信息采用線性插值形函數映射到相應背景網格節點上,在背景網格節點上施加本質邊界條件．而在B樣條物質點法中,物質點信息采用B樣條基函數映射到節點自由度空間的各節點自由度上,在節點自由度空間施加本質邊界條件．其具體算法實現過程如下:

1) 劃分參數網格,將研究對象離散成一組物質點．

2) 初始化物質點的位置、質量、動量等物質信息．

3) 采用B樣條基函數,在第k個計算時間步,將物質點上的質量和動量信息映射到節點自由度空間的節點自由度上:

(5)

(6)

(7)

(8)

式中,?S(i,j),(q,p)表示B樣條插值基函數的梯度;σ為Cauchy應力張量;b和t分別為體力和表面力矢量;V為體積;h表示面力邊界厚度．

5) 計算節點自由度上的加速度和速度:

(9)

(10)

式中,a表示加速度矢量．

6) 在節點自由度空間,施加本質邊界條件．

7) 將各節點自由度的加速度和速度信息,通過B樣條基函數映射回各物質點,得到第k+1個時間步上物質點的速度和位置:

(11)

(12)

8) 將更新后的各物質點速度重新映射回各節點自由度上,并計算物質點第k+1個時間步上的應變增量Δε:

(13)

(14)

9) 更新物質點的密度、應力、應變、位置等信息,開始下一計算步．

1.2 弱可壓縮流體本構模型

潰壩流問題中流體的本構模型可表示為[27]

(15)

(16)

2 計算模型和結果分析

2.1 計算模型

模型布置如圖3所示,圖3(a)中h0是上游水庫水位高度,l0是大壩水庫區域的長度,l和h分別是水箱的長度和高度,v是閥門上升速度．圖3(b)中L(t)是t時刻流體前緣位置,h(t)是t時刻上游水位高度．邊界條件取液柱的左右兩側和上下底面都為可滑移邊界條件, 初始下游河床無水．在模擬計算中取重力加速度g=9.81 m/s2,水的密度ρ=1 000 kg/m3．閥門上升速度v=2.0 m/s,人工聲速c和時間步長Δt分別為100 m/s和1×10-5s．

(a) 初始時刻(b) t時的水流位置(a) The initial setup (b) The flow position at time t圖3 潰壩流問題示意圖Fig. 3 Illustration of the dam break flow problem

本文基于三維程序模擬二維平面應變問題, 即在z方向僅設置一個背景網格, 單元網格大小為0.02 m,每個背景單元網格內布置2×2×2個物質點粒子．為了方便進行對比,采用無量綱化的L*-T*,H*-T*曲線表示．L*代表無量綱化后流體波前到達位置,H*代表無量綱化后給定位置水位高程,T*代表流體流動時間:

(17)

(18)

(19)

2.2 模擬結果分析

圖4 潰壩流體波前位置隨時間的變化Fig. 4 Variation of the dam-break fluid wavefront position with time

為了說明B樣條插值基函數階數對結果的影響,圖5(a)、(b)、(c)和(d)分別給出了給定位置h1,h2,h3和h4水位高程隨時間的變化．從圖中可以明顯看出: 1) 模擬得到的結果與實驗結果曲線單調性相近,模擬結果與實驗結果基本吻合．2) 隨著B樣條插值基函數階數的增加模擬結果曲線變得愈發平穩和光滑,在B樣條插值基函數階數為3階時,模擬實驗所得到的給定位置h1,h2,h3和h4水位高程隨時間變化關系曲線最為平穩和光滑．3) 與實驗結果相比,不同B樣條插值基函數階數下,模擬給定位置h1,h2,h3和h4水位高程隨時間變化關系曲線過程存在較為明顯的鋸齒變化．B樣條插值基函數階數的增加產生影響和模擬水位高程隨時間變化關系曲線結果存在明顯鋸齒變化的主要原因是: 1) 在B樣條物質點法中,隨著B樣條插值基函數階數的增加,節點矢量內的節點和基函數個數將會相應地增加．2) 在樣條物質點法中,隨著B樣條插值基函數階數的增加,基函數更為光滑,同時基函數在邊界處具有更大的梯度,提高了B樣條物質點法的求解精度和收斂性．3) 在數值模擬中,對水位高度的變化,是通過捕捉固定水平位置處一個背景網格單元范圍內物質點y坐標的最大值得到的,然而試驗數據是基于100組試驗的平均值,因此圖中試驗結果較光滑,而數值模擬結果呈鋸齒變化,但兩者吻合較好．由圖4和圖5可知,本文的模擬方法可以準確地模擬潰壩流體的傳播過程以及水庫區和下游區的水位高程變化規律．

圖5 初始水位高度h0=0.3 m,(a)、(b)、(c)和(d)分別為給定位置h1,h2,h3和h4水位高程隨時間的變化Fig. 5 Initial water level height h0=0.3 m, (a),(b),(c) and (d) denoting the changes of the water level elevation at given positions h1,h2,h3 and h4 with time, respectively

為了進一步說明潰壩流的演進過程和B樣條插值基函數階數對結果的影響,圖6給出了不同時刻速度剖面下B樣條物質點法1階、2階、3階的模擬結果和實驗結果的對比．從圖中可以看出: 1) 模擬結果的速度剖面與實驗結果有良好的一致性,且隨著B樣條插值基函數階數的增加,速度剖面愈發一致; 2) 隨著時間的推移,潰壩水流的波前位置的移動速度加快,右側水庫內水位逐漸降低．

圖6 初始水位高度h0=0.3 m,(a)、(b)、(c)和(d)分別為BSMPM 1階、2階、3階的模擬結果和實驗結果在t=0.32 s(T*=1.83),t=0.41 s(T*=2.34),t=0.46 s(T*=2.63)下的速度云圖Fig. 6 Initial water level height h0=0.3 m, (a),(b),(c) and (d) denoting the velocity profiles of the BSMPM 1st-order, 2nd-order and 3rd-order simulation results and experimental results at t=0.32 s(T*=1.83),t=0.41 s(T*=2.34),t=0.46 s(T*=2.63), respectively

2.3 計算效率

通過改變網格尺寸,對2.1小節所述的潰壩流問題,1階、2階和3階基函數下B樣條物質點法求解耗時進行分析．表1給出了1階、2階和3階基函數下B樣條物質點法,在網格尺寸為0.01 m,0.02 m和0.04 m下的單步CPU計算耗時,圖7繪制了其變化關系曲線．本文程序運算所用計算機為64位CentOS Linux 7系統、Inter Xeon Gold 6226R @ 2.90 GHz×64 CPU、128 G內存;程序基于InterFortran90編輯器,串行計算,CPU耗時通過CPU—TIME命令得到．

表1 不同網格尺寸下,1階、2階和3階基函數下B樣條物質點法的求解耗時

由表1和圖7可以看出: 1) 隨著基函數階數的增加,B樣條物質點法求解計算耗時呈約1.5倍增長; 2) 不同階次B樣條物質點法的計算耗時隨背景網格尺寸的增長率基本一致,約呈線性增長．

圖7 單步CPU計算耗時隨網格尺寸和基函數階數的變化Fig. 7 Variation of the single-step CPU computation time with the grid size and the order of basis functions

3 結 論

本文基于B樣條物質點法,通過引入人工狀態方程,開展了弱可壓縮B樣條物質點法模擬潰壩流問題的研究,分析了B樣條插值基函數階數對模擬結果的影響．從模擬結果可以得出:

1) B樣條物質點法所得潰壩流問題模擬結果與物理實驗結果基本吻合．

2) B樣條物質點法可以很好地模擬潰壩流體的流動特性．潰壩水流的波前速度,隨著時間的推進越來越快．對于給定位置的h1高程,隨著時間的推進逐漸下降,而h2,h3和h4的高程,隨著時間的推進逐漸上升,驗證了B樣條物質點法模擬潰壩流問題的可行性．

3) 在改變B樣條插值基函數階數的條件下,通過對比1階、2階和3階基函數下的潰壩流體波前位置、波前速度及給定位置的高程變化的模擬結果和物理實驗結果,得出隨著基函數階數的增加,模擬結果與實驗結果擬合度提高．

4) 在改變背景網格尺寸條件下,通過對比1階、2階和3階基函數B樣條物質點法計算耗時,得出相較物質點法,更高階B樣條物質點法的計算耗時約呈1.5倍增長;隨著背景網格數目的增加,各階次B樣條物質點法的計算耗時增長率與傳統物質點法基本一致,約呈線性增長．