潮流作用下風電塔單樁基礎沖刷數值模擬

唐冬玥 趙 鳴

（同濟大學土木工程學院，上海200092）

0 引 言

風能已經成為當今最具規模化發展和商業化前景的可再生能源，海上風力發電受到了廣泛的關注。單樁基礎因其結構簡單且承載力高，是近海風電塔的主要基礎形式［1］。其主要通過上層土體傳遞側向力，但近海潮流往復流動會在單樁基礎周圍，形成沖刷坑，影響其動力特性［2］。

基于經驗方法和半經驗半理論方法，學者們建立了單向流下的基礎周圍最大沖刷深度公式［3］。國內主要使用的是公路工程水文勘測設計規范中建議的65-1 修正式和65-2 式［4］，國外比較常用的有美國聯邦公路署推薦的HEC-18 公式［5］和Sheppard-Melville方程（S-M方程）［6］等。

世界上首次提出潮流水流下橋墩局部沖刷的應用規范是2000 年的美國行業規范《潮汐河道沖刷》［7］。關于雙向流下海洋結構物基礎周圍沖刷最大深度，學者們主要是通過與單向流計算結果對比［8-10］，或者是基于橋墩沖刷深度測量，如韓海騫的潮流作用下橋墩局部沖刷深度計算公式［11］。

目前的研究中對樁-土-潮流相互作用的沖刷機理認識不清［12］，因此本文對雙向流下海上單樁基礎風電塔沖刷進行研究。在研究方法選擇上，雖然學者們可以通過實驗方法獲得可與規范對比的結果［13-14］，但是無法避免縮尺效應，相比之下數值模擬可以原型參數模擬［3］，現場觀測可以得到工程實際數據。因此本文采用了數值模擬與現場觀測相結合的方法進行研究，并通過參數分析推導雙向流下大直徑單樁基礎風電塔基礎周圍的局部沖刷深度公式。

1 理論模型

1.1 水沙接觸面模擬

當流動的水流受到樁基礎阻擋，在樁周圍形成湍流，帶走了樁基礎周圍的土體。隨著沖刷坑的發展，水沙接觸面的高度、形狀以及位置是不斷變化的。為追蹤多相混合體系中運動相之間界面，采用了自由表面追蹤法（VOF）［15］。對于結構化網格，引入流體體積分數VF，即：對于一個網格單元，若VF在（0，1）范圍內，則表示自由液面所在網格單元。

1.2 湍流模型

重整化群湍流模型（RNG）可以有效模擬旋轉均勻剪切流，帶分離的流動，適合模擬受樁阻擋后的水流繞流和渦旋，引入VF后的RNG湍流模型方程為

式中：xi為坐標軸方向，i=1，2，3為流體時均流速分量；Axi為xi方向上的可流動面積分數；ρ 為流體密度；k 為湍動能；ε 為湍動能耗散率；PT、GT為湍動能產生項；Dk、Dε為耗散項；C1、C2、C3為模型常數。

1.3 泥沙運動方程

沖刷過程中水流對泥沙的作用主要包括侵蝕、搬運和沉積三個部分。

1.3.1 泥沙起動

泥沙侵蝕即泥沙起動過程，采用時間和空間平均的臨界希爾茲切應力理論計算泥沙的起動，計算公式為［16］

式中：ulift為泥沙起動速度；ρ 為流體密度；ρs為泥沙密度；d 為粒徑參數；d*為無量綱參數；為考慮斜坡修正后的臨界希爾茲數［17］；α 取值為夾帶系數；θ 為局部希爾茲數，當局部希爾茲數超過臨界希爾茲數，泥沙起動出現侵蝕。

1.3.2 泥沙搬運

泥沙搬運分為推移質搬運和懸移質搬運兩種。

粒徑較粗的泥沙沿床面作滾動、滑動或躍移被稱為推移質。基于水流搬運泥沙所做的功與水流消耗的能量守恒的推移質輸沙模型，推移質體積輸沙率qb計算公式為

式中，Φ為無量綱的推移質輸沙率。

計算公式為［18］

粒徑較細的泥沙，由于紊動作用懸浮在水中并隨水流一起運動被稱為懸移質。懸移質泥沙濃度Cs可以根據泥沙擴散方程求解：

式中，D為擴散率。

從而根據濃度與流速乘積的垂線積分得到懸移質輸沙率。

1.3.3 泥沙沉積

當推移質停止搬運或者懸移質沉降則發生淤積，泥沙沉降速度usetting計算公式為［17］

