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

劉章軍，董心宇，張文遠

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

引 言

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

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

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

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

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

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

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

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

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

其中：

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

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

式中m與n均為整數(shù)。

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

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

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

式中ρ0為海水密度。

這樣，群樁結(jié)構(gòu)中樁柱j任意坐標(biāo)z處的波浪力功率譜密度函數(shù)為：

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

特別地，當(dāng)結(jié)構(gòu)為大尺度單樁（即N=1）時［2］，待定系數(shù)Aj，n=-in，此時傳遞函數(shù)為：

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

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

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

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

式中Θ為基本隨機變量，服從(0，2π]上的均勻分布；α為任意的確定性參數(shù)，在本文中取α=π/4?？梢则炞C，式（12）是滿足式（11）的基本條件（必要條件）。

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

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

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

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

3 樁柱尺度的判別方法

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

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

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

根據(jù)特征頻率，利用色散關(guān)系，即可得到與特征頻率對應(yīng)的特征波長為：

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

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

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

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

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

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

4 工程算例

4.1 數(shù)值模擬分析

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

表1 模擬參數(shù)的取值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 處隨機波浪力的自相關(guān)函數(shù)比較，圖5為z=12 m 和z=15 m之間的互相關(guān)函數(shù)比較?？梢姡捎帽疚姆椒M的隨機波浪力自（互）相關(guān)函數(shù)擬合效果要好于Monte Carlo 方法，初步驗證了降維模擬方法的有效性。圖6，7 為模擬的均值，標(biāo)準(zhǔn)差與目標(biāo)值比較（以z=12 m 為例），圖8為z=12 m 處隨機波浪力的功率譜密度函數(shù)在有效頻率0.9～2.5 rad/s 范圍內(nèi)與目標(biāo)值的比較。其中，有效頻率指提供頻譜總能量99.9%的頻譜范圍。表2為隨機波浪力樣本均值誤差，標(biāo)準(zhǔn)差誤差及功率譜誤差的具體結(jié)果，從表2可以看出，本文方法的均值、標(biāo)準(zhǔn)差及功率譜密度函數(shù)均與目標(biāo)值擬合一致，進一步驗證了本文方法的有效性。且其均值誤差、標(biāo)準(zhǔn)差誤差及功率譜誤差均小于Monte Carlo 方法。在模擬效率方面，采用本文方法生成單條樣本耗時6.467 s，Monte Carlo 方法為6.623 s。可見，本文方法在模擬效率方面與Monte Carlo 方法大致相當(dāng)，但本文方法的模擬精度更高，體現(xiàn)了其優(yōu)越性。

圖4 自相關(guān)函數(shù)比較Fig.4 Comparison of auto?correlation functions

圖5 互相關(guān)函數(shù)比較Fig.5 Comparison of cross?correlation functions

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

圖7 標(biāo)準(zhǔn)差的比較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 方程帶來的誤差?，F(xiàn)以4.1 中樁柱模型和波浪力參數(shù)為例，分別基于繞射理論和Morison 方程對大尺度樁柱進行計算，并將兩者的結(jié)果進行對比，分析基于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 條隨機波浪力代表性樣本的統(tǒng)計分析，圖10 給出了樁柱不同位置處基于兩種計算公式得到的最大波浪力均值對比。

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

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

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

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

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

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

4.3 群樁效應(yīng)影響分析

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

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

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

圖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 號樁進行類似的群樁影響分析，可知群樁效應(yīng)會使2 號樁波浪力在總體上略偏大，而使3 號樁波浪力在總體上偏小。

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

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

5 結(jié) 論

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

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

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

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