水流荷載作用下海洋群樁受力特性研究

時 健，蔣宗南，孔德森,3,*

(1.山東科技大學 土木工程與建筑學院, 青島 266590；2.北京中巖大地科技股份有限公司, 北京 100089；3.山東科技大學 山東省土木工程防災減災重點實驗室, 青島 266590)

近海風電基礎多采用大直徑單樁和群樁的形式，其中大直徑單樁工藝成熟，使用廣泛，但會有溜樁[1]等問題的存在，導致施工困難.而群樁雖然施工方便，但在海洋環境中的受力情況十分復雜，有上部風機自重荷載、結構自重荷載、風荷載、波浪荷載、水流荷載等，其中群樁的水流荷載由于其結構的特殊性，受力十分復雜.對于單樁結構在特定的雷諾數范圍內，會在樁后產生左右兩側成對、交替出現、方向相反的漩渦，即卡門渦街現象[2].

在計算流體力學領域，針對圓柱繞流的研究，主要以雷諾數為控制參數，研究單圓柱的漩渦脫落情況[3-5].近年來，針對高雷諾數[6]和多圓柱繞流[7-9]的研究開始增多.現在，關于變截面圓柱繞流的研究[10]也逐漸展開.

但海洋環境中雷諾數都比較高，并且群樁的排列方式多種多樣，水動力特性極為復雜，對于2×2群樁在海洋環境下的水流荷載受力特性研究并不深入.

程友良等[11]采用大渦模型，對雷諾數Re=3.5×106的情況下，對多個比轉速的旋轉圓柱繞流進行了數值模擬，并分析了圓柱的升、阻力系數和斯特勞哈爾數(St).結果表明，隨著比轉速的增大，樁柱的橫向力逐漸增大，阻力系數也隨之增大，而升力系數呈周期性減弱.

崔宜梁等[12]對串列雙圓柱進行三維數值模擬，分析了雙圓柱受到的平均阻力、升力，后柱周向壓力以及斯特勞哈爾數等水動力特性.結果表明，在Re=20 000的串列雙圓柱繞流中，兩圓柱中心間距L與圓柱直徑D的比值為4時，后柱受前柱尾流影響大.

李聰洲等[13]采用改進的延遲分離漩渦方法模擬了高雷諾數下的樁體繞流，包括單圓柱、單方柱、串列雙圓柱和串列雙方柱四種情況，研究了不同雷諾數下圓柱與方柱的水動力特性.研究表明，單方柱繞流在2000

綜上所述，為較好地反映群樁結構在海洋洋流荷載作用下的受力和振動特性，本文采用Fluent有限元軟件，對海洋群樁進行建模和計算，分析各樁的升力系數和阻力系數的時程曲線以及樁周環向壓力分布，以期得到海洋群樁受力分布規律和振動特性，對海洋工程中群樁基礎的設計和優化提供參考.

1 數值計算模型的建立

1.1 模型尺寸

海洋群樁的工程環境為四面廣闊的流場，而數值模型無法完全模擬，因此模型假設流場為一維流場，從左側進口邊界流入，在右側出口邊界流出.同時，為了模擬真實的海洋環境，將上下邊界設為對稱邊界.海洋群樁數值計算模型如圖1所示.

圖1 海洋群樁數值計算模型(單位：m)

目前，海洋石油鉆井平臺以及大型風電機組都已從淺海走向深海[14].在海洋平臺的基礎中，四樁導管架的應用很廣泛，故進行四樁結構的模擬計算.樁直徑D為0.5 m，樁間距為2D.為避免邊界對流場產生影響，選擇進口邊界距前排樁15D，出口邊界距后排樁25D，上下邊界距離各側樁10D.流場進出口總長22 m，上下邊界相距12 m，樁長12 m.為便于分析，對各樁體進行編號，如圖1所示.各樁對應的數據用樁體編號作為下標來區分.

1.2 控制方程