式中：υf為流體運動黏度；ζ0為用于考慮粒子間相互作用Richardon-Zaki相關系數。

2 雙向流下單樁基礎周圍沖刷三維數值模擬

2.1 模型建立

基于上述沖刷理論方程，對江蘇龍源蔣家沙風電場中風電塔建立模型。該風電場主要由樁徑5.3 m 的4.0 MW 單樁基礎風電塔組成，所在海域大部分水深在8 m 以內，為正規半日潮海區，潮流流速強勁，垂向平均流速最大達到2.36m∕s。地勘報告顯示勘探深度范圍內上部粉土、粉砂，下部為沉積物。得到模型計算參數如表1所示。

海床由沉積沙組成，其中夾帶系數α 為0.018［16］，推移質系數β 取8［18］，初始泥沙厚度取8 m，流體高度為水深H，海床中掏空建立單樁，單樁為非變形實體類型。當受到單樁擾動的水流流速超過ulift，則泥沙侵蝕，沉積沙以推移質或懸移質搬運，當流速減小時，泥沙沉積。

對其建立結構化計算網格，以樁中心為原點，水流流動方向為X 方向，垂直流動方向為Y 方向，高度方向為Z 方向，三個方向計算范圍分別為106 m、72 m、18 m。為提高計算速度，根據結構對稱性采用半網格法。為提高計算精度，采用嵌套網格法對樁周圍3D范圍內網格加密，嵌套網格與外網格的尺寸比為1∶2。網格劃分如圖1所示。

表1 模型計算參數Table 1 Model calculation parameters

圖1 網格劃分Fig.1 Meshing block

水流初始速度為0，初始壓強為Z向的靜水壓強。水流的邊界條件如圖2 所示，其中流速邊界：設置X方向速度，其余方向速度為0；壓力邊界：各水流速度分量、紊動能和水壓力的法向梯度均為零；對稱邊界：各通量均為零，僅有切向信息。需要注意的是，當流動方向改變時，入口與出口邊界需要互換。

圖2 邊界條件Fig.2 Boundary conditions

周期用T 表示，每經過0.5T，流動方向改變。以2 小時40 分鐘為一個周期，取計算周期為2T，則水流變向3次，樁側參考點位分布如圖3所示。

圖3 樁周圍參考點位置示意圖Fig.3 Schematic diagram of the position of the reference points around the pile

2.2 仿真結果

取以樁為中心邊長為60 m 的方形泥面區域，則風電塔單樁基礎周圍從0 至2T 時間段內的沖刷坑形狀變化過程如圖4所示。

在0.5T 和1.5T 時，最大沖刷深度位于樁左側，在1T 和2T 時，最大沖刷深度位于樁右側，雙向流下最大沖刷深度的位置隨著流速方向的改變而改變。最大沖刷深度為3.16 m。

從水流變向后一段時間內沖刷坑的形狀變化看，可以發現，在雙向流流體變向后的一段時間內，沖刷的主要作用是重塑沖刷坑的形狀，轉移最大沖刷深度的位置。

對圖3 中樁側各參考點位置處的沖刷深度繪制其隨著時間的變化曲線，如圖5所示。

根據圖5可以看出，A至D的沖刷深度曲線與E 至H 點的沖刷深度曲線表現為兩股交織的曲線，A 至D 點和E 至H 點隨時間變化分別呈現相同的規律。0 至0.5T 和1T 至1.5T 時間范圍內，流動沿X 軸正向，因此A 至D 點的沖刷深度大幅增大，而E 至H 點由于堆積的影響，沖刷深度有減小，在0.5T 至1T 和1.5T 至2T 時間范圍內，流動沿X 軸負向流動，因此E 至H 點的沖刷深度大幅增大，而A至D點沖刷深度有所減小。

最大沖刷深度的位置隨流動方向改變而改變，0 至0.5T 和1T 至1.5T 時間范圍內，沖刷深度最大值位于C 點與D 點與樁圓心組成的扇形位于樁上游45°至70°范圍內。0.5T 至1T 以及1.5T 至2T 時間范圍內，沖刷深度最大值位于E 點與F 點與樁圓心組成的扇形位于樁上游45°至70°范圍內。

圖4 沖刷坑形狀變化過程（單位：m）Fig.4 Change process of scouring pit shape（Unit：m）

2.3 與實測結果對比

