李宛玲，張 琪，周香蓮

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2.中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071；3.上海交通大學 土木工程系，上海 200240)

近年來,隨著沿海城市人口總數(shù)快速增長，沿海地區(qū)電力需求增加.相比陸上風電機組，因具有不占用土地面積、風能資源豐富、綠色清潔、低碳排放、風速高且湍流強度小等優(yōu)點，海上風力渦輪機在世界范圍內得到廣泛的應用和發(fā)展[1-2].由于海上風電機基礎長期受風荷載、波浪荷載及其他荷載的共同作用，威脅到風機工程的結構安全，所以，海上風機的受力及影響等相關問題引起了人們的廣泛關注[3-5].

在海洋環(huán)境中，風荷載往往會影響波浪的水力特性，很多學者曾對風在海洋環(huán)境中的作用進行研究[6-8].Phillips[9]提出風荷載通過在水面上產(chǎn)生剪切力的理論來探究風浪間的相互作用,該理論雖然解釋了海洋環(huán)境中風荷載的產(chǎn)生和作用機制，但并未對風荷載作用后的效應做進一步說明.Liu等[10]用演化方程和數(shù)值模型相結合的方法分析了風對Boussinesq波浪的非線性影響，推導出在風荷載作用與理想的水平海底條件下Boussinesq波浪的演化方程.文獻[11]將試驗結果與數(shù)值模擬對比后發(fā)現(xiàn)風荷載會加速波浪破碎，延長破碎時間并增大波浪的最大振幅，且這種影響在深水區(qū)更加穩(wěn)定.

單樁是最常見的海上風機基礎結構，在多種荷載共同作用下，海上風機基礎會產(chǎn)生一定程度的水平位移與彎矩.許多學者通過理論分析、試驗模擬以及數(shù)值模擬等方法對海上風機在波浪荷載作用下的動態(tài)響應進行了研究[12-15].例如：Arany等[16]將現(xiàn)場試驗中的渦輪機數(shù)據(jù)轉換為彎矩矩陣，評估風荷載、波浪荷載、轉動頻率及風機葉片運動頻率這4個載荷的相對大小，建立了交界面處的彎矩譜公式來預測發(fā)生疲勞損壞的位置.Abhinav等[17]討論了在波浪荷載的作用下樁-土作用對海上風機動態(tài)響應的影響，研究發(fā)現(xiàn)，忽略樁-土作用會導致海上風機的極限強度高出3%～60%.需要說明的是，以上研究均局限于將風與波浪作為獨立等效荷載作用于海床與樁基礎，忽略了風-浪-樁之間的連續(xù)作用.

總之，在海洋環(huán)境下樁基動力問題的研究上，國內外學者已進行了一定的研究，但海洋環(huán)境荷載作用下的樁基響應機理和規(guī)律至今尚未得到充分揭示，因此有必要對風-浪-樁的連續(xù)相互作用展開研究.本文采用有限元方法建立風浪-海床-單樁三維數(shù)值模型，通過討論不同的風與波浪特性參數(shù)影響，研究了風-浪相互作用以及單樁基礎在以上兩種荷載共同作用下產(chǎn)生的水平位移與彎矩，分析風浪荷載在海洋工程中的作用與影響.

1 數(shù)值模型與計算方法

圖1所示為風浪-海床-單樁數(shù)值模型示意圖.圖中：H為波高；L為波長；D為樁的直徑；d為水深；h為海床深度；l為樁的埋入深度.

圖1 風浪-海床-單樁數(shù)值模型示意圖

1.1 控制方程

1.1.1風-浪荷載控制方程 風-浪模型將流體體積(VOF)法與部分面積體積障礙模擬(FAVOR)技術相結合，描述波浪模型的自由液面.為了探究風荷載與波浪荷載共同作用下的水動力特性，采用風荷載作用下的波浪動量方程[18]：

(1)

