波浪作用下滲透率各向異性的海床液化分析

榮 富， 廖晨聰， 童大貴， 周香蓮

(上海交通大學(xué) 土木工程系， 上海 200240)

波浪荷載會(huì)在海床表面產(chǎn)生周期性波浪壓力，這種循環(huán)的波浪壓力作用會(huì)引起海床土體內(nèi)的孔隙水壓力和有效應(yīng)力隨時(shí)間和空間而發(fā)生變化，進(jìn)而導(dǎo)致海床土體液化失穩(wěn)，對(duì)海洋建筑物的地基穩(wěn)定性構(gòu)成極大威脅.因此，研究波浪荷載作用下滲透率各向異性的海床液化問題，對(duì)于近海和離岸結(jié)構(gòu)物的建設(shè)與設(shè)計(jì)具有重要意義.

關(guān)于波浪與海床之間相互作用的問題，自1941年Biot[1]提出Biot固結(jié)理論以來，國(guó)內(nèi)外學(xué)者已開展了大量研究[2-8].Madsen[9]根據(jù)Biot固結(jié)理論，應(yīng)用復(fù)變量表達(dá)式推導(dǎo)了線性波作用下海床土體為多孔彈性介質(zhì)時(shí)位移-孔隙水壓的解析解；Yamamoto等[10]利用半解析的方法求解線性推進(jìn)波作用下有限或無限深度多孔彈性介質(zhì)海床的動(dòng)力響應(yīng)；Mei等[11]采用邊界層近似法進(jìn)一步研究了此類問題并提出邊界層理論；王棟等[12]根據(jù)二維Biot固結(jié)理論，采用有限元數(shù)值模擬的方法研究了飽和多孔彈性介質(zhì)海床在線性波作用下的動(dòng)力響應(yīng).從20世紀(jì)90年代開始，國(guó)內(nèi)外學(xué)者對(duì)波浪導(dǎo)致海床液化失穩(wěn)的現(xiàn)象進(jìn)行了深入研究.Zen等[13]通過對(duì)海床動(dòng)力響應(yīng)的分析，探討了波浪作用下海床的致密化機(jī)制，并對(duì)振蕩孔隙水壓力引起的海床瞬時(shí)液化進(jìn)行分析；Sakai等[14]在分析波浪作用下海床瞬時(shí)液化時(shí)發(fā)現(xiàn)，利用邊界層孔隙水壓力的近似解能夠估計(jì)波谷區(qū)海床的瞬時(shí)液化深度.20世紀(jì)初期，Lin等[15]提出三維短峰波作用下各向異性的海床動(dòng)力響應(yīng)模型，并用于估計(jì)波浪導(dǎo)致的海床液化； Suzuki 等[16]采用室內(nèi)試驗(yàn)的方法研究了波浪作用下細(xì)砂海床的液化，并結(jié)合有限元模擬方法計(jì)算土體的最大液化深度.黃光爵等[17]在考慮振蕩孔隙水壓力和殘余孔隙水壓力的條件下推導(dǎo)出波浪作用下海床的最大液化深度的解析表達(dá)式.

但是，以上研究大多針對(duì)滲透率各向同性的海床的動(dòng)力響應(yīng)，對(duì)波浪荷載作用下滲透率各向異性的海床的動(dòng)力響應(yīng)及液化分析較少.Jeng等[18]考慮海床土體的各向異性，推導(dǎo)了波浪作用下有限深度海床動(dòng)力響應(yīng)的解析解；張金鳳等[19]采用數(shù)值計(jì)算方法分析了在非均質(zhì)各向異性海床中，土體的滲透系數(shù)對(duì)海床孔隙水壓力和有效應(yīng)力的影響；Wen等[20]采用有限單元法研究了不同滲透系數(shù)的2層土體在波浪荷載作用下的動(dòng)力響應(yīng).

為了研究滲透率各向異性的海床的瞬態(tài)液化特性，本文采用雷諾平均Navier-Stokes(RANS)方程作為波浪的控制方程，采用動(dòng)量源函數(shù)進(jìn)行造波，通過LSM(Level Set Method)法對(duì)波浪的自由表面進(jìn)行精確跟蹤.海床響應(yīng)以Biot多孔彈性方程為控制方程，并且考慮海床土體的滲透率各向異性.通過 COMSOL Multiphysics軟件來實(shí)現(xiàn)波浪模型和海床土體模型在接觸面上信息(位移、壓力)的交互耦合.同時(shí)，在驗(yàn)證數(shù)值模型合理性的基礎(chǔ)上，分析了波浪參數(shù)、海床的土體飽和度和滲透系數(shù)對(duì)海床瞬態(tài)最大液化深度的影響.

1 數(shù)值模型

