夏冰，苗得勝，張敏，劉懷西，葛文澎

（明陽智慧能源集團股份公司，廣東 中山 528400）

隨著國家政策對新能源行業(yè)的傾斜和風能行業(yè)的蓬勃發(fā)展，海上風電由于風資源充足，地形影響小，距離電網(wǎng)消納端近等特點，受到日益廣泛的關注，加之技術發(fā)展降低了海上風電的度電成本，海上風電場建設正在進入快車道。海上風電快速發(fā)展的同時，仍然面臨一些技術性問題有待解決，其中，環(huán)境載荷對風機基礎結構的沖擊作用，及其引發(fā)的風機結構安全性問題需要給予一定的重視。根據(jù)歐洲海上風電場的基礎沖刷實測數(shù)據(jù)，包括Barrow、Scroby Sands和Otzumer Balje在內(nèi)的多個風電場場址內(nèi)，部分機位的基礎沖刷深度達到或超過1.2倍樁徑[1]。對于我國南方部分省份的沿海風電場，由于位于太平洋熱帶氣旋和南海熱帶氣旋登陸中國的主要路徑上，極端天氣引起的惡劣海況對基礎附近海床的沖刷作用十分明顯，隨之形成的沖刷坑對當?shù)仫L電場的穩(wěn)定運行造成了嚴峻挑戰(zhàn)。張冬冬[2]、楊博[3]等人通過不同范圍和深度的沖刷坑對單樁風機結構進行了結構模態(tài)和響應模擬，發(fā)現(xiàn)沖刷坑的存在和發(fā)展會降低風機結構的自振頻率，增加塔頂結構的水平位移，從而額外增加基礎結構的負荷，降低結構整體的安全性。目前的風機基礎設計中，對于由海流、波浪等環(huán)境載荷導致的基礎沖刷問題，一般通過解析公式、縮比試驗和數(shù)值模擬等方法對沖刷深度和沖刷坑形態(tài)進行預估。其中，解析公式能夠?qū)崿F(xiàn)對沖刷深度的快速預測，但計算結果難以通過精度分析的方式進行衡量，且計算結果無法給出沖刷坑的具體形態(tài)，因而也難以預測沖刷坑的存在對基礎結構的動力學影響?？s比試驗能夠在給定的環(huán)境條件下，得到基礎沖刷的真實數(shù)據(jù)，同時能夠結合擬采用的防護措施進行防沖刷試驗，可信度較高。但縮比試驗所采用的環(huán)境條件較為單一，難以復現(xiàn)實際工程中的復雜流動環(huán)境，且由于試驗條件的限制，實驗結果往往受到邊界條件的影響，因此縮比試驗結果與工程實測數(shù)據(jù)仍然存在一定的差異。數(shù)值模擬能夠克服縮比試驗中的場地條件限制，通過數(shù)學模型實現(xiàn)復雜環(huán)境條件下的沖刷模擬，計算數(shù)據(jù)的可視化程度高，且計算結果具有較高的精度，因此能夠彌補解析公式和縮比試驗的不足。數(shù)值方法目前被廣泛用于對沖刷問題開展細致地研究，以了解基礎周圍的流動環(huán)境、沖刷產(chǎn)生的機理和沖刷坑的發(fā)展規(guī)律，進而提出更合適的沖刷預測方法和防沖刷方案。

