：本發(fā)明屬于堤防地質(zhì)災(zāi)害模擬，具體涉及一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法、裝置及介質(zhì)。

背景技術(shù)

0、
背景技術(shù)：

1、在堤防建設(shè)及地質(zhì)災(zāi)害預(yù)警分析中，滑坡涌浪次生災(zāi)害是十分重要的研究?jī)?nèi)容。臨水土質(zhì)岸坡失穩(wěn)破壞產(chǎn)生滑坡、滑坡體入水誘發(fā)涌浪的過(guò)程涉及到復(fù)雜的流固耦合作用，探明此作用機(jī)理對(duì)于有效反演滑坡涌浪演進(jìn)過(guò)程、建立地質(zhì)災(zāi)害預(yù)警系統(tǒng)有著重要意義。

2、現(xiàn)有數(shù)值模擬研究中，受制土體、水體的大變形特性及兩者間不同控制方程限制，開展相關(guān)多相耦合計(jì)算存在一定技術(shù)難題，現(xiàn)有的耦合算法對(duì)于局部水體壓強(qiáng)、土體輪廓發(fā)育過(guò)程的反演并不可靠，所得涌浪演進(jìn)波幅、土體輪廓發(fā)育軌跡均存在較大誤差。因此，現(xiàn)有的數(shù)值模擬方法無(wú)法準(zhǔn)確模擬土體滑坡涌浪的全過(guò)程。

技術(shù)實(shí)現(xiàn)思路

0、
技術(shù)實(shí)現(xiàn)要素：

1、本發(fā)明是為了解決上述不足而提出的，目的在于提供一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法、裝置及介質(zhì)，能夠解決現(xiàn)有數(shù)值模擬方法無(wú)法準(zhǔn)確模擬土體滑坡涌浪全過(guò)程的技術(shù)問(wèn)題。

2、為了實(shí)現(xiàn)以上目的，本發(fā)明采用了以下方案：

3、本發(fā)明設(shè)計(jì)一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，包括以下步驟：

4、步驟1、對(duì)所研究工程區(qū)域進(jìn)行精確建模，確定場(chǎng)地的滑坡、水體形狀特征；

5、步驟2、對(duì)計(jì)算區(qū)域進(jìn)行離散化處理，針對(duì)不規(guī)則的土體、水體邊界粒子進(jìn)行二次修正，確定計(jì)算所需相關(guān)材料特性參數(shù)；

6、步驟3、計(jì)算滑坡失穩(wěn)土體區(qū)域受力，在保證滑坡模擬精度的前提下，選擇合適的離散元接觸模型與離散元尺度，采用離散元方法進(jìn)行計(jì)算，得到滑坡體各區(qū)域加速度；

7、步驟4、水體及涌浪的模擬計(jì)算采用光滑粒子動(dòng)力學(xué)方法，河道水體區(qū)域離散為一系列流體粒子，通過(guò)積分格式的navier-stokes方程進(jìn)行計(jì)算，保證流體粒子尺寸與離散元接近；

8、步驟5、進(jìn)行土體、水體之間的法向耦合力計(jì)算，計(jì)算過(guò)程參照優(yōu)化全解析格式算法開展，基于以下公式進(jìn)行法向耦合加速度計(jì)算：

9、

10、式中：dui/dt代表耦合加速度，i代表離散粒子i，j代表離散粒子j，ui代表計(jì)算粒子i速度，t為時(shí)間，pi和pj分別為粒子i和j的流體壓強(qiáng)(耦合計(jì)算過(guò)程中土體顆粒壓強(qiáng)計(jì)算參考公式(3))，ρi為粒子i的密度，w為核函數(shù)，為核函數(shù)的位移偏微分形式，wij表示粒子i和j相對(duì)位置下的核函數(shù)取值，vj為粒子j的體積，μi和μj分別為粒子i和j的粘性系數(shù)，rij為粒子i和j之間距離，h為光滑長(zhǎng)度，g為重力加速度，取值-9.81m/s2，γ＝7，c0為人工聲速，ρ0χ為相χ的參考密度，fcoupling為耦合力，mi為粒子i的質(zhì)量，γ為與粒子i發(fā)生耦合作用的粒子集合；土體、水體間切向耦合力計(jì)算過(guò)程參考非解析經(jīng)驗(yàn)公式的離散格式；