波浪-海床相互作用的幾何模型如圖1所示.其中：Ls為模型的長(zhǎng)度；d為水深；h為海床土體的厚度；F為海床表面一點(diǎn).本文的數(shù)值模型包含2個(gè)模塊，即波浪模型和海床模型，以準(zhǔn)確模擬波浪與滲透率各向異性的海床的相互作用.采用RANS方程描述黏性流體流動(dòng)中波浪的產(chǎn)生和運(yùn)動(dòng)，采用Biot固結(jié)理論分析海床的動(dòng)力響應(yīng).波浪模型和海床模型之間采用單向耦合，即在模型計(jì)算的同一個(gè)時(shí)間步內(nèi)，先以波浪模型計(jì)算所得海床表面的波浪壓力作為海床模型的邊界條件，以此計(jì)算海床的位移及孔隙水壓力.

圖1 波浪-海床相互作用幾何模型Fig.1 Geometrical model of wave-seabed interactions

1.1 波浪模型

如圖1所示，在水槽的中間部分進(jìn)行數(shù)值造波，并采用LSM法追蹤波浪自由液面，而在水槽的兩端添加消波區(qū).以RANS方程作為二維不可壓縮流體運(yùn)動(dòng)的控制方程，即

μf為流體動(dòng)力黏度，k為湍流動(dòng)能，δij為Kronecker Delta記號(hào)，μt=ρfCuk2/ε，為湍動(dòng)流體黏度，Cu=0.09，ε為湍流動(dòng)能耗散率.

本文的流體力學(xué)性能模擬采用k-ε湍流模型，標(biāo)準(zhǔn)湍流模型的k方程和ε方程分別為

Gk-ρfε

(3)

(4)

式中：參數(shù)σk、σε、Cε1、和Cε2均為相關(guān)經(jīng)驗(yàn)系數(shù)，一般取σk=1.00，σε=1.30，Cε1=1.44，Cε2=1.92.

1.2 海床模型

在建立海床模型的過程中，假設(shè)波浪荷載引起的砂質(zhì)海床孔隙水滲流為層流運(yùn)動(dòng)，且滿足達(dá)西定律，則波浪荷載將引起多孔彈性海床土體內(nèi)產(chǎn)生超孔隙水壓力，使土骨架發(fā)生變形，并服從虎克定律.本文研究滲透率各向異性的海床土體，其連續(xù)性方程為

(5)

式中：Kx與Kz分別為海床土體在x和z軸方向的滲透率；ps為波浪作用下的孔隙水壓力；γw為孔隙水容重；ns為海床土體的孔隙率；us為海床土體的水平位移；ws為海床土體的豎向位移；β為孔隙水的壓縮性系數(shù)，其表達(dá)式為

Sr為土體的飽和度；Ew為土體孔隙的流體體積壓縮模量，其取值為Ew=2.24 GPa；pw0為絕對(duì)水壓力.

(6)

(7)

(8)

海床中超孔隙水壓力與有效應(yīng)力的關(guān)系為

忽略慣性力，將式(6)～(8)代入(9)和(10)，所得以超孔隙水壓力和位移表示的平衡微分方程為

(11)

(12)

1.3 邊界條件

1.3.1海床表面 在海床表面z=-d處，土體應(yīng)力相比于波浪壓力可忽略不計(jì)，故設(shè)土體垂直方向的正應(yīng)力和剪應(yīng)力等于0，并且海床表面的波浪壓力與孔隙水壓力相等，即有

式中：p0為作用于海床表面的波浪壓力.

1.3.2海床底面 在海床底面z=-(d+h)處，土體可視為不透水的剛性基巖，即海床土體的水平方向位移和豎直方向位移為0，故有

(15)

1.3.3海床側(cè)邊 假設(shè)在海床左右兩側(cè)邊界處，土體的水平方向位移為0且為不透水邊界，則有

(16)

2 模型驗(yàn)證

圖2 不同網(wǎng)格數(shù)目時(shí)點(diǎn)F的孔隙水壓力時(shí)程曲線Fig.2 Variations of pore pressure with time at point F for various mesh systems

由于滲透率各向異性的海床瞬態(tài)響應(yīng)的實(shí)測(cè)數(shù)據(jù)較少，所以本文選用波浪作用下滲透率各向同性的海床動(dòng)態(tài)響應(yīng)的試驗(yàn)數(shù)據(jù)和解析結(jié)果來驗(yàn)證模型.首先，驗(yàn)證波浪模型的網(wǎng)格收斂性.所用波浪的相關(guān)參數(shù)分別為波浪高度H=2 m，水深d=10 m，波浪周期T=12 s；海床的相關(guān)參數(shù)分別為海床厚度h=20 m，ns=0.4，滲透系數(shù)K=1.0×10-7m/s，Sr=0.97，νs=0.333，G=0.10 MPa.圖2所示為在不同網(wǎng)格密度下圖1中模型表面點(diǎn)F的孔隙水壓力ps的時(shí)程曲線.其中，N表示波浪模型的網(wǎng)格數(shù)目.由圖2可見，海床表面點(diǎn)F的孔隙水壓力不隨網(wǎng)格數(shù)目的變化而變化，經(jīng)過網(wǎng)格收斂性檢驗(yàn)，確定波浪模型網(wǎng)格數(shù)目取N=10 034.

