大尺度群樁隨機波浪力場的降維模擬

劉章軍，董心宇，張文遠

（武漢工程大學土木工程與建筑學院，湖北武漢 430074）

引 言

近年來，隨著中國海洋戰略的提出，大型離岸海工建筑物的荷載設計逐漸成為工程界的研究熱點。由于海工結構長期處于復雜海況中，因此，合理地模擬海工基礎性結構在服役期間所承受的波浪力是進行結構荷載設計的基礎和前提。其中，探討海工基礎性結構大小尺度樁柱波浪力模擬的異同對工程結構設計有著重要的參考意義。

在實際工程中，樁柱的大小尺度由直徑與波長之比決定，比值大于0.2 的樁柱被稱為大尺度樁柱，反之，則被稱為小尺度樁柱［1］。計算不同尺度樁柱的波浪力時所采用的原理和公式不盡相同。對于波浪力作用下的大尺度樁柱，國內外學者多采用繞射理論進行研究。在大尺度單樁方面，MacCamy 等［2］發展了以繞射理論為基礎的波浪力計算方法。邱大洪［3］基于一階近似橢圓余弦波理論研究了淺水區非線性波浪對圓柱樁的繞射作用，改進了作用于大尺度樁柱的波浪力計算公式。王俊杰等［4］應用波浪力計算公式，建立了大尺度樁柱波浪力的Nystr?m 數值模型。鄧莎莎等［5］以大尺度橋墩為研究對象，對該橋墩基礎所受的波浪載荷進行了數值模擬。在大尺度群樁方面，Spring 等［6］首先提出了針對群樁結構的矩陣算法。Linton 等［7］對該方法進行了改進，簡化了群樁波浪力的計算。陶建華等［8］基于矩陣算法，由譜分析法推導得到群樁系數。季新然等［9］通過物理模型實驗模擬了多向不規則波浪與群樁結構的相互作用。小尺度樁柱的波浪力計算則多以Morison 方程作為理論基礎［10］，該方程由Morison 于1950年提出。毛鴻飛等［11］基于Morison 方程，對大波幅波浪作用下，非完全淹沒水平圓柱的波浪力進行了數值模擬。程永舟等［12］基于波浪力的數值模擬，通過波浪水槽試驗，對孤立波作用下垂直樁柱周圍的局部沖刷特性進行了分析研究。

然而，上述研究僅對給定尺度的樁柱進行了受力分析，并未對樁柱尺度進行判別。實際上，在工程設計中通常依據設計波長進行大小尺度樁柱的判定，而設計波長則通常根據工程經驗來確定，這無疑會對工程應用造成困惑。此外，大小尺度樁柱的判別并非一成不變，其將隨著海況變化而發生改變。因此，如何選擇波浪力計算公式將非常困難。通常，隨機波浪力的模擬多采用傳統Monte Carlo 方法。然而，傳統Monte Carlo 方法所需的隨機變量與樣本數量巨大，導致結構動力計算量大，且生成的樣本無法在概率密度層面上精細地描述隨機過程的概率特性，這為結構精細化的動力反應分析與動力可靠度評價帶來挑戰。

針對上述研究現狀，本文首先提出了隨機波浪力作用下大小尺度樁柱的判別方法；其次，引入隨機函數降維思想［13?14］，實現了僅用一個基本隨機變量即可精確地模擬大尺度樁柱波浪力隨機場的目的；最后，考察了基于Morison 方程計算大尺度樁柱時所帶來的誤差，以及群樁效應對樁柱波浪力的影響。值得指出的是，由于本文方法生成的代表性樣本集合的概率信息完備，因此可與概率密度演化方法［15?16］相結合，進而實現海工結構在波浪力作用下的精細化動力反應分析與動力可靠性評價。

1 基于繞射理論的大尺度樁柱波浪力譜模型

在海工結構中，當樁柱之間的間距較小時，在波浪荷載的作用下，其受力與單樁有很大差別。樁中心距與樁直徑的比值不大于4 時，需要考慮周邊樁柱對波浪力的影響［17］。同時，大尺度樁柱對波動場的影響不可忽略，即，當波浪接觸到樁柱后，將在柱面產生一個向外擴散的波，在入射波與散射波的疊加達到穩態時，會形成一個新的波動場。