式中：u、v及w分別為流體在x、y及z方向的流速；t為時間；VF為可流動流體的體積分數(shù);Ai(i=x,y,z)為i方向的可流動流體面積分數(shù);ρ為流體密度;p為孔隙水壓力;gi(i=x,y,z)為i方向上的重力加速度;τs,x、τs,y及τs,z分別為x、y及z方向上風在流體上作用的切應力；τb,x、τb,y及τb,z分別為x、y及z方向上波浪模型底部作用的切應力；fi(i=x,y,z)為i方向上流體的黏滯力加速度.τs,x、τs,y及τs,z的表達式如下[9]：

(2)

為了消除出波邊界產(chǎn)生的反射波的影響，在出波邊界前設置了波浪吸收層，上式修正為

ρksp(〈ufi〉-〈ufi〉str)

ρksp(〈ufi〉-〈ufi〉str)

ρksp(〈ufi〉-〈ufi〉str)

式中：ksp(〈ufi〉-〈ufi〉str)為波浪吸收層的人造阻力;ksp為以s-1為單位的阻尼系數(shù)；ufi為波浪吸收層的流體速度；〈ufi〉str為流體原始流速.

1.1.2海床控制方程 采用Biot方程的部分動力方程“u-p方程”作為海床土的控制方程，假設海床為各向同性的彈性多孔介質，各個方向的滲透系數(shù)相同，土骨架的應力應變服從胡克定律，孔隙流體為可壓縮流體且滿足達西定律.連續(xù)性方程和平衡方程可表示為[18]

(3)

(4)

(5)

(6)

式中：Kw為孔隙水的體積模量(一般取Kw=2×109N/m2);Sr為海床土飽和度；pw0為孔隙水的絕對壓力.

海床土的有效正應力和切應力可表示為

(7)

式中：Gs為海床土的切變模量;μ為海床土的泊松比.

1.2 數(shù)值模型與邊界條件

1.2.1數(shù)值模型 為考慮風-浪-樁的相互作用影響，在Flow3D中建立風-浪模型，該模型包含風荷載、波浪以及水-土交界面以上的單樁，整個流體區(qū)域長度為160 m，寬度為40 m，高度為25 m，在經(jīng)過不同尺寸規(guī)格流體域的測試和驗證后，當前尺寸既能保證計算結果的精度，消除邊界效應的影響，同時也能提高計算效率.波浪吸收層的長度為80 m，單樁位于坐標原點處，距離入流邊界40 m.從風-浪模型中得到樁周波形的變化情況與流體傳播速度，樁身與海床所受動態(tài)水壓，將計算水壓作為外部荷載加載到海床-樁基模型上.海床模型的控制方程和邊界條件是由COMSOL Multiphysics軟件通過偏微分方程模塊設置，長度為80 m，寬度為40 m，高度為75 m，包含 262 079 個自由度和 31 106 個單位，最大網(wǎng)格單元尺寸為2 m，并且將單樁基礎附近的網(wǎng)格局部細化至最大單元網(wǎng)格尺寸為1.5 m，如圖2所示，圖中g為重力加速度.

圖2 數(shù)值模型網(wǎng)格

1.2.2造波模型邊界條件 風-浪模型的左側為入流邊界,在入流邊界進行線性波的造波；右側為設有波浪吸收層的出流邊界；模型前后兩側設置為對稱邊界，即在該邊界法線方向上流體可以沿該邊界自由運動而不向外滲透和溢出；模型底部的水-土交界面設置為壁面邊界；上部的水-空氣交界面設定為自由表面，其大氣壓力為標準大氣壓力；樁-水交界面設置為不透水邊界.

1.2.3海床模型邊界條件 在海床表面，將從風-浪模型中得到的波浪對海床土體的水壓力pb設置為海床的孔隙水壓力p，忽略流體的黏性和摩擦力，豎向有效正應力和切應力均為0，即在z=0處，

(8)

海床底部是不透水的剛性邊界，其海床位移為0且流體在此邊界無豎向流動，即在z=-h處，

(9)

