近?；聸_擊樁體的數(shù)值模擬*

單治鋼，李 陽，王 棟**，董友扣

(1.浙江華東建設(shè)工程有限公司，浙江 杭州 310006；2.中國海洋大學(xué)環(huán)境科學(xué)與工程學(xué)院，山東 青島 266100；3.中國地質(zhì)大學(xué)(武漢)海洋學(xué)院，湖北 武漢 430074)

近海工程建設(shè)中廣泛采用樁基礎(chǔ)支撐橋梁、油氣平臺和風(fēng)機(jī)，但樁基礎(chǔ)可能受到海底滑坡威脅：水下邊坡失穩(wěn)后形成的流滑體沖擊滑動路線上的樁身，造成樁基乃至上部結(jié)構(gòu)的破壞。近年發(fā)展起來的剪切帶漸進(jìn)擴(kuò)展理論能夠較好地描述海底土質(zhì)邊坡的失穩(wěn)機(jī)制[1]，但追蹤失穩(wěn)后的流滑形態(tài)及其對下構(gòu)筑物的沖擊主要還是依靠模型試驗與數(shù)值模擬。流滑體由土水混合組成，力學(xué)性質(zhì)接近非牛頓流體[2-4]。當(dāng)流滑體以較高速度沖擊構(gòu)筑物，法向(垂直于構(gòu)筑物軸線方向)沖擊壓力可借助流體力學(xué)中的繞流拖曳力形式表達(dá)[5]：

(1)

式中：沖擊壓力p指單位投影面積上的沖擊力；CD為拖曳力系數(shù)，是雷諾數(shù)Re的函數(shù)[6-7]；ρ為流滑體密度；v為沖擊速度。Zakeri等[8-9]和Liu等[10-11]分別結(jié)合非牛頓流體的Herschel-Bulkley模型和冪律模型，按照式(1)給出了高速流滑體對懸浮或放置在海床表面的管線的沖擊壓力。當(dāng)流滑體以較低速度沖擊構(gòu)筑物時，式(1)中CD隨Re劇烈變化。如果沖擊速度接近于零，更適合采用土力學(xué)中的地基基礎(chǔ)承載力公式估算沖擊壓力[12]：

p=NCsu。

(2)

式中:NC為承載力系數(shù)；su為流滑體擬靜力不排水抗剪強(qiáng)度，即流體力學(xué)中的剪應(yīng)力。

為統(tǒng)一描述低速和高速流滑的沖擊作用，Randolph和White[13]建議將流滑體對有埋深的海底管線沖擊壓力分解為拖曳力項和承載力項兩部分：

(3)

如果流滑體自重形成的上覆壓力也影響沖擊壓力，式(3)的承載力項可進(jìn)一步細(xì)分為兩部分，分別與流滑體的不排水抗剪強(qiáng)度和上覆壓力相關(guān)[14]。但在以往的研究中，公式(3)僅被用于流滑體沖擊海底管線問題，不清楚是否適合流滑體沖擊樁身問題?，F(xiàn)有的沖擊樁基研究仍多采用式(1)總結(jié)沖擊力[15-16]，低雷諾數(shù)時預(yù)測沖擊力易導(dǎo)致較大誤差。

本文采用計算流體力學(xué)軟件Fluent，模擬不同速度的海底流滑體對樁身的沖擊過程，關(guān)注低雷諾數(shù)下的沖擊形態(tài)。在驗證數(shù)值模型后，進(jìn)行大量變動參數(shù)分析，考察影響沖擊壓力的主要因素，最終提出適應(yīng)寬廣雷諾數(shù)范圍的沖擊力表達(dá)式。

1 數(shù)值模型構(gòu)建

1.1 問題定義

薄層滑坡厚度較小[17]，坡腳侵蝕和快速沉積誘發(fā)的近?；乱?guī)模不大，流滑體高度小，不考慮其自重的影響，將流滑體對樁的沖擊簡化為圖1所示的平面應(yīng)變問題。假定單樁為位置固定的完全粗糙剛體，樁徑為D，在中低速沖擊情況下不考慮流滑體周圍水阻力的影響。采用基于有限體積法的商業(yè)軟件Fluent，實(shí)施“流滑體-氣”兩相流模擬。

圖1的計算域左邊界為流體恒速進(jìn)口，速度等于流滑體初始沖擊速度v，用于表征流滑體到達(dá)單樁前的速度；右邊界為壓力出口，壓力等于大氣壓；上下為自由滑移壁面。各邊界與樁的距離必須足夠遠(yuǎn)，才能保證沖擊力達(dá)到穩(wěn)態(tài)時不受邊界影響。經(jīng)過試算，確定與上下邊界的距離為8D。為減弱初始計算的不穩(wěn)定問題，假定流滑體與樁之間的初始距離為0.01 m，以模擬二者在計算開始時即將接觸。