國內(nèi)外對于樁基沖刷的數(shù)值研究開展較晚，主要集中于圓柱形單樁基礎的沖刷模擬。Brilay和McDonald[4]首先通過三維不可壓縮N-S方程研究了硬化床面上樁基礎周圍的流動結構。Olsen和Melaaen等人[5]使用穩(wěn)態(tài)流動方程研究了無粘沙粒床面上，圓柱單樁的清水沖刷問題，床面切應力使用k－ε湍流模型求解，沖刷計算則采用泥沙輸運方程。Roulund等人[6]實施了穩(wěn)流沖刷水槽實驗，在穩(wěn)流條件下開展了硬化海床上圓柱形樁的定常和非定常動床沖刷數(shù)值模擬，對比了兩種方法的計算結果；并將沙粒床面上的單樁基礎沖刷數(shù)值模擬結果與實驗數(shù)據(jù)進行了對比，在沖刷坑上游得到了較好的一致性，但數(shù)值模擬時間較長，需要在計算過程中逐步調(diào)整迭代時間步。Baykal[7]等人研究了波浪運動對單樁基礎沖刷和海床回填作用的影響，在計算方法上，針對Roulund方法進行了改進，加入了懸移質(zhì)沖刷算法，改善了樁后沖刷深度計算結果與實驗值的一致性，但同樣存在數(shù)值模擬時間較長的問題。

國內(nèi)研究學者中，陳兵[8]等人使用切割單元法對海底管道沖刷進行了二維數(shù)值模擬，該方法在計算過程中不需要進行網(wǎng)格運動和重構操作，只需要施加合適的邊界近似。數(shù)值計算結果與實驗數(shù)據(jù)的對比顯示，切割單元法的模型準確度高，數(shù)值結果符合較好。李紹武[9]等人使用Flow-3D軟件對不同直徑的圓樁沖刷進行了數(shù)值模擬，總結了沖刷深度與樁徑的關系并研究了基礎近壁面網(wǎng)格寬度對計算結果的影響。齊梅蘭[10]等人通過三維數(shù)值模擬，研究了不同流速下，推移質(zhì)和懸移質(zhì)泥沙輸運對沖刷坑發(fā)展的影響，提高了樁下游沖刷深度的計算精度。

本文基于海上風機單樁基礎沖刷問題，在有限體積法求解框架下構建了流動-沖刷耦合計算模型。模型使用StarCCM+進行流動模擬，將床面切應力計算結果作為沖刷計算的輸入，在自編程序中進行床面沖刷計算和計算域內(nèi)網(wǎng)格的運動控制，隨后將更新后的計算域返回到流動計算中進行下一時間步的計算，從而通過流動和沖刷的交叉迭代實現(xiàn)耦合問題的實時模擬。為解決沖刷計算周期過長的問題，在沖刷計算中引入了加速算法，在保證求解過程魯棒性的基礎上大幅縮短了計算周期。本文創(chuàng)新點在于構建了便于工程應用的樁基礎沖刷數(shù)值模擬流程，在計算精度和資源消耗上實現(xiàn)了一定平衡，并通過算例驗證了流程的可靠性。

1 數(shù)值模型

1.1 流動控制方程

流動計算采用三維非定常不可壓縮N-S方程，湍流模型為SSTk-ω模型，控制方程表達式如下：

連續(xù)性方程：

動量方程：

湍動能方程：

耗散率方程：

式中：

ρ ?清水密度（kg/m3）；

u?流速矢量（m/s）；

μ ?水的動力粘度（kg/m·s）；

p+?包括靜壓在內(nèi)的總壓（Pa），p+=p+ρgh；

Pk=Gk+Gb+Gnl與Pω=Gω+Dω為結果項。湍流方程中各項的表達式和常數(shù)取值詳見StarCCM+軟件用戶說明書。

1.2 海床沖刷模型

沖刷運動主要受到床面切應力和泥沙粒徑的影響。沖刷過程的每個時間步內(nèi)，可以假設單元邊界的床面切應力與單元內(nèi)泥沙輸運量存在平衡關系，泥沙輸運處于守恒狀態(tài)。因此海床高程的變化可以通過泥沙輸運方程表示為：

式中：

n?泥沙孔隙率，取0.4；

i?坐標系內(nèi)的方向分量；

qbi?推移質(zhì)沿坐標系各方向的單寬輸沙率矢量，由Engelund和Freds?e公式得到：

式中：

μd?動摩擦系數(shù)，取0.51。