2019 年3 月31 日對該風電場中某單樁基礎風電塔測量其沖刷深度。首先記錄水位高度h1=7.60 m；其次記錄樁附近的水深h2=11.50 m，見圖6，由此得到風電塔附近的沖刷坑深度為

根據上一節的數值模擬計算結果顯示：在水深為8 m、水流平均流速1.89 m∕s 的雙向流作用下，直徑為5.3 m的單樁基礎周圍沖刷深度計算值為3.16 m；根據現場實測結果，5.3 m 樁徑的周圍沖刷深度的測量值為3.90 m。將在雙向流下的數值模擬與現場實測的沖刷深度相比較，兩者的誤差在可以接受的范圍內，因此本文對于近海風電塔在雙向流下的單樁基礎周圍沖刷深度的數值模擬方法可行。

圖5 樁周圍參考點沖刷深度的時間變化曲線Fig.5 Time variation curve of scouring depth at the reference points around the pile

圖6 測量單樁基礎風電塔的基礎沖刷深度Fig.6 Measuring the scouring depth around monopile foundation of one wind turbine

3 沖刷機理研究

以上一節工況作為工況1，分別改變變量流速、反向水深、樁徑和土顆粒粒徑數值，以及對于不改變流向的情況重新計算沖刷，其中區分正向和反向水深是考慮到退潮時水深小于漲潮水深，各工況的參數表見表2。

表2 各工況參數表Table 2 Parameters of multiple cases

各工況2T 對應時刻的沖刷坑計算結果見圖7。

將各工況以及兩個周期后的最大沖刷深度總結在表3 中，結合表2 可以得到，沖刷時長相同的情況下，流速增大，土顆粒粒徑減小、樁徑增大，最大沖刷深度增大，單向流的最大沖深大于雙向流下，減小反向水深后最大沖刷深度減小。

圖7 各工況沖刷坑形狀（單位：m）Fig.7 Scour pit shape under different working conditions（Unit：m）

表3 各工況最大沖刷深度表Table 3 Maximum scouring depth for multiple cases

4 雙向流下最大沖刷深度

將數值模擬得到的沖刷深度計算值與現行沖刷深度計算公式進行對比，包括了國內主要使用的65-1修正式和65-2式［4］，國外主要使用的Sheppard-Melville方程（S-M方程）［6］，計算潮流作用沖刷深度的韓海騫公式［11］，對比結果如表4所示。

表4 數值模擬與公式計算的最大沖刷深度比較Table 4 Comparison of maximum scouring depth calculated by numerical simulation and existing formula

從表4 可以看出，與數值模擬結果相比，各公式計算結果均偏大，其中65-1 修正式計算結果偏差最大，其次計算橋墩群樁基礎在潮流作用下的沖刷深度的韓海騫公式，相較之下，65-2公式結果與計算結果最為接近，平均誤差為17%，但是也可以發現該計算公式對各參數的改變不夠敏感。因此需要重新推導雙向流下風電塔樁基礎的沖刷深度計算公式。雙向流下沖刷深度S 主要取決于樁徑D、來流流速V、正向水深H1、反向水深H2，泥沙粒徑d50等，因此可以得出：

根據量綱和諧法則，可以得到：

利用已有的數值模擬對式（10）參數進行率定，為了便于進行相關分析，對式（10）兩邊取對數，得

求解多元線性方程組，可以得到雙向流下大直徑單樁基礎風電塔基礎周圍的局部沖刷深度公式：

5 結 論

本文基于蔣家沙的風電塔，結合數值模擬和現場觀測研究了雙向流下單樁基礎風電塔基礎沖刷。計算了不同水利條件下的雙向流下單樁基礎風電塔沖刷，得到如下結論：

（1）本文通過數值模擬得到的沖刷深度與現場實測結果吻合良好，這驗證了本文的雙向流下單樁基礎風電塔的沖刷深度的數值模擬方法。

（2）雙向流下最大沖刷深度的位置隨著流速方向的改變而改變，在流體變向后的一段時間內，沖刷的主要作用是重塑沖刷坑的形狀，轉移最大沖刷深度的位置。

（3）沖刷時長相同的情況下，流速增大，土顆粒粒徑減小、樁徑增大，最大沖刷深度增大。考慮雙向流的往返水深的不對稱性，若正向水深相同，減小反向水深后最大沖刷深度減小。

（4）基于數值模擬結果并根據量綱分析和多元線性方程求解得到了雙向流下大直徑單樁基礎風電塔基礎周圍的局部沖刷深度公式。