為驗(yàn)證所建立的海床模型的合理性，本文將Liu等[21]的室內(nèi)模型試驗(yàn)值和Hsu等[22]的解析解與本文模型的數(shù)值計(jì)算結(jié)果進(jìn)行比較.根據(jù)Liu等的試驗(yàn)設(shè)置，計(jì)算所用參數(shù)分別為H=1.23 m，T=9 s，h=1.8 m，d=5 m，G=0.127 MPa，νs=0.3，K=1.8×10-4m/s，ns=0.425，Sr=0.975，計(jì)算所得沿不同海床深度(z/h)分布的歸一化孔隙水壓力(|ps|/p0)如圖3所示.由圖3可見，3種方法所得孔隙水壓力的變化趨勢(shì)基本一致，且吻合較好，說明本文模型可以有效模擬波浪作用下的海床響應(yīng).

圖3 沿不同海床深度分布的歸一化孔隙水壓力Fig.3 Vertical distributions of unified pore pressure versus seabed depth

3 海床瞬態(tài)液化分析

波浪荷載對(duì)海床作用時(shí)，波谷區(qū)水質(zhì)點(diǎn)的向上運(yùn)動(dòng)會(huì)對(duì)海床產(chǎn)生向上的吸力.本文采用Zen等[13]提出的海床瞬態(tài)液化判別標(biāo)準(zhǔn)進(jìn)行海床液化情況分析(見圖4)，即

(17)

pε=|p0|-|pw|

pw為波浪作用下海床中某點(diǎn)的孔隙水壓力.

圖4 海床瞬態(tài)液化判別示意圖Fig.4 Schematic graph of criteria for transient seabed liquefaction

3.1 波浪參數(shù)和海床土體飽和度的影響

波浪高度直接影響波浪作用在海床表面的波浪壓力，波浪周期通過影響波長(zhǎng)而影響海床表面的波浪壓力，海床土體飽和度通過影響土體的垂直孔壓梯度進(jìn)而影響液化深度.對(duì)于滲透率各向同性的海床，Chowdhury等[23]通過試驗(yàn)研究了波浪與海床相互作用下波浪周期和土體飽和度對(duì)海床液化的影響；Liu等[21]采用一維圓筒試驗(yàn)分析了波浪高度、波浪周期和土體飽和度對(duì)海床最大液化深度的影響.本文進(jìn)一步探討波浪高度、波浪周期以及海床土體飽和度對(duì)滲透率各向異性的海床液化情況.在研究過程中，分別選取H=2.0，3.0，4.0，5.0 m，T=10，12，14，16 s及Sr=0.95，0.96，0.97，0.98，滲透系數(shù)取Kx=1.0×10-7m/s，Kz=1.0×10-6m/s，其他參數(shù)與驗(yàn)證波浪模型網(wǎng)格收斂性時(shí)的相同.

海床瞬態(tài)最大液化深度dm隨波浪高度和波浪周期的變化規(guī)律如圖5(a)和(b)所示.由圖可見，波浪高度和波浪周期對(duì)dm都具有顯著影響.隨著波浪高度、波浪周期增大，作用在海床表面的波浪壓力不斷增大，從而導(dǎo)致海床較大深度范圍內(nèi)的土體有效應(yīng)力小于pε，滲透各向異性的海床的dm增加.結(jié)合圖4可知，有效應(yīng)力隨海床深度增加而呈現(xiàn)出線性增長(zhǎng)，而pε隨海床深度的增加有限，從而使得dm的增幅逐漸減緩.圖5(c)所示為不同海床土體飽和度時(shí)海床瞬態(tài)最大液化深度的變化規(guī)律.由圖可見，隨著海床土體飽和度增大，海床土體的垂直孔壓梯度減小，從而使得dm隨著土體飽和度的增加而減小，并且最大液化深度的減小速率基本保持不變.

圖5 海床瞬態(tài)最大液化深度隨波浪高度、波浪周期和海床土體飽和度變化的規(guī)律Fig.5 Curves of maximum liquefaction depth with wave height, wave period and saturation of soil

3.2 海床土體滲透性的影響