不難看出，單寬輸沙率qbi主要受到沙粒起動幾率pEF和平均輸沙速度矢量ub的 影響，其中pEF可表達為：

式中：

θ=uf2/(s?1)gd希爾茲數(shù)，表征靜止海床上推移質(zhì)抵抗剪切起動?的能力；

s= ρ√s/ρ ?推移質(zhì)與清水的密度比；

uf=基于壁面切應力 τ的海床摩擦速度（cm/s）。

任意床面坡度下的臨界希爾茲數(shù)可通過下式計算：

式中：

μs?靜摩擦系數(shù)，取0.63；

θc0?水平海床條件下的臨界希爾茲數(shù)，取0.05；

α?坡向和摩擦速度矢量的夾角（°），初始狀態(tài)下為0；

β?床面坡度（°）。

泥沙起動時，ub的大小可由下式得到：

泥沙起動后，α的大小發(fā)生變化，無法通過解析方式求解，此時需要通過牛頓迭代法求解床面泥沙平衡方程組求解相關參數(shù)，具體的計算過程可參考文獻[6]。

1.3 網(wǎng)格運動

沖刷過程在每個時間段上，可以看作是海床的垂向運動，因此網(wǎng)格變化集中于垂直方向，在來流平面上無變化。本文根據(jù)單元邊界處的輸沙率矢量計算單元格心處的高度變化，格點處的高度變化則通過周圍單元格心值的反距離加權得到：

式中：

dHn?格點高度變化（m）；

dHm?格心高度變化（m）；

wm=1/∥rmn∥?格心到格點距離的倒數(shù)。

得到床面邊界格點高度變化后，更新計算域底部邊界，通過底部邊界修改計算域內(nèi)的格點位置，輸入到流動計算模塊中開始下一個迭代步的計算。

1.4 沙滑模型

基礎周圍的床面遭受沖刷后形成沖刷坑，上游坑道的坡度往往大于下游坑道。對于無黏沙粒構成的海床，上游坑道坡度范圍一般為28°～32°，一般情況下當?shù)仄露炔粫^該范圍，因此也稱為臨界坡度。文獻[6]在實驗過程中發(fā)現(xiàn)，上游沖刷坑的部分區(qū)域內(nèi)，坑道坡度存在短時間內(nèi)超過臨界值的情況，隨后出現(xiàn)局部滑坡，且坑道坡度在重新穩(wěn)定后回歸正常值范圍。數(shù)值模擬中同樣需要考慮這一過程，因此采用韋雁機[11]等提出的方法，在完成全局沖刷計算后，搜索床面邊界內(nèi)存在局部坡度超過34°的單元，將單元節(jié)點沿水平軸翻轉，直至當?shù)仄露刃∮?0°，隨后重新形成完整的海床地貌。如圖1所示，計算過程中，單元節(jié)點的坐標變化遵循下式：

圖1 沙滑模型Fig.1 Sand slide model

式中：

? ??滑坡后的穩(wěn)定坡度（°），取 30°；

O?單元形心；

Ni滑坡前的單元格點；

?滑坡后的單元格點；

A?形心與格點連線向水平軸的垂線交點。

1.5 沖刷加速

沖刷問題是流動和海床形變互相作用的耦合問題，本文采用弱耦合分步交叉進行求解，在數(shù)值模擬過程中需要考慮流動迭代時間步和沖刷時間步長的適配性。海床形變相較于水流流動是個緩慢發(fā)展的過程，如果沖刷時間步長與流動步長相同，會導致迭代次數(shù)大幅增加，計算周期十分漫長，從而失去工程應用的價值；若沖刷步長過大，不僅導致沖刷和流動各自的計算結果喪失相關性，還會因為沖刷深度變化過大而造成計算的發(fā)散。文獻[6]在沖刷開始時限定沖刷計算步長，并隨著沖刷的進行逐步增大步長。Liang[12]等人通過數(shù)值模擬發(fā)現(xiàn)，要兼顧耦合計算的穩(wěn)定性和計算成本的經(jīng)濟性，需要每次迭代計算的深度變化不超過0.001D。本文采用文獻[12]所述方法，首先給定沖刷計算的初始時間步長，在計算得到全場深度變化后，將全場最大坑深變化其放大至限定值并反求對應的時間步長，以此縮短耦合計算的次數(shù)，節(jié)約計算資源。