流體運動需要遵守質量守恒定律.對于流場范圍內任意選定的幾何空間，單位時間內在此空間的流體質量保持守恒.因此，流體的速度及密度都是時間和空間坐標的函數，流體的運動屬于不可壓縮黏性流體的運動，其控制方程為連續性方程，即

(1)

式中：u，v，w為3個速度分量；ρ為流體密度，g/cm3.

1.3 湍流模型

湍流流動是一種高度非線性的復雜流動，處在湍流狀態下的流體，流場的各種參數如水質點的速度、壓力等都隨時間與空間發生不規律的變化.針對湍流流動的特點，按產生漩渦的尺寸，將流場分為大尺寸漩渦和小尺寸漩渦.

大渦模型(LES)的基本思想是：大漩渦對流場的影響比較大，湍流擴散、質量、動量、能量等的變換和雷諾應力都是通過大漩渦實現的，而小漩渦主要起耗散作用.研究表明，在尺度足夠小的情況下，幾乎所有的漩渦都具有一定的相似性.

因此，將湍流中的大漩渦和小漩渦分開處理，大漩渦通過N-S方程直接求解計算，小漩渦通過亞格子模型建立與大漩渦之間的關系，對其進行模擬，而大小漩渦的區分是通過濾波函數進行的.目前較為常用的濾波函數主要有三種：盒式濾波函數、富式截斷濾波函數和高斯濾波函數[15].

(2)

式中：V為計算單元的體積；x為濾波后的空間坐標；x′為實際流場中的空間坐標.

用此理論來解決不可壓縮流體，過濾N-S方程，將得到以下方程[16]：

(3)

(4)

其中，亞格子尺度應力τij為

(5)

以上公式為N-S方程的張良表示法，其中u，v為速度分量；p為壓力；上劃線符號標記的量代表該量為時均值；系數i，j的取值范圍為(1，2，3).

亞格子模型采用壁掛式本地渦黏度(WALE)模型[17].在Fluent中，默認的WALE常數Cω為0.325，研究表明將它用于流量范圍跨度大的區域能取得準確的結果.

1.4 初始條件與材料屬性

模型進口邊界采用inlet屬性，進口流速為1 m/s，進口壓力為0 Pa；出口邊界采用outflow屬性；上下邊界采用symmetry屬性，用于模擬無限流域；樁邊界采用wall屬性，即無滑移固壁邊界.

模型按照15 ℃時的海水屬性賦值，流體選用液體密度為998.2 kg/m3，動力黏度為0.001 021 Pa·s.

模型的雷諾數[18]為

(6)

式中：U為來流速度；L為特征長度；μ為動力黏度.

斯特勞哈爾數為

(7)

漩渦的脫落頻率f[19]為

(8)

漩渦的脫落周期t為

(9)

根據漩渦脫落周期，設置時間分析步為0.05 s，時長25 s.

1.5 網格劃分

采用ANSYS ICEM CFD軟件對模型進行網格劃分.將模型分為10個計算區域，從外到內逐漸加密，總體網格尺寸為0.2 m.由于流場模擬的雷諾數Re=488 834，為湍流流動，在采用大渦模型計算時要借助亞格子模型進行計算，故需要對樁附近的網格進行加密.加密區域分為內部區域和外部區域，外部區域寬度為2D，網格尺寸為0.05 m；內部區域寬度為1D，采用等比加密，第一層網格節點到樁體表面的距離為0.005 m，按1.076倍尺寸過渡.網格類型采用結構性網格，網格總數量為28 714.網格劃分及樁附近網格加密如圖2所示.

圖2 網格劃分及細部放大

檢測網格質量的參數y+值的計算公式[20]如下：

(10)

(11)

(12)

式中：Δy為第一層網格節點到壁面的距離；u*為近壁面摩擦速度；τω為壁面剪切應力；Cf為壁面摩擦系數.