(10)

(11)

單樁和海床的交界面為不透水邊界，即

(12)

式中：a0為風-浪模型的長；b0為風-浪模型的寬；n為樁表面沿樁徑方向.

2 數(shù)值模型驗證

本文所研究波浪荷載為線性波，因此可將通過本文數(shù)值造波得到的波形與解析解進行對比，以驗證所提出模型的準確性.線性波的波浪參數(shù)為：波高H為2 m，水深d為70 m，周期T為15 s，波長L為311.59 m.結果對比如圖3所示.由圖3可見，本文所采用造波模型波面高度z的結果與理論值較為吻合.

圖3 線性波理論值與本文數(shù)值造波模型波面高度結果對比

為驗證海床模型的準確性，將本文模型結果與文獻中海床相應解析解[19]進行對比，所得土體有效正應力值與波浪壓力值如圖4所示，圖中p0表示波浪作用于海床表面的壓力幅值.所采用的海床模型深度為25 m，長度為 1 000 m，泊松比為0.333, 孔隙率為0.3, 土體密度為 1 850 kg/m3, 切變模量為1×107N/m2, 滲透率為1×10-2m/s, 飽和度為1.從圖4可以看到，本文計算的海床模型結果與解析解結果基本一致，本文采用的u-p模型(考慮慣性力)與文獻[19]所使用Q-S模型(沒有考慮慣性力)的最大相對誤差為3.6%p0，證明本文模型所采用的數(shù)值模型對風荷載作用下波浪與海床樁基響應問題的研究結果是合理的.

圖4 本文數(shù)值解與文獻[19]解析解的對比

3 結果與分析

對風速、風剪切系數(shù)及波高等參數(shù)進行討論，分析風荷載對樁周流體變化、樁身水平位移與彎矩的影響.選取的數(shù)值模型計算參數(shù)見表1.

表1 本文采用數(shù)值模型的計算參數(shù)

3.1 樁周局部流體變化

圖5所示為當波峰傳播到樁前位置時(t=55.9 s)不同風浪參數(shù)下單樁周圍流體形態(tài)及波浪沿x方向流速vwave示意圖，圖中：W為風速.由圖5(a)、5(b)可以看出,在同一水深及波浪條件下，樁周流體的傳播速度隨風速與風剪切系數(shù)的增大而增大，且波浪在樁前產(chǎn)生的爬升更高，樁側波形變化明顯.該現(xiàn)象與文獻[9，20]的研究結果一致.當風速與風剪切系數(shù)增大時作用在流體上的風剪切力增大，風浪相互作用時流體能量增強，使得波浪加速向前傳播，更多的動能轉化為勢能，造成波浪沿樁基的爬升幅度增大(尤其在樁的迎浪面附近)，因而對樁造成更強的沖擊力.圖5(c)中，在相同風荷載條件下，由于波高增加，風荷載與波浪、波浪與樁身的作用面積擴大，樁身表面受到的動水壓力幅值上升，因而樁周流體流速幅值也隨之增大，加劇樁基周圍波面的非線性擾動和變形，甚至導致波浪破碎的形成.

圖5 不同的風浪參數(shù)下樁周流體波形與在x方向的流速

3.2 樁身水平位移