現假設水體無黏性且不可壓縮，波浪做有勢運動，群樁結構由N=p×q根樁柱構成，如圖1（a）所示，其中p代表行數，q代表列數。為簡便計，本文僅考慮直立圓柱樁的情況。對于任意的兩個樁柱i和j，以各自樁柱中心為原點建立局部的極坐標系，樁柱的相對位置及相關參數如圖1（b）所示。

圖1 群樁波浪力的計算簡圖Fig.1 Diagram for waves interact with multiple circular cyl?inders

假設樁柱i的中心坐標為(xi，yi)，樁柱j的中心坐標為(xj，yj)，根據繞射理論，波動場中任意一點的總速度勢空間分量為［7］：

其中：

式中 i 為虛數單位，g為重力加速度，H為波高，z是以海底為坐標原點的豎向坐標，ω為圓頻率，d為靜水面處水深，β為入射波方向；κ表示波數，由色散關系ω2=gκtanh(κd)求得；Ai，n為待定系數；Hn(κri)表示n階第一類漢克爾函數；Jn(κri)表示n階第一類貝塞爾函數；Yn(κri)表示n階第二類貝塞爾函數；J′n(κai)表示函數Jn(κai)對變量κai的一階導數；H′n(κai)表示函數Hn(κai)對變量κai的一階導數，ai為樁柱i的半徑。

當該點位于圓柱i的表面時，利用歐拉公式和Graf 加法定理［18］，式（1）可以改寫為：

式中m與n均為整數。

在式（5）中，通過對等號左邊的無窮項進行截斷，如取2M+1 項，則可得到N×(2M+1)個未知數Ai，n的聯立方程組，從而求解待定系數Ai，n。根據文獻［7］，本文中直接取M=6。

于是，圓柱j周圍速度勢的空間分量為：

在式（6）中，對貝塞爾函數應用Wronskian 關系式，并對通過速度勢求得的樁柱表面壓強進行積分和無因次化，得到波浪作用于圓柱j時波浪荷載的傳遞函數［7］：

式中ρ0為海水密度。

這樣，群樁結構中樁柱j任意坐標z處的波浪力功率譜密度函數為：

式中Sη(ω)為波高過程的單邊功率譜密度函數。

特別地，當結構為大尺度單樁（即N=1）時［2］，待定系數Aj，n=-in，此時傳遞函數為：

2 大尺度樁柱波浪力場的降維模擬

根據隨機過程模擬的降維方法，作用于大尺度群樁中第j個樁柱上的水平波浪力可表示為［19?21］：

式中 Δω為頻率離散步長；Nω為頻域離散點數；ωu為上限截止頻率，且有為一組標準的正交隨機變量集，滿足以下基本條件（必要條件）：

式中k，l=1，2，…，Nω；δkl為Kronecker 符號。

式中Θ為基本隨機變量，服從(0，2π]上的均勻分布；α為任意的確定性參數，在本文中取α=π/4。可以驗證，式（12）是滿足式（11）的基本條件（必要條件）。

在式（12）中，(=1，2，…，Nω) 與k(k=1，2，…，Nω)是一種確定性的一一映射關系，這種一一映射關系可由MATLAB 工具箱的內置函數，即rand（‘state’，0）和temp = randperm（Nω）來實現。(=1，2，…，Nω)與k(k=1，2，…，Nω)的一一映射關系即為kˉ=temp（k），這恰是降維模擬方法的一個充分條件。

于是，將式（12）代入式（10），即可得到隨機波浪力場的降維模擬公式：

式中?kˉ(Θ)為基本隨機變量Θ的函數，其表達式為：

通過選取基本隨機變量Θ的代表性點集，即可得到大尺度樁柱隨機波浪力的代表性時程集合，進而結合概率密度演化理論，進行樁柱結構隨機波浪力作用下的動力反應分析與可靠度精細化計算。

3 樁柱尺度的判別方法

繞射理論僅適用于直徑與波長之比大于0.2 的大尺度樁柱波浪力計算，但對于隨機波浪而言，其波長是一個變量，這將導致樁柱尺度的判別較為困難。因此，可以通過確定其特征頻率的方式來確定其特征波長。這樣，便可通過特征波長來判別其是否屬于大尺度樁柱。

對于窄譜波浪，其能量主要集中于波浪譜中心位置處。因此，本文采用波浪譜重心處的頻率作為特征頻率ˉ，其計算如下：

式中m0和m1分別為0 階與1 階譜矩，即：

根據特征頻率，利用色散關系，即可得到與特征頻率對應的特征波長為：