2 算例驗證

2.1 算例設置

本文使用文獻[6]中的實驗作為算例，驗證計算模型的真實性與可靠性。實驗采用無黏沙粒，沙粒的50%中位數(shù)粒徑d50=2.6 mm，計算時不考慮粒徑分級，采用該數(shù)值作為代表。海床邊界層厚度δ=0.2 m，圓樁直徑D=0.1 m，近底平均流速Vref=0.46 m/s，邊界層內(nèi)速度廓線如圖2所示，表示為：

圖2 流速廓線Fig.2 Velocity profile

式中：

z??水線高度（m）；

ks?尼古拉斯平衡沙粒粗糙度（cm），??；

uf海床摩擦速度（cm/s），取4.8。

為保證算例中計算域入口處的湍流充分發(fā)展，需要先在相同粗糙度的無擾動水平床面上，以廓線指定的入流速度分布進行試算，得到出口處流速垂向分布及湍流耗散情況，并將出口數(shù)據(jù)作為沖刷算例的入口條件加入求解器。

流動計算求解前，需要在StarCCM+中設置計算參數(shù)和計算條件。計算模型采用隱式非定常，柱面邊界粗糙度設為2.5d50，壁面函數(shù)采用高Y+壁面處理，入口處的流速分布，湍動能和耗散率均從入口場函數(shù)中提取。迭代時間設為0.02 s，最大內(nèi)迭代步設為10。

2.2 網(wǎng)格信息

本文采用多塊分區(qū)六面體網(wǎng)格的劃分方式，計算域為沿流向18D，垂直流向15D，高度2D，域內(nèi)包含12個網(wǎng)格分區(qū)，每個分區(qū)內(nèi)的網(wǎng)格數(shù)量為323。樁和海床邊界附近網(wǎng)格加密，邊界附近的網(wǎng)格進行稀疏處理，網(wǎng)格具體形式如圖3所示。對于邊界條件的設置，入口條件設置為速度進口，出口條件指定為壓力出口；考慮邊界層頂部垂向運動對算例的影響較小，因此頂部及兩側邊界條件設為對稱邊界條件，柱面及底面邊界指定為無滑移壁面邊界條件。

圖3 網(wǎng)格劃分Fig.3 Detail of mesh for scour calculation Meshing

2.3 結果分析

初始時刻的床面切應力計算結果如圖4所示。通過與文獻[6]中計算結果進行對比可以看出，最大區(qū)別在于本文流動計算中，海床切應力核心區(qū)沿樁體后移，靠近樁體沿流向的最大厚度位置，同時核心區(qū)的后移導致樁后渦街的運動范圍受到擠壓，因此樁側剪切主導了樁周圍的切應力分布。而文獻[6]中的切應力核心區(qū)則較為靠前，接近樁前45°方向，且影響區(qū)域較為集中，與樁后渦街運動區(qū)域互不影響。這種區(qū)別的主要原因來自于軟件和文獻關于床面切應力計算方法和湍流相關計算項在公式選用上的差別。不同的海床切應力分布數(shù)據(jù)會導致沖刷計算結果產(chǎn)生一定差異，影響沖刷演化的進程，但由于沖刷演化是極為緩慢的過程，中間結果的微小擾動并不會對沖刷坑平衡形態(tài)產(chǎn)生太大的影響。

圖4 海床切應力云圖Fig.4 Seabed shear stress diagram