圖6 所示為波峰傳播到樁前時，不同風浪參數(shù)下單樁基礎樁身水平位移up的豎向分布以及風浪參數(shù)對up的影響規(guī)律.圖6反映了樁身水平位移隨風速、風剪切系數(shù)和波高增大明顯增大的情況，同時還給出了up在樁頂(z=25 m)、水-土交界面(z=0)、海床中(z=-25 m)及樁底(z=-50 m)4個位置隨風浪參數(shù)的變化規(guī)律.需要指出的是，為了可以更直觀地表現(xiàn)出參數(shù)變化對樁身水平位移的影響，圖中將產(chǎn)生的樁身水平位移按比例擴大為實際的 1 500 倍表達.圖中海床以上沿樁身水平位移近似線性增大，這是由于單樁與模型土體的彈性模量相比遠大于后者，所以單樁基礎近似剛性.由圖6(a)與圖6(c)可以看出，隨著風速與波高增加，樁身同一高度處的水平位移逐漸增大，但是樁身水平位移的增大速率越來越小，即風速增大對樁身的水平位移響應影響減小.而當風剪切系數(shù)增大時，水平位移增大的同時增長速率越來越大，說明風剪切系數(shù)增大對樁身的水平位移響應影響更加明顯.這與上文中對流體的響應分析現(xiàn)象一致.受風荷載影響后的波浪荷載在樁基結構上的作用力變化，進而影響樁身的水平位移.表2所示為不同條件下樁身最大水平位移ux,max，表中ux0,max(標準位移)為風速為15 m/s，風剪切系數(shù)為0.5，波高為3 m時樁身的最大水平位移.因此，在結構設計過程中必須考慮風浪荷載的共同作用對樁基水平位移的影響.

表2 不同條件下ux,max和ux,max/ux0,max

圖6 up的豎向分布以及風浪參數(shù)對up的影響

3.3 樁身水平受力

圖7所示為風浪荷載共同作用下考慮不同風浪參數(shù)時沿樁身彎矩Mx的豎向分布曲線，圖8所示為風速、風剪切系數(shù)及波高對Mx的影響規(guī)律.由圖7和8可以看出，隨著風速、風剪切系數(shù)及波高的增加，樁身同一位置所受彎矩增加，且由于土抗力的存在，彎矩的最大值均發(fā)生在泥面位置以下，同時風速與波高增大對樁身所受彎矩的影響趨勢變小.風剪切系數(shù)增大時樁身所受彎矩的增大速率越來越大，即影響趨勢變大.表3所示為不同風速、風剪切系數(shù)及波高時樁身最大水平彎矩Mx,max，表中Mx0,max(標準彎矩)為風速為15 m/s，風剪切系數(shù)為0.5，波高為3 m時樁身的最大水平彎矩.樁身彎矩結果與上文中的風荷載對波浪的作用、風浪共同作用下的樁基水平位移結論一致，可見風浪荷載共同作用對海上風電機樁基荷載影響較大，所以優(yōu)化過程中需要考慮風荷載對波浪的作用以及以上兩種荷載共同作用下整個樁基荷載的影響，防止設計荷載偏小.

圖7 不同風和波浪下Mx的豎向分布

圖8 不同風和波浪參數(shù)對Mx的影響

表3 不同條件下Mx,max和Mx,max/Mx0,max

4 結論

本文在考慮風荷載影響的情況下建立了風浪-海床-單樁三維單向耦合數(shù)值模型，研究海洋環(huán)境中風浪荷載共同作用下樁周流體的變化及其對樁身水平變形和受力特性的影響，得到以下結論：

(1)對風荷載作用下樁周流體的分析結果表明，風荷載的存在會加劇樁周的流場擾動.風速、風剪切系數(shù)及波高的增大會加快波浪流速，使波浪在到達樁前時有更多的動能可以轉化為勢能，從而增大波浪在樁前產(chǎn)生的爬升幅度，單樁周圍波面的非線性擾動與變形加劇，甚至導致波浪破碎.

(2)對風浪荷載作用下樁基響應的計算結果表明，風荷載使波浪作用下的樁基響應進一步增大.風速、風剪切系數(shù)與波高的增大會加劇樁周的波浪擾動，進一步增大樁基樁身的水平位移與彎矩.樁身水平位移與彎矩的增幅隨風速與波高的增大而衰減，隨風剪切系數(shù)的增大而增大.

(3)在海洋環(huán)境中，風荷載對波浪的傳播與樁基的響應均有較大的影響.因此在設計海上風機結構計算樁基承載力的過程中，應在考慮風荷載對樁周流場作用的基礎上，綜合考慮風浪荷載共同作用對樁基基礎的影響.