這是一個超越方程，可以用MATLAB 中fsolve函數進行數值求解。

這樣，可以給出判別大尺度樁柱的量化公式如下：

為了進一步說明式（18）的合理性，本文利用式（18）對文獻［22?23］中樁柱的尺度進行判別，所得結果與相應的文獻一致。此外，對于同一樁柱，其尺度并非是一成不變的，也會隨著海況的變化而發生改變。圖2給出了不同直徑樁柱隨風速v19.5的尺度變化情況。 其中，波高譜采用Pierson?Moskow?itz 譜［24］：

圖2 樁柱尺度隨風速變化情況Fig.2 Cylinder scale varying with wind speed

式中 無因次常數α=8.1×10-3，β=0.74，v19.5為海面上19.5 m 高度處的平均風速。

由圖2可知，樁柱直徑越大，使樁柱尺度發生改變所需的臨界風速也越大，即風速的變化會影響樁柱尺度的判別。因此，在應用中，首先需要根據樁柱直徑和風速對樁柱尺度進行判別，在此基礎上，利用繞射理論或Morison 方程進行隨機波浪力的降維模擬。

4 工程算例

4.1 數值模擬分析

為了說明本文建議的降維模擬方法的有效性，首先考慮群樁的特例，即大尺度單樁的隨機波浪力場模擬，其參數取值如表1所示。

表1 模擬參數的取值Tab.1 Calculation parameters and corresponding values

利用本文所提降維模擬方法生成610 條隨機波浪力代表性樣本集合。圖3給出了z=20 m 處隨機波浪力的第400 條和第500 條代表性樣本，它們具有隨機波浪力的典型特征。

圖3 隨機波浪力代表性樣本Fig.3 Representative samples of random wave force

在模擬精度方面，圖4為z=12 m 處隨機波浪力的自相關函數比較，圖5為z=12 m 和z=15 m之間的互相關函數比較。可見，采用本文方法模擬的隨機波浪力自（互）相關函數擬合效果要好于Monte Carlo 方法，初步驗證了降維模擬方法的有效性。圖6，7 為模擬的均值，標準差與目標值比較（以z=12 m 為例），圖8為z=12 m 處隨機波浪力的功率譜密度函數在有效頻率0.9～2.5 rad/s 范圍內與目標值的比較。其中，有效頻率指提供頻譜總能量99.9%的頻譜范圍。表2為隨機波浪力樣本均值誤差，標準差誤差及功率譜誤差的具體結果，從表2可以看出，本文方法的均值、標準差及功率譜密度函數均與目標值擬合一致，進一步驗證了本文方法的有效性。且其均值誤差、標準差誤差及功率譜誤差均小于Monte Carlo 方法。在模擬效率方面，采用本文方法生成單條樣本耗時6.467 s，Monte Carlo 方法為6.623 s。可見，本文方法在模擬效率方面與Monte Carlo 方法大致相當，但本文方法的模擬精度更高，體現了其優越性。

圖4 自相關函數比較Fig.4 Comparison of auto?correlation functions

圖5 互相關函數比較Fig.5 Comparison of cross?correlation functions

圖6 均值的比較Fig.6 Comparison between simulated mean values and target value

圖7 標準差的比較Fig.7 Comparison between simulated standard deviations and target value

圖8 功率譜的比較Fig.8 Comparison between simulated power spectrum and target value

表2 模擬精度比較Tab.2 Comparison of simulation accuracy

4.2 大尺度樁柱計算公式對比分析

為了考察大尺度樁柱基于Morison 方程帶來的誤差。現以4.1 中樁柱模型和波浪力參數為例，分別基于繞射理論和Morison 方程對大尺度樁柱進行計算，并將兩者的結果進行對比，分析基于Morison方程計算大尺度樁柱所帶來的誤差。圖9為z=18 m 處基于繞射理論和Morison 方程所生成的第400 條隨機波浪力代表性樣本比較。

圖9 z=18 m 處隨機波浪力代表性樣本的比較Fig.9 Comparison of representative sample of random wave force at z=18 m

從圖9可知，在樁柱z=18 m 處，基于Morison方程計算得到的樁柱波浪力明顯偏大。為進一步分析兩種計算公式的差異性，通過對610 條隨機波浪力代表性樣本的統計分析，圖10 給出了樁柱不同位置處基于兩種計算公式得到的最大波浪力均值對比。