圖1 流滑體沖擊單樁俯視圖Fig.1 Plan view of sliding impact on a pile

樁身附近的網(wǎng)格劃分如圖2所示，樁周0.2D范圍為加密區(qū)，采用四邊形單元，其它區(qū)域為三角形單元。樁身附近單元必須足夠小才能保證計算結(jié)果的收斂性，單元大小的確定將隨后討論。采用中心差分法顯式求解瞬態(tài)流動方程，為保證計算精度，規(guī)定時間增量步為：

圖2 圍繞樁身的細(xì)網(wǎng)格Fig.2 Fine mesh around the pile

(4)

式中:Lmin為單元最小長度；vmax為節(jié)點(diǎn)速度的最大值。

1.2 冪律流體及其雷諾數(shù)

(5)

式中:n為冪律指數(shù)，當(dāng)n=1時，式(5)退化為牛頓流體，高塑性無機(jī)黏土的n=0.06～0.17[18]；K為稠度系數(shù)，對于黏土質(zhì)流滑體，K的范圍在10～1 300 Pa·sn之間[11]。參照滑坡沖擊管線問題[9]，引入名義應(yīng)變速率：

(6)

對應(yīng)的名義剪應(yīng)力為：

(7)

牛頓流體中雷諾數(shù)定義為：

(8)

(9)

2 數(shù)值模型驗證

2.1 網(wǎng)格敏感性

對于圖1所示的流滑體沖擊單樁，首先進(jìn)行試算，以確定合適的網(wǎng)格密度。以典型冪律流體參數(shù)為例：n=0.1，K=330 Pa·s0.1，ρ=1 600 kg/m3，沖擊速度v=2 m/s，樁徑D=0.5 m。采用三種網(wǎng)格劃分方案：緊鄰樁身的單元尺寸分別取0.005D，0.01D和0.02D，沿徑向的加密區(qū)內(nèi)單元尺寸按照1.05倍的比例逐漸增大(見圖2)，計算域邊界上的最大單元尺寸分別為0.1D，0.2D和0.4D。得到歸一化的沖擊壓力p/τapp隨時間t的變化如圖3。沖擊壓力達(dá)到峰值后迅速下降至穩(wěn)態(tài)值，三種網(wǎng)格劃分對應(yīng)的穩(wěn)態(tài)沖擊壓力差別在3%以內(nèi)，這表明樁周圍采用0.01D的單元尺寸足以消除網(wǎng)格敏感性。更多典型流體參數(shù)的計算也得到類似規(guī)律。為平衡計算精度與耗時，在后面的模擬中統(tǒng)一取樁周單元尺寸為0.01D，計算域邊界上的最大單元尺寸不超過0.2D。

圖3 最小單元尺寸為0.02D(綠色)，0.01D(紅色)和0.005D(黑色)的模型計算結(jié)果Fig.3 Resultsofmodels with minimum element size varied 0.02D (Green),0.01D (Red)and 0.005D (Black)

2.2 流滑體沖擊半埋管線問題

圖4 滑坡沖擊半埋管線問題Fig.4 Problemof semi-buried pipeline impacted by landslide

在Fluent中建立二維模型，圖5比較了所得沖擊壓力時程曲線與文獻(xiàn)[14]物質(zhì)點(diǎn)法結(jié)果。兩種方法差別在5%以內(nèi)，初步證明了Fluent模擬非牛頓流體沖擊樁身是可行的。

圖5 Fluent與物質(zhì)點(diǎn)法預(yù)測的沖擊管線壓力Fig.5 Normalised impact pressures(on the pipeline)predicted by Fluent and the material point method

3 結(jié)果和分析

3.1 沖擊過程

采用建立的Fluent數(shù)值模型，考察冪律型非牛頓流滑體對樁身的沖擊。如圖3所示，沖擊壓力在0.015 s內(nèi)增大到峰值(OA段)，然后由峰值迅速減小(AB段)，0.05 s之后變化速率放緩(BC段)，0.15 s后沖擊壓力基本達(dá)到穩(wěn)定值(CD段)。各典型時刻的流滑體形態(tài)如圖6所示：沖擊初始階段(t=0.02 s)，流滑體沿著樁表面流動，然后流滑體在樁后展開(t=0.035 s)，達(dá)到穩(wěn)態(tài)沖擊壓力時(t=0.2 s)流滑體與樁背側(cè)保持脫離。因峰值持續(xù)時間較短，對樁身威脅更大的是持續(xù)時間較長的穩(wěn)態(tài)沖擊，以下將集中討論穩(wěn)態(tài)沖擊壓力。

圖6 沖擊過程中的流滑體形態(tài)Fig.6 Morphological evolution of slide mass during impacting

3.2 變動參數(shù)分析

