劉章軍，董心宇，張文遠(yuǎn)

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

引 言

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

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

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

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

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

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

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

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

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

其中：

式中 i 為虛數(shù)單位，g為重力加速度，H為波高，z是以海底為坐標(biāo)原點(diǎn)的豎向坐標(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)該點(diǎn)位于圓柱i的表面時(shí)，利用歐拉公式和Graf 加法定理［18］，式（1）可以改寫為：

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

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

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

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

式中ρ0為海水密度。

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

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

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

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

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

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

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

式中Θ為基本隨機(jī)變量，服從(0，2π]上的均勻分布；α為任意的確定性參數(shù)，在本文中取α=π/4?？梢则?yàn)證，式（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ω）來實(shí)現(xiàn)。(=1，2，…，Nω)與k(k=1，2，…，Nω)的一一映射關(guān)系即為kˉ=temp（k），這恰是降維模擬方法的一個(gè)充分條件。

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

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

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

3 樁柱尺度的判別方法

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

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

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

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

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

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

為了進(jìn)一步說明式（18）的合理性，本文利用式（18）對文獻(xiàn)［22?23］中樁柱的尺度進(jìn)行判別，所得結(jié)果與相應(yīng)的文獻(xiàn)一致。此外，對于同一樁柱，其尺度并非是一成不變的，也會(huì)隨著海況的變化而發(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)速的變化會(huì)影響樁柱尺度的判別。因此，在應(yīng)用中，首先需要根據(jù)樁柱直徑和風(fēng)速對樁柱尺度進(jìn)行判別，在此基礎(chǔ)上，利用繞射理論或Morison 方程進(jìn)行隨機(jī)波浪力的降維模擬。

4 工程算例

4.1 數(shù)值模擬分析

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

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

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

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

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

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

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

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

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

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

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

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

由圖11 可知，當(dāng)Smax小于Lmax時(shí)，隨著z的增大，相對誤差ε的絕對值將隨之減小，當(dāng)Smax大于Lmax后，隨著z的增大，相對誤差ε也將隨之增大。同時(shí)，隨著特征波長Lˉ的增大，相對誤差ε=0 時(shí)所對應(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)對隨機(jī)波浪力場模擬的影響，現(xiàn)基于繞射理論對平行于波向的串列三樁（N=3×1）和其中某一單樁（N=1）分別進(jìn)行隨機(jī)波浪力場的降維模擬。群樁的排列及波浪的作用方向如圖12 所示，其中樁柱模型和波浪參數(shù)均與上文相同，樁的中心距為8 m。

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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