根據計算，模型中各個壁面中y+值最大為197，壁面各網格點y+值隨網格點距壁面上游邊界最小值點順流方向距離的分布規律如圖3所示.研究表明，只要y+值在300以內，壁面函數就可以準確地計算出結果.

圖3 壁面各網格點y+值分布

2 模型驗證

模擬結果的準確程度由參數的選取和網格質量兩方面決定，而網格的質量又由網格的邊界層劃分、網格密度、光滑性以及單元類型決定.計算得到的升、阻力系數是否準確，漩渦的脫落和泄放是否能真實的模擬，很大程度上取決于網格質量的好壞.為了驗證模型的參數選取和網格劃分的準確性，以及與群樁模型做對比，建立了相同尺寸、相同參數的單樁模型，并采用試驗對比和斯特勞哈爾數驗算兩種方法進行分析.單樁模型計算結果如圖4所示.

圖4 單樁阻力系數和升力系數時程曲線

對比算例采用了吳國雄教授Re=500圓柱繞流試驗[21]，試驗結果如圖5所示.由于試驗無法實現與實際工程一樣大小的足尺模型，所以試驗環境下的雷諾數都相對較小.故采用文獻[22]中的弗勞德相似原理進行縮尺對比.計算結果表明，單樁模型結果與吳國雄教授試驗結果基本一致.

圖5 單圓柱繞流水動力系數試驗結果[21]

斯特勞哈爾數(St)是研究圓柱繞流的一個重要特征參數，可以衡量數值模擬結果準確性，是流體力學中相似準數的一個，為量綱為1常數[23].

斯特勞哈爾數的計算式為

(13)

式中：fs為漩渦脫落頻率.

如前文所述，穩定后升力系數的頻率即該樁漩渦脫落的頻率.因此，將樁體的升力系數曲線進行快速傅里葉變換(FFT)，生成功率譜密度分布，如圖6所示.

由圖6可以看出，1號樁與2號樁升力系數功率譜密度都在0.404 Hz時達到峰值.將漩渦脫落頻率fs=0.404 Hz代入式(13)，得到斯特勞哈爾數為0.202，與式(11)計算結果相近，說明數值模擬結果有效.

3 計算結果分析

3.1 阻力系數與升力系數分析

流體運動時對細長結構物的作用力主要有沿流向的阻力和垂直流向的升力.

當流體接近物體前緣時，流動受到阻礙而產生壓力，并且隨著流動，這個壓力會沿物體表面轉移到后方.但當Re數較大時，這一壓力會在物體截面最寬點產生分離，分離點即為物體表面流速從正到負的變速點，在分離點后面沿物體表面倒流，產生漩渦.

漩渦在樁體兩側交替產生和脫落，當樁體一側產生漩渦時，該側流速小于另一側，產生壓力差，也就是升力(FL).此漩渦脫落的同時，另一側漩渦也在產生，升力就因此改變方向.因此，漩渦脫落的周期就是升力時程曲線的周期，即升力周期.與地震荷載不同[24]，水流荷載通常是具有規則性的.

與此同時，漩渦的脫落會在樁體后面產生低壓區，與樁體前部流體產生壓力差，即阻力(FD).當漩渦脫落以后，低壓區壓力變大，阻力變小.阻力周期為升力周期的一半，即每個漩渦的產生和脫落，都會構成阻力系數的一個周期.

圓柱繞流過程中，樁受到來流方向平行和垂直的力，分別為阻力(FD)和升力(FL).其計算公式[25]為

(14)

(15)

式中：CD為阻力系數；CL為升力系數；U為無窮遠處的流速，m/s；D為樁徑，m.

在阻力和升力的計算公式中，密度、流速和樁徑都是定值，因此，阻力系數和升力系數的變化，即阻力和升力的變化.并且阻力系數和升力系數都是量綱為1系數，對樁的受力有代表性.