對(duì)于滲透率各向異性的海床， 在研究海床不同方向的滲透性對(duì)瞬態(tài)液化的影響時(shí)，分別選取Kx=1.0×10-7m/s與Kz=1.0×10-7m/s作為基準(zhǔn)進(jìn)行數(shù)值模擬，并且在分析一個(gè)方向的滲透系數(shù)對(duì)海床瞬態(tài)液化影響的過程中，固定另一個(gè)方向的滲透系數(shù).

3.2.1垂直方向滲透系數(shù)的影響 本文在Kx=1.0×10-7m/s，Kz處于 1.0×10-6～1.0×10-2m/s的范圍進(jìn)行分析.數(shù)值計(jì)算結(jié)果顯示，在Kz處于 1.0×10-6～1.0×10-2m/s時(shí)，海床土體的dm的變化劇烈，因此，在該范圍內(nèi)應(yīng)盡可能得多選取一些控制點(diǎn).

圖6(a)～(c)分別給出了在4組不同的H、T和Sr情況下，滲透率各向異性的海床瞬態(tài)最大液化深度隨著垂直方向滲透系數(shù)的變化規(guī)律.由圖可見，當(dāng)Kz<1.0×10-4m/s時(shí)，不同波浪條件下，dm隨著Kz的變化不明顯，較小的Kz使得海床中波浪所引起的淺層孔壓梯度較大，導(dǎo)致一定厚度的海床土體發(fā)生液化.當(dāng)Kz>1.0×10-4m/s時(shí)，對(duì)于H=2，3 m的海床，其波浪引起的海床土體的孔壓梯度迅速減小，液化難以向海床深部發(fā)展，dm隨著Kz的增大而急劇衰減；對(duì)于H=5 m的海床(如圖6(a)中點(diǎn)A、B和C所示)，隨著Kz不斷增大，所對(duì)應(yīng)的點(diǎn)A、B和C的垂直孔壓梯度先增后減(見圖7)，因此，dm先增后減；在不同的波浪周期和海床土體飽和度條件下，隨著Kz增大，dm急劇衰減，直至為0.

3.2.2水平方向滲透系數(shù)的影響 本文在Kz=1.0×10-7m/s，Kx處于 1.0×10-6～1.0×10-3m/s的條件下進(jìn)行分析.圖8示出了在不同的H、T和Sr情況下，滲透率各向異性的海床瞬態(tài)最大液化深度隨著Kx變化的規(guī)律.由圖8(a)可見，對(duì)于不同的H，隨著Kx增大，dm基本不變化.由圖8(b)可見：當(dāng)T=12，14，16 s時(shí)，dm隨Kx的變化不顯著；當(dāng)T=10 s，Kx>1.0×10-5m/s時(shí)，dm隨Kx增大而呈現(xiàn)出緩慢減小.由圖8(c)可見：當(dāng)Sr=0.95，0.96， 0.97 時(shí)，Kx對(duì)dm的影響不大；當(dāng)Sr=0.98，Kx>1.0×10-4m/s時(shí)，dm隨Kx的增大而略有減小.

圖6 海床瞬態(tài)最大液化深度隨著垂直方向滲透系數(shù)的變化規(guī)律Fig.6 Curves of maximum transient liquefaction depth with vertical permeability coefficient

圖7 不同垂直方向滲透系數(shù)下海床的孔隙水壓力示意圖Fig.7 Schematic graph of pore water pressure under different vertical permeabilities

圖8 海床瞬態(tài)最大液化深度隨著水平方向滲透系數(shù)的變化規(guī)律Fig.8 Curves of maximum transient liquefaction depth with horizontal permeability coefficient

4 結(jié)語

(1) 所建立的波浪-海床相互作用的數(shù)值模型與以往提出的解析解和相應(yīng)的室內(nèi)試驗(yàn)結(jié)果吻合較好，從而驗(yàn)證了其合理性，并可以有效模擬波浪荷載作用的海床動(dòng)力效應(yīng).

(2) 波浪高度、波浪周期及海床土體飽和度對(duì)滲透率各向異性的海床瞬態(tài)液化影響顯著，海床的瞬態(tài)最大液化深度隨著波浪高度和波浪周期的增大而增大，且其最大液化深度的增長(zhǎng)速度都逐漸減緩；海床瞬態(tài)最大液化深度隨海床土體飽和度的增大而減小，但其減小速率基本保持不變.

(3) 滲透率各向異性的海床垂直方向的滲透性對(duì)其瞬態(tài)液化影響顯著.當(dāng)海床土體垂直方向的滲透系數(shù)較小時(shí)，海床瞬態(tài)最大液化深度隨著垂直方向的滲透系數(shù)的變化不明顯；而當(dāng)滲透系數(shù)增大到一定數(shù)值時(shí)，海床瞬態(tài)最大液化深度隨著垂直方向的滲透系數(shù)的增大而急劇衰減.

(4) 海床水平方向的滲透系數(shù)對(duì)海床瞬態(tài)最大液化深度影響較小.