11、步驟6、計(jì)算離散計(jì)算域內(nèi)各粒子的合加速度，進(jìn)行單時(shí)間步內(nèi)速度與位置更新、目標(biāo)時(shí)刻數(shù)據(jù)存儲(chǔ)，重復(fù)步驟1～步驟3，直至達(dá)到計(jì)算時(shí)長(zhǎng)要求。

12、步驟7、讀取導(dǎo)出數(shù)據(jù)，進(jìn)行可視化處理，進(jìn)行滿足計(jì)算要求的后續(xù)圖像分析。

13、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟2中，土體滑坡與水體的離散可借助現(xiàn)有前處理工具，旨在生成一系列距離相近、均勻分布的坐標(biāo)節(jié)點(diǎn)，以表征土體滑坡、水體實(shí)際輪廓，并為后續(xù)計(jì)算提供初始條件。

14、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟3中，具體包括以下子步驟：

15、步驟3.1、位于土體區(qū)域內(nèi)的系列離散節(jié)點(diǎn)設(shè)置為離散元(discrete?elementmethod，dem)顆粒，半徑為離散節(jié)點(diǎn)距離的0.5倍，通過(guò)計(jì)算顆粒受力累計(jì)細(xì)觀位移，獲取宏觀土體變形過(guò)程；

16、步驟3.2、顆粒離散元間接觸力計(jì)算公式如下所示：

17、fn＝knδnn+fndamp?????????????????????????????(4)

18、ft＝max|(ft)t-δt+ktδt+ftdamp，μ||fn||??????????????????????(5)式中：fn為法向接觸力，ft為切向接觸力，δn為粒子間法向重疊量，n為兩顆粒圓心連線方向，δt為粒子間切向重疊量，δt為單時(shí)間步內(nèi)接觸位置間的相對(duì)位移，fndamp與ftdamp分別表示法向與切向的阻尼力，(ft)t-δt為上個(gè)時(shí)間步結(jié)束時(shí)顆粒接觸面所受切向接觸力，μ為顆粒間動(dòng)摩擦系數(shù)，取兩者動(dòng)摩擦系數(shù)較小值，kn與kt分別為接觸界面法向剛度與切向剛度，其計(jì)算選擇hertz-mindlin接觸模型：

19、

20、式中：g*＝0.5(gi+gj)，r*＝2rirj/(ri+rj)，v*＝0.5(vi+vj)，i和j分別表示顆粒i和顆粒j；g為剪切模量，gi和gj分別表示顆粒i和j的剪切模量，r為顆粒半徑，ri和rj分別表示顆粒i和j的顆粒半徑，v為泊松比，vi和vj分別表示顆粒i和j的泊松比，soverlap表示顆粒i與顆粒j接觸重合距離；

21、步驟3.3、單顆粒所受接觸的合力計(jì)算參考以下公式：

22、

23、式中：fntotal，i為單顆粒所受接觸的法向合力，fttotal,i為單顆粒所受接觸的切向合力，fn，ij為顆粒i與顆粒j之間法向接觸力，ft，ij為顆粒i與顆粒j之間切向接觸力，n表示與顆粒i接觸作用的顆粒集合。

24、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟4中，具體包括以下子步驟：

25、步驟4.1、水體離散方法選擇光滑粒子流體動(dòng)力學(xué)(smoothed?particlehydrodynamics，sph)方法，位于水體區(qū)域內(nèi)的系列離散節(jié)點(diǎn)設(shè)置為sph粒子，半徑為離散節(jié)點(diǎn)間距，通過(guò)計(jì)算粒子受力累計(jì)細(xì)觀位移，獲取宏觀水體演進(jìn)過(guò)程；

26、步驟4.2、流體粒子受力計(jì)算參考以下公式：

27、

28、式中：dρi/dt為粒子i的密度梯度，dui/dt為粒子i的耦合加速度，ρi和ρj分別為粒子i和j的密度，t為時(shí)間，uj和ui分別為粒子j和i的速度矢量，為核函數(shù)的位移偏微分形式，wij表示粒子i和j相對(duì)位置下的核函數(shù)取值，mj為粒子j的質(zhì)量，pi和pj分別為粒子i和j的流體壓強(qiáng)，c0為人工聲速，為保證流體弱可壓假設(shè)，sph粒子密度變化梯度需小于1％，因此c0取值需不小于流場(chǎng)最大速度的10倍，ri和rj分別為粒子i和j位移矢量，ρ0為參考密度(20℃一個(gè)標(biāo)準(zhǔn)大氣壓下，水密度取值為1000kg/m3)，uij為粒子i和j的速度矢量差，δ為密度擴(kuò)散項(xiàng)，μ為粘性系數(shù)，γ＝7，g為重力加速度，取-9.81m/s2，w為核函數(shù)，wij＝w(ri-rj,h)，格式為：