完成模型設置和網(wǎng)格劃分后，對模型進行了用時兩周的數(shù)值模擬，模擬結果如圖5所示。模型計算時間達到16 min時停止計算，此時的沖刷深度已達到0.81D。

結合圖5所示的沖刷過程和圖6所示的樁前坑深時變曲線，不難看出對于本算例中的樁基礎，30 s時上游沖刷坑已經(jīng)展現(xiàn)基本形態(tài)，且由于泥沙向下游堆積，樁兩側出現(xiàn)沙丘，沙丘的位置隨著沖刷的進行也在不斷后移，而沙丘背后伴隨產(chǎn)生的凹坑則是由于沙丘對水流的遮擋導致水流在沙丘上游積蓄能量，到達下游后加劇了對局部床面的沖擊。沖刷的前2 min內(nèi)，沖刷速度最快，但隨著坑深的增加和坑道范圍的擴展，沖刷速度從5 min開始逐漸變慢，這是因為樁前沖刷坑在沖刷初始階段尚未形成，此時樁前下洗渦和側向馬蹄渦的強度很大，對當?shù)卮裁娴臎_刷效應十分顯著。樁前沖刷坑逐漸形成后，坑深的增加和坑道坡度共同削弱了下洗渦和馬蹄渦對床面的掏空作用，使得沖刷速度放緩。但坑道整體由于沖刷的積累仍然處于擴展狀態(tài)，沖刷深度在16 min的計算時間內(nèi)達到0.81D，說明樁土結構遭受了十分顯著的破壞。

圖5 沖刷坑發(fā)展Fig.5 Scour pit development

如圖7所示的沖刷坑縱向?qū)ΨQ軸剖面可以看出，坑道下游在樁后脫落渦的影響下逐漸成形，但由于上游泥沙的持續(xù)補給，下游坑深的變化始終小于上游坑深，坑道向下游發(fā)展，坑深緩慢增加，但變化速度遠小于上游坑道。根據(jù)實驗結果[6]，沖刷平衡狀態(tài)下的樁下游坑深達到0.9D，稍低于上游1.2D的坑深。本文數(shù)值計算雖未達到平衡狀態(tài)，但根據(jù)下游坑深發(fā)展趨勢判斷，樁下游的沖刷計算速度偏小。這一問題也出現(xiàn)在文獻[13]中，原因是樁后沖刷由懸移質(zhì)沖刷主導，泥沙在渦流起動作用下離開床面，隨水流運動至下游。而對于推移質(zhì)沖刷，由于樁體阻擋了水流的沖擊，對樁后泥沙產(chǎn)生了遮擋效應。兩方面原因共同造成了推移質(zhì)模型下，樁后沖刷速度的計算結果偏小。

圖7 縱向?qū)ΨQ軸剖面Fig.7 Longitudinal symmetry axis profile

3 結論

本文通過結合CFD方法和泥沙輸運模型，以不可壓縮雷諾時均應力模型作為流動模塊，以自編程序為海床沖刷及網(wǎng)格運動模塊，建立了一套圓柱形單樁的實時沖刷動網(wǎng)格數(shù)值模擬流程。

1）通過本文模型與參考實驗進行對比，驗證了計算方法的可行性。

2）針對沖刷計算和流動計算在時間尺度上的不同，使用數(shù)值優(yōu)化算法加快了迭代速度，并保證了算法的魯棒性。

3）模擬過程中發(fā)現(xiàn)，樁基沖刷的最初階段是沖刷發(fā)展最快，預防效果最好的時段，實際施工過程中應當考慮改進基礎施工流程，在樁基入土后就開始沖刷的預防，從而更好地保護樁土結構的完整性。

4）后續(xù)工作將進行算法的改進，優(yōu)化流動計算和網(wǎng)格運動算法，改進沖刷計算模塊的程序結構，從而提高計算模型的整體精度，縮短預測時間。此外還將嘗試擴展程序功能，使之具備多樁基礎在浪、流荷載復雜流動條件下的數(shù)值模擬能力。