從圖7和圖8中各樁的阻力和升力時程曲線的對比中可以發現，1號樁和3號樁，即前排樁的受力相近；2號樁和4號樁的受力相近.說明2×2群樁中，面對來流方向的兩列樁會發生幾乎等效的相互影響.圖9所展示的渦量也很好地印證了這一點.

圖9 渦量等值線(t=25 s)

以1號和2號樁為例，對比分析群樁中各樁受力的不同.

1號樁的阻力系數穩定后維持在1.0～2.0，而2號樁阻力系數穩定后在-0.5～2.0，兩者到達波峰與波谷的時刻接近.

1號樁的阻力系數時程曲線十分具有代表性，其工況與單樁受力(圖4(a))相似.但2號樁的曲線與單樁的受力狀態很不一樣，其曲線的上限值與1號樁接近，但其穩定后的下限值卻遠低于1.0，甚至為負數，即阻力改變了方向.結合渦量等值線圖(圖9)分析可知，造成這一現象的原因在于1號樁后脫落的漩渦泄放至2號樁前時，會產生低壓區，其壓力甚至低于2號樁后產生的低壓區，故使2號樁曲線穩定后的下限值低于1號樁，甚至產生方向向前的阻力.

1號樁的升力系數穩定后以CL=0為對稱軸，上下對稱波動，幅值為1.5左右，與圖4(b)所示的單樁升力系數曲線比較接近，只是達到漩渦穩定釋放的時間，群樁要比單樁短；2號樁的升力系數穩定后以CL=0為對稱軸，上下對稱波動，幅值為3左右，兩者每次達到峰值的時間接近，波峰與波谷的位置也相近.

但2號樁的升力系數時程曲線每一個波峰的峰值都不相同，圖8(b)中t取25 s并不能完全揭示該現象的全貌，故取t=100 s進行計算，計算結果如圖10所示.

從圖10中可以看出，1號樁的升力系數時程曲線相對穩定，但2號樁曲線中各個峰值的大小波動較大，且上升周期與下降周期均不確定，只呈現出大致的波動趨勢，但其趨勢與1號樁曲線的峰值波動趨勢類似.同時，圖10還呈現出另一個現象，即2號樁所受的升力，負向峰值大于正向峰值，說明1號樁產生的漩渦大多經過2號樁的同一側，因此會減小2號樁該側的壓力，從而增大指向該側的升力.

3.2 樁周壓力分析

取24.4 s時，即1號樁升力系數CL1=0和單樁CL=0時，樁周環向的絕對壓力(取參考壓力為大氣壓101 325 Pa)如圖11所示.

由圖11可知，1號樁與單樁的樁周壓力的趨勢和數值都比較接近，大致呈現出沿0～180°線上下對稱，但1號樁明顯呈現出內側壓力大于外側的情況，經分析是由于群樁對流體的阻礙，降低了群樁范圍內的流體速度，增大了壓力，因此1號樁內側壓力大于外側.

相比較1號樁，2號樁的樁周壓力分布卻有很大的不同.2號樁的曲線也呈現出了內側壓力大于外側的規律，但峰值并沒出現在樁前(0°)和樁后(180°)的位置，而是出現在了樁前內測(300°)的位置，即表現出了峰值的滯后性，這一點與在圖9中所表現出的特征相互印證.

4 結論

1) 2×2群樁中，上下兩排樁的受力是上下對稱的，發生干擾的效果也是相互的.

2) 群樁中后排樁會受到前排樁所釋放的漩渦的影響，主要表現在后排樁漩渦釋放穩定后，阻力系數時程曲線的波谷極值對比前排樁有極大的降低，甚至出現負值，但波峰極值卻幾乎沒有影響.

3) 后排樁的升力系數對比前排樁，所受的升力波峰峰值不穩定，主要表現在2號樁的升力系數曲線波谷極值相比于1號樁增大了接近一倍.

4) 分析了群樁中各樁的樁周環向壓力，并與單樁進行對比，發現群樁的存在使其范圍內的流體流速降低，壓力增大，因此各樁會受到由群樁內部向外的壓力.