29、

30、式中：q＝||ri-rj||/h，αdim是一個(gè)與計(jì)算維度相關(guān)的參數(shù)，二維工況下取7/(4πh2)，三維工況下取21/(16πh3)，光滑長(zhǎng)度h在二維工況下通常取1.5～2倍粒子初始間距，三維工況下取1.5倍粒子間距，能夠保證計(jì)算的精度；

31、步驟4.3、單流體粒子所受內(nèi)力計(jì)算參考以下公式：

32、

33、式中：ui為計(jì)算所得單顆粒速度矢量。

34、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟5中，具體包括以下子步驟：

35、步驟5.1，進(jìn)行耦合界面固、液相粒子檢索，對(duì)流體粒子影響域內(nèi)固相粒子進(jìn)行屬性修正，將其賦予相同密度流體sph粒子特性；

36、步驟5.2、計(jì)算土體-水體間耦合力，粒子間法向、切向計(jì)算力按照不同解析格式分別進(jìn)行計(jì)算，其中法向力計(jì)算參考全解析格式算法，具體計(jì)算見公式(1)～公式(3)；

37、步驟5.3、切向力計(jì)算公式為非解析經(jīng)驗(yàn)格式，具體計(jì)算公式如下：

38、

39、式中：i和k分別表示流體粒子i和土體顆粒k，為顆粒k所受切向力，為顆粒i所受切向力，vi和vj分別為粒子i和j的體積，mi為流體粒子i的質(zhì)量，pi為流體粒子壓強(qiáng)，w為核函數(shù)，εk為顆粒k的局部孔隙率，uk為顆粒k的平均速度，為插值后顆粒k周圍流體粒子的平均速度，ω和ω’均表示與目標(biāo)粒子作用的粒子集合，wik和wjk分別表示粒子i、k和j、k相對(duì)位置下的核函數(shù)取值，η為相間動(dòng)量傳遞系數(shù)，該參數(shù)與局部平均孔隙率和相對(duì)速度有關(guān)：

40、

41、式中：μf為流體粘性系數(shù)(20℃、一個(gè)標(biāo)準(zhǔn)大氣壓條件下，水取1.0×10-3pa·s)，rk為dem顆粒半徑，ρf為流體密度(標(biāo)準(zhǔn)工況下，水取1.0×103kg/m3)，cd為阻力系數(shù)，更多計(jì)算細(xì)節(jié)如下：

42、

43、式中：εk為顆粒k的局部孔隙率，εi為顆粒i的局部孔隙率，wik表示粒子i和k相對(duì)位置下的核函數(shù)取值，vk為粒子k的體積，mi為流體粒子i的質(zhì)量，rk為dem顆粒半徑，rek為粒子附近流場(chǎng)雷諾數(shù)。

44、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟6中，具體包括以下子步驟：

45、步驟6.1、進(jìn)行單時(shí)間步內(nèi)的預(yù)測(cè)步與修正步更新：

46、

47、式中：為粒子i在第(n+1/2)步后的速度矢量，為粒子i在第(n)步后的速度矢量，δt為時(shí)間步長(zhǎng)，為粒子i在第(n)步后的加速度矢量，為粒子i在第(n+1/2)步后的角速度矢量，為粒子i在第(n)步后的角速度矢量，為粒子i在第(n)步后的角加速度矢量，為粒子i在第(n+1/2)步后的轉(zhuǎn)角，為粒子i在第(n)步后的轉(zhuǎn)角，為粒子i在第(n+1/2)步后的密度，為粒子i在第(n)步后的密度，為粒子i在第(n)步后的密度變化梯度，為粒子i在第(n+1/2)步后的位移矢量，rin為粒子i在第(n)步后的位移矢量，為粒子i在第(n)步后的速度矢量，為第(n+1/2)步后的流體壓強(qiáng)，為第(n+1/2)步后的密度，f為比例系數(shù)，計(jì)算格式參考公式(12)；

