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

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

0 引 言

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

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

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

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

1 計 算 方 法

1.1 弱可壓縮B樣條物質(zhì)點法

1.1.1 B樣條基函數(shù)

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

(1)

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

(2)

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

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

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

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

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

(3)

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

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

(4)

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

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

1.1.2 B樣條物質(zhì)點法

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

1) 劃分參數(shù)網(wǎng)格,將研究對象離散成一組物質(zhì)點．

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

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

(5)

(6)

(7)

(8)

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

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

(9)

(10)

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

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

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

(11)

(12)

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

(13)

(14)

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

1.2 弱可壓縮流體本構(gòu)模型

潰壩流問題中流體的本構(gòu)模型可表示為[27]

(15)

(16)

2 計算模型和結(jié)果分析

2.1 計算模型

模型布置如圖3所示,圖3(a)中h0是上游水庫水位高度,l0是大壩水庫區(qū)域的長度,l和h分別是水箱的長度和高度,v是閥門上升速度．圖3(b)中L(t)是t時刻流體前緣位置,h(t)是t時刻上游水位高度．邊界條件取液柱的左右兩側(cè)和上下底面都為可滑移邊界條件, 初始下游河床無水．在模擬計算中取重力加速度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方向僅設(shè)置一個背景網(wǎng)格, 單元網(wǎng)格大小為0.02 m,每個背景單元網(wǎng)格內(nèi)布置2×2×2個物質(zhì)點粒子．為了方便進行對比,采用無量綱化的L*-T*,H*-T*曲線表示．L*代表無量綱化后流體波前到達位置,H*代表無量綱化后給定位置水位高程,T*代表流體流動時間:

(17)

(18)

(19)

2.2 模擬結(jié)果分析

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

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

圖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樣條插值基函數(shù)階數(shù)對結(jié)果的影響,圖6給出了不同時刻速度剖面下B樣條物質(zhì)點法1階、2階、3階的模擬結(jié)果和實驗結(jié)果的對比．從圖中可以看出: 1) 模擬結(jié)果的速度剖面與實驗結(jié)果有良好的一致性,且隨著B樣條插值基函數(shù)階數(shù)的增加,速度剖面愈發(fā)一致; 2) 隨著時間的推移,潰壩水流的波前位置的移動速度加快,右側(cè)水庫內(nèi)水位逐漸降低．

圖6 初始水位高度h0=0.3 m,(a)、(b)、(c)和(d)分別為BSMPM 1階、2階、3階的模擬結(jié)果和實驗結(jié)果在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 計算效率

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

表1 不同網(wǎng)格尺寸下,1階、2階和3階基函數(shù)下B樣條物質(zhì)點法的求解耗時

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

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

3 結(jié) 論

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

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

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

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

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