從圖10 可以看出，在大尺度樁柱豎直方向上存在某一臨界點，在該臨界點以下，基于繞射理論模擬的最大波浪力均值大于基于Morison 方程模擬的結果，在該臨界點以上，則小于基于Morison 方程模擬的結果。

圖10 樁柱不同位置處最大波浪力均值的比較Fig.10 Comparison of maximum wave force at different po?sitions of the cylinder

圖11 為兩種計算公式得到最大波浪力均值的相對誤差隨海況（即特征波長Lˉ）的變化情況，其中，相對誤差ε計算如下：

式中Smax為基于Morison 方程得到的最大波浪力均值，Lmax為基于繞射理論得到的最大波浪力均值。

由圖11 可知，當Smax小于Lmax時，隨著z的增大，相對誤差ε的絕對值將隨之減小，當Smax大于Lmax后，隨著z的增大，相對誤差ε也將隨之增大。同時，隨著特征波長Lˉ的增大，相對誤差ε=0 時所對應的z值也將隨之增大。

圖11 在3 種海況下樁柱不同位置的差值比Fig.11 The difference ratio of cylinder at different positions under three sea conditions

4.3 群樁效應影響分析

為了說明大尺度群樁效應對隨機波浪力場模擬的影響，現基于繞射理論對平行于波向的串列三樁（N=3×1）和其中某一單樁（N=1）分別進行隨機波浪力場的降維模擬。群樁的排列及波浪的作用方向如圖12 所示，其中樁柱模型和波浪參數均與上文相同，樁的中心距為8 m。

圖12 串列三樁位置示意圖Fig.12 Schematic diagram of three cylinder positions

圖13 給出了z=12 m 處1 號樁按單樁和群樁情況下第400 條隨機波浪力代表性樣本的比較。從圖中可見，1 號樁在群樁情況下的波浪力最大值大于單樁情況。事實上，通過對610 條代表性樣本的統計分析，可得1 號樁在群樁情況下的最大波浪力均值也大于單樁情況。因此，對于1 號樁，群樁效應將會使波浪力在總體上偏大。

圖13 1 號樁在單樁情況與群樁情況下隨機波浪力代表性樣本比較Fig.13 Comparison of representative samples of random wave force of 1st cylinder in single pile case and pile group case

同樣地，對2 號和3 號樁進行類似的群樁影響分析，可知群樁效應會使2 號樁波浪力在總體上略偏大，而使3 號樁波浪力在總體上偏小。

圖14 給出了在群樁情況下，樁柱z=12 m 處1 號樁與3 號樁第400 條隨機波浪力代表性樣本比較，顯然1 號樁受到的波浪力明顯大于3 號樁，且通過對610 條代表性樣本的統計分析，可知1 號樁受到的最大波浪力均值比3 號樁大10.35%，該結論與文獻［8］譜分析法得到的群樁系數結果相似，進一步驗證了本文方法的正確性。

圖14 1 號樁與3 號樁隨機波浪力代表性樣本比較Fig.14 Comparison of representative samples of wave force between 1st cylinder and 3rd cylinder

5 結 論

本文首先根據波高過程的功率譜密度函數，建議了不規則波浪作用下樁柱大小尺度的判別方法。然后，基于繞射理論推導了大尺度樁柱隨機波浪力的功率譜密度函數，并引入隨機函數的降維思想，實現了隨機波浪力場的降維模擬。最后，探討了大尺度樁柱基于Morison 方程所帶來的誤差，以及群樁效應對樁柱波浪力的影響。通過數值分析，可得以下結論：

（1）本文方法實現了僅用一個基本隨機變量即可精細地模擬大尺度群樁隨機波浪力場，極大地減少了隨機變量的數量。而且，相對于傳統Monte Carlo 方法，本文方法具有更高的精度，并能夠與概率密度演化方法結合來實現海工結構的動力反應與可靠度精細化分析，因此具有明顯的優勢。

（2）大尺度樁柱在豎直方向上存在某一臨界點，在該臨界點以下，采用Morison 方程計算樁柱時會導致波浪力偏小，在該臨界點以上，則會導致波浪力偏大，且該臨界點隨著特征波長Lˉ的增加而向樁頂移動。

（3）當單向不規則波浪作用于群樁結構時，沿著波浪傳播方向，群樁效應將使第1 樁柱的波浪力在總體上比單樁情況偏大。后面的樁柱由于受到前面樁柱的遮蔽，其波浪力逐漸減小。