48、其中，公式(23)為狀態(tài)方程的簡(jiǎn)化表達(dá)式，預(yù)測(cè)步后可以獲取第(n+1/2)步的速度u，角速度ω，角位移θ，密度ρ，線位移r，流體壓強(qiáng)p等參數(shù)，這些參數(shù)用于sph及dem的連續(xù)性方程、動(dòng)量方程的計(jì)算，以獲取第(n+1/2)步的加速度ain+1/2、密度梯度(dρ/dt)in+1/2等參數(shù)，為進(jìn)行下一步修正做準(zhǔn)備；

49、步驟6.2、修正步：

50、

51、式中：為粒子i在第(n+1/2)步后的速度矢量，為粒子i在第(n)步后的速度矢量，δt為時(shí)間步長(zhǎng)，為粒子i在第(n+1/2)步后的加速度矢量，為粒子i在第(n+1/2)步后的角速度矢量，為粒子i在第(n)步后的角速度矢量，為粒子i在第(n+1/2)步后的角加速度矢量，為粒子i在第(n+1/2)步后的轉(zhuǎn)角，為粒子i在第(n)步后的轉(zhuǎn)角，為粒子i在第(n+1/2)步后的密度，為粒子i在第(n)步后的密度，為粒子i在第(n+1/2)步后的密度變化梯度，為粒子i在第(n+1/2)步后的位移矢量，ri2為粒子i在第(n)步后的位移矢量；

52、步驟6.3、參數(shù)更新：

53、

54、pn+1＝f(ρn+1)????????????????????(27)式中：為粒子i在第(n+1)步后的速度矢量，為粒子i在第(n+1/2)步后的速度矢量，為粒子i在第(n)步后的速度矢量，為粒子i在第(n+1)步后的角速度矢量，為粒子i在第(n+1/2)步后的角速度矢量，為粒子i在第(n)步后的角速度矢量，為粒子i在第(n+1)步后的轉(zhuǎn)角，為粒子i在第(n+1/2)步后的轉(zhuǎn)角，為粒子i在第(n)步后的轉(zhuǎn)角，為粒子i在第(n+1)步后的密度，為粒子i在第(n+1/2)步后的密度，為粒子i在第(n)步后的密度，rin+1為粒子i在第(n+1)步后的位移矢量，為粒子i在第(n+1/2)步后的位移矢量，rin為粒子i在第(n)步后的位移矢量；pn+1為第(n+1)步后的流體壓強(qiáng)，ρn+1為第(n+1)步后的密度，f為比例系數(shù)，計(jì)算格式參考公式(12)；

55、dem-sph框架中，兩個(gè)模塊的時(shí)間步長(zhǎng)保持一致，并按以下公式(28)進(jìn)行確定：

56、

57、式中：δt為單步時(shí)長(zhǎng)，δtdem為dem顆粒的單步時(shí)長(zhǎng)，δtsph為sph顆粒的單步時(shí)長(zhǎng)，mi為顆粒離散元i質(zhì)量，kni為離散單元i最大接觸剛度，hi為sph粒子光滑長(zhǎng)度，v為流體運(yùn)動(dòng)粘性系數(shù)，ai為粒子加速度，c0為人工聲速，umax為流場(chǎng)最大流速。

58、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟6中，整體模擬方法中土體、水體框架的計(jì)算相互獨(dú)立，步驟5計(jì)算所得耦合力作為附加項(xiàng)在每個(gè)時(shí)間步添加到兩個(gè)框架中。整體而言，所提的耦合模型為半解析模型，浮力計(jì)算采用多相流模型直接求解，而拖曳力計(jì)算則基于非解析經(jīng)驗(yàn)公式。

59、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法，還可以具有這樣的特征：在步驟7中，導(dǎo)出數(shù)據(jù)的處理可借助現(xiàn)有的開源后處理軟件實(shí)現(xiàn)，以特定的局部、平面、立體、動(dòng)畫等形式開展后續(xù)分析研究，將模型更新進(jìn)行可視化。

60、本發(fā)明還設(shè)計(jì)了一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬裝置，能夠自動(dòng)實(shí)現(xiàn)上述<方法>，包括：

61、建模模塊，對(duì)所研究工程區(qū)域進(jìn)行精確建模，確定場(chǎng)地的滑坡、水體形狀特征；

62、前處理模塊，對(duì)計(jì)算區(qū)域進(jìn)行離散化處理，針對(duì)不規(guī)則的土體、水體邊界粒子進(jìn)行二次修正，并確定計(jì)算所需相關(guān)材料特性參數(shù)；

