基于B樣條物質(zhì)點(diǎn)法的潰壩流模擬研究

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

(1. 江西理工大學(xué) 土木與測(cè)繪工程學(xué)院, 江西 贛州 341000;2. 江西理工大學(xué) 江西省環(huán)境巖土與工程災(zāi)害控制重點(diǎn)實(shí)驗(yàn)室, 江西 贛州 341000)

0 引 言

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

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

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

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

1 計(jì) 算 方 法

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

1.1.1 B樣條基函數(shù)

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

(1)

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

(2)

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

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

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

(a) 參數(shù)網(wǎng)格空間 (b) 節(jié)點(diǎn)自由度空間(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(ξ)的導(dǎo)數(shù)也可以用遞歸公式定義,即

(3)

本文研究的是二維潰壩流模擬實(shí)驗(yàn),基函數(shù)是雙變量基函數(shù),假定節(jié)點(diǎn)向量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ù)推導(dǎo)得到．圖2(a)參數(shù)網(wǎng)格中的空心圓和實(shí)心圓分別表示節(jié)點(diǎn)向量和物質(zhì)點(diǎn)．

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

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

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

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

2) 初始化物質(zhì)點(diǎn)的位置、質(zhì)量、動(dòng)量等物質(zhì)信息．

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

(5)

(6)

(7)

(8)

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

5) 計(jì)算節(jié)點(diǎn)自由度上的加速度和速度:

(9)

(10)

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

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

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

(11)

(12)

8) 將更新后的各物質(zhì)點(diǎn)速度重新映射回各節(jié)點(diǎn)自由度上,并計(jì)算物質(zhì)點(diǎn)第k+1個(gè)時(shí)間步上的應(yīng)變?cè)隽喀う?

(13)

(14)

9) 更新物質(zhì)點(diǎn)的密度、應(yīng)力、應(yīng)變、位置等信息,開(kāi)始下一計(jì)算步．

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

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

(15)

(16)

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

2.1 計(jì)算模型

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

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

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

(17)

(18)

(19)

2.2 模擬結(jié)果分析

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

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

圖5 初始水位高度h0=0.3 m,(a)、(b)、(c)和(d)分別為給定位置h1,h2,h3和h4水位高程隨時(shí)間的變化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

為了進(jìn)一步說(shuō)明潰壩流的演進(jìn)過(guò)程和B樣條插值基函數(shù)階數(shù)對(duì)結(jié)果的影響,圖6給出了不同時(shí)刻速度剖面下B樣條物質(zhì)點(diǎn)法1階、2階、3階的模擬結(jié)果和實(shí)驗(yàn)結(jié)果的對(duì)比．從圖中可以看出: 1) 模擬結(jié)果的速度剖面與實(shí)驗(yàn)結(jié)果有良好的一致性,且隨著B(niǎo)樣條插值基函數(shù)階數(shù)的增加,速度剖面愈發(fā)一致; 2) 隨著時(shí)間的推移,潰壩水流的波前位置的移動(dòng)速度加快,右側(cè)水庫(kù)內(nèi)水位逐漸降低．

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

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

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

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

圖7 單步CPU計(jì)算耗時(shí)隨網(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ì)點(diǎn)法,通過(guò)引入人工狀態(tài)方程,開(kāi)展了弱可壓縮B樣條物質(zhì)點(diǎn)法模擬潰壩流問(wèn)題的研究,分析了B樣條插值基函數(shù)階數(shù)對(duì)模擬結(jié)果的影響．從模擬結(jié)果可以得出:

1) B樣條物質(zhì)點(diǎn)法所得潰壩流問(wèn)題模擬結(jié)果與物理實(shí)驗(yàn)結(jié)果基本吻合．

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

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

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