改變流滑體的密度和強(qiáng)度特性、沖擊速度及樁徑，實(shí)施變動參數(shù)分析，探索穩(wěn)態(tài)沖擊壓力與上述參數(shù)之間的關(guān)系。如表1所示，根據(jù)海底滑坡的實(shí)際條件，選擇ρ=1 300、1 600、1 800、2 000 kg/m3；n=0.08、0.1、0.15；模型表達(dá)式分別為τ=3100.08、τ=3300.1、和τ=3800.15Pa；v=0.1～6 m/s；D=0.5、1、2、4 m。表1中的48個算例涵蓋了Re=0.1～140之間的范圍。

表1 計算模型參數(shù)取值Table 1 Variable parameters of calculation models

采用傳統(tǒng)形式的單純拖曳力公式(1)，總結(jié)表1中所有算例的Fluent模擬結(jié)果，拖曳力系數(shù)與雷諾數(shù)關(guān)系擬合為CD=0.95+12.4/Re1.1，見圖7(a)、(b)更清楚地展示了低雷諾數(shù)下數(shù)值計算與擬合的CD?？梢钥闯?，當(dāng)CD-Re曲線的坡度趨于平緩即Re>5時，擬合公式效果良好；但當(dāng)Re<2時，如圖7(b)所示，CD-Re曲線的坡度急劇變化，數(shù)值模擬與擬合公式所得結(jié)果的相對誤差可能高達(dá)20%～50%，這表明低雷諾數(shù)的慢速滑動沖擊問題不適合采用單純拖曳力表達(dá)形式。

((a)寬廣雷諾數(shù)范圍；(b)低雷諾數(shù)范圍。(a)Wide range of Reynolds numbers;(b)Low values of Reynolds numbers.)圖7 采用單純拖曳力公式總結(jié)的數(shù)值模擬結(jié)果Fig.7 Numerical results fitted as drag-force equation under

為解決低雷諾數(shù)下的沖擊問題，需要同時考慮流滑體的擬靜力不排水強(qiáng)度和慣性拖曳力。參考公式(3)，擬合表1算例的沖擊壓力值，見圖8(a)，得到NC=7.2和CD=0.76，即

(10)

低雷諾數(shù)下的沖擊壓力與雷諾數(shù)關(guān)系如圖8(b)所示。當(dāng)Re<4時，F(xiàn)luent模擬給出的p/τapp趨向于定值，也就是與公式(3)中的NC相關(guān)。圖8表明，公式(10)能夠合理預(yù)測Re=0.1～140范圍內(nèi)流滑體對樁身的沖擊壓力。與公式(1)相比，公式(3)不僅提供了更高的擬合精度，而且具有清晰的物理意義，低雷諾數(shù)(低流速)時的沖擊力實(shí)質(zhì)上就是土力學(xué)擬靜力條件下的地基承載力。

((a)寬廣雷諾數(shù)范圍；(b)低雷諾數(shù)范圍。(a)Wide range of Reynolds numbers;(b)Low values of Reynolds numbers.)圖8 歸一化沖擊壓力與雷諾數(shù)擬合曲線Fig.8 Fitting curve of normalised impact pressure

式(10)給出了沖擊壓力中拖曳力項和承載力項各自的貢獻(xiàn)，每一項所占比重與流滑體速度和名義剪應(yīng)力相關(guān)。如圖9所示，當(dāng)Re<10時，承載力占據(jù)主導(dǎo)地位，尤其是Re<5時，拖曳力占比超過70%；隨著Re增大，拖曳力在沖擊壓力中的占比增大；當(dāng)Re>75時，拖曳力占比大于80%，此時可以用單純拖曳力形式的式(1)估計沖擊壓力。

圖9 承載力項和拖曳力項在歸一化沖擊壓力中的占比Fig.9 Proportion of bearing capacity and drag force terms in normalized impact pressure

4 結(jié)語

將近海滑坡流滑體視為冪律型非牛頓流體，采用流體力學(xué)軟件Fluent，探索流滑體對樁身的沖擊?？偨Y(jié)了中低雷諾數(shù)下，流滑體的變形形態(tài)和基礎(chǔ)承受的沖擊力變化趨勢。給出了穩(wěn)態(tài)沖擊壓力表達(dá)式。針對流滑體密度、粘度和沖擊速度以及樁體直徑等影響因素，建議采用同時包含承載力項與拖曳力項的穩(wěn)態(tài)沖擊壓力表達(dá)式，即公式(3)，并給出相應(yīng)系數(shù)，可以覆蓋雷諾數(shù)在0.1～140范圍內(nèi)的工況。

當(dāng)雷諾數(shù)較小時，樁身承受的沖擊壓力主要來自流滑體的不排水抗剪強(qiáng)度，可以通過承載力項估計沖擊力。大雷諾數(shù)情況下，拖曳力項的貢獻(xiàn)將趨于主導(dǎo)，拖曳力項可以代表沖擊力的數(shù)量級。