63、土體滑坡模擬模塊，計(jì)算滑坡失穩(wěn)土體區(qū)域受力，在保證滑坡模擬精度的前提下，選擇合適的離散元接觸模型與離散元尺度，采用離散元方法進(jìn)行計(jì)算，得到滑坡體各區(qū)域加速度；

64、涌浪模擬模塊，水體及涌浪的模擬計(jì)算采用光滑粒子動(dòng)力學(xué)方法，河道水體區(qū)域離散為一系列流體粒子，通過(guò)積分格式的navier-stokes方程進(jìn)行計(jì)算，保證流體粒子尺寸與離散元接近；

65、滑坡涌浪耦合計(jì)算模塊，進(jìn)行土體、水體之間的法向耦合力計(jì)算，計(jì)算過(guò)程參照優(yōu)化全解析格式算法開展，基于以下公式進(jìn)行法向耦合加速度計(jì)算：

66、

67、式中：dui/dt代表粒子i的耦合加速度，i代表離散粒子i，j代表離散粒子j，ui代表計(jì)算粒子i速度，t為時(shí)間，pi和pj分別為粒子i和j的流體壓強(qiáng)，ρi為粒子i的密度，w為核函數(shù)，為核函數(shù)的位移偏微分形式，wij表示粒子i和j相對(duì)位置下的核函數(shù)取值，vj為粒子j的體積，μi和μj分別為粒子i和j的粘性系數(shù)，rij為粒子i和j之間距離，h為光滑長(zhǎng)度，g為重力加速度，取值-9.81m/s2，γ＝7，c0為人工聲速，ρ0χ為相χ的參考密度，fcoupling為耦合力，mi為粒子i的質(zhì)量，γ為與粒子i發(fā)生耦合作用的粒子集合；

68、時(shí)間步更新模塊，計(jì)算離散計(jì)算域內(nèi)各粒子的合加速度，進(jìn)行單時(shí)間步內(nèi)速度與位置更新、目標(biāo)時(shí)刻數(shù)據(jù)存儲(chǔ)，重復(fù)步驟1～步驟3，直至達(dá)到計(jì)算時(shí)長(zhǎng)要求；

69、后處理模塊，讀取導(dǎo)出數(shù)據(jù)，進(jìn)行可視化處理，進(jìn)行滿足計(jì)算要求的后續(xù)圖像分析；

70、控制模塊，與所述建模模塊、所述前處理模塊、所述土體滑坡模擬模塊、所述涌浪模擬模塊、所述滑坡涌浪耦合計(jì)算模塊、所述時(shí)間步更新模塊、所述后處理模塊均通信相連，控制它們的運(yùn)行。

71、優(yōu)選地，本發(fā)明提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬裝置，還可以具有這樣的特征：對(duì)實(shí)際工程進(jìn)行概化建模，并劃分重點(diǎn)研究區(qū)域，以stl格式或者dxf格式生成數(shù)值模型，導(dǎo)入前處理軟件中，進(jìn)行網(wǎng)格劃分，生成節(jié)點(diǎn)坐標(biāo)導(dǎo)入后續(xù)計(jì)算模塊。

72、本發(fā)明的有益效果：

73、本發(fā)明所提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬方法、裝置及介質(zhì)，解決了耦合界面局部水體壓強(qiáng)、土體輪廓發(fā)育過(guò)程的反演不精確的問(wèn)題，能夠有效表征滑坡體、水體間雙相耦合作用，并準(zhǔn)確模擬土體滑坡涌浪誘發(fā)傳播過(guò)程。能合理模擬耦合界面的壓強(qiáng)分布，消除非物理高壓間隙，對(duì)土體大變形運(yùn)動(dòng)過(guò)程、涌浪傳播過(guò)程均能進(jìn)行準(zhǔn)確預(yù)測(cè)，對(duì)于局部的水體沖擊翻卷流態(tài)也能進(jìn)行精確反演，并提供可視化分析材料。本發(fā)明所提供的一種計(jì)算土體滑坡涌浪的無(wú)網(wǎng)格模擬裝置能夠精確模擬快速、慢速土體滑坡涌浪等多種工況，應(yīng)用于典型斷面、區(qū)域內(nèi)的相關(guān)研究分析，為厘清土體滑坡涌浪演化規(guī)律及涌浪災(zāi)害防治提供指導(dǎo)。因此，本發(fā)明能夠準(zhǔn)確模擬土體滑坡涌浪全過(guò)程。