橢圓形柱體地震動水壓力的簡化分析方法1

王丕光 黃義銘 趙 密 杜修力 張麗華

1）清華大學土木工程系，北京 100084

2）北京工業大學城市與工程安全減災教育部重點實驗室，北京 100124

3）華北科技學院，北京 101601

引言

近年來，隨著世界經濟的發展，中國近海結構的研究得到了快速發展，如跨海橋梁、海上風電、人工島和石油平臺等。然而，中國福建、廣東沿海和中國臺灣等地區的近海結構面臨著地震作用的威脅。在地震作用下，橋梁與周圍水體的相互作用會對橋墩產生動水壓力。地震作用下水體與結構的動力相互作用會對結構產生附加動水壓力，其不僅會改變結構的動力特性，還會對水中結構的動力反應產生較大影響（Liaw等，1974；Han等，1996；黃信等，2011a；江輝等，2014）。因此，充分認識和研究地震激勵下近海結構的動水壓力，對于近海結構的抗震設計具有重要意義。

在實際近海工程中存在各種截面形式的墩柱結構，如圓形和橢圓形等。目前，國內外學者對圓柱結構地震動水壓力的研究已取得大量的研究成果。Liaw等（1974）基于輻射波浪理論推導了圓柱結構的動水壓力公式，結果表明，對于細長結構可以忽略水體壓縮性的影響；表面波僅在低頻時對動水壓力的影響比較明顯。忽略水體壓縮性和表面波時，該動水壓力可視為部分水體質量與結構加速度的乘積，這部分水體稱為“附加質量”。Williams（1986）采用邊界積分方法研究了地震作用下水中圓柱結構的動力反應。Tanaka等（1988）分析了水平地震作用下彈性圓柱體動水壓力的附加質量系數和阻尼系數。Han等（1996）提出了計算水中圓柱結構自振頻率的簡化公式。賴偉等（2004）提出了一種圓形橋墩上地震動水壓力的半解析半數值方法。黃信等（2011b，2012）討論了水體壓縮性、表面波和水底吸收邊界對圓形橋墩地震動水壓力的影響。杜修力等（2012）、Du等（2014）提出了可壓縮水體條件下圓柱結構地震動水壓力的時域算法和時域簡化公式。Goyal等（1989）、Li等（2013）和Jiang等（2017）提出了圓柱形橋墩上地震動水壓力的附加質量簡化計算公式。另外，Liao（1985）研究了水中多個圓柱結構的動力相互作用；Sun等（1991）研究了軸對稱結構的地震動水壓力；Avilés等（2001）討論了海底吸收條件對軸對稱結構地震動水壓力的影響；Wei等（2015）提出了水中軸對稱結構的地震設計和分析的簡化方法。

上述分析表明，地震作用下圓柱結構與水體的動力相互作用問題已取得了可供實際應用的研究成果，但對橢圓柱體結構動水壓力方面仍鮮有研究。本文基于橢圓坐標系，根據線性輻射波浪理論，利用分離變量法推導了橢圓柱體結構動水壓力的解析解，并采用有限元方法建立了地震作用下水體與橢圓柱體結構動力相互作用的簡化分析方法。

1 橢圓形柱體動水壓力求解

地震作用下橢圓柱體結構與水體動力相互作用的分析模型如圖1所示，a和b分別表示橢圓外徑的半長軸和半短軸，a1和b1分別表示橢圓內徑的半長軸和半短軸，h為水深，H為柱體高度。直角坐標系下，z軸沿柱體軸線向上，坐標原點位于柱體底部；地基為剛性，地面運動加速度為。水體假定為無旋、無粘、可壓縮的小擾動流體，并忽略表面重力波的影響。

圖1 橢圓柱體與水體相互作用分析模型 Fig.1 Analytical model of the interactiin of elliptic cylinder and water

圖2 橢圓坐標系 Fig.2 Elliptical cylindrical coordinate

1.1 控制方程和邊界條件

由于柱體截面呈橢圓形，故在橢圓坐標系下求解橢圓柱體的地震動水壓力。橢圓坐標系如圖2所示。直角坐標系與橢圓坐標系直接的變化關系為：

在橢圓柱坐標系下，以動水壓力p表示的流體控制方程為（Bhatta，2005）：

水體底部、靜水表面和無窮遠邊界條件為：

地震動沿長軸方向時，η=0和η=0.5π面的對稱邊界條件為：

地震動沿短軸方向時，η=0和η=0.5π面的對稱邊界條件為：

地震動沿長軸和短軸方向時，水體與結構交界面邊界條件分別為：

式中，ρ表示水體密度，表示結構的加速度。

1.2 分離變量求解

在橢圓柱坐標系下，p(ξ,η,z)可分離變量為：

將式（12）代入式（4）整理得到3個解耦的方程：

式中，λj和a0是分離變量常數，q=μ2λ2j/4為一無量綱參數。

由方程（13）和邊界條件式（5）、（6）可得：

式中，dj為待定系數，λj= (2j- 1)π /2h，j=1，2，…。通過正交歸一化可得

方程（14）為修正的角向馬蒂厄方程（熊天信，2014），其解為第一類角向馬蒂厄函數cen(η,-q)和sen(η,-q)，即：

函數cen(η,q)和sen(η,q)的歸一化正交關系為：

方程（15）為修正的徑向馬蒂厄方程（熊天信，2014）。整數階徑向馬蒂厄方程的完全解為：

式中，是任意常數；函數Ien(ξ,-q)和Ion(ξ,-q)稱為第一類變形貝塞爾型徑向馬蒂厄函數，是單調遞增函數；函數Ken(ξ,-q)和Kon(ξ,-q)稱為第二類變形貝塞爾型徑 向馬蒂厄函數，是單調遞減函數。

1.3 柱體的動水壓力

當地震動沿長軸方向時，根據邊界條件式（5）—（8）和（10），并利用式（18）的正交性可得式（12）的解為：

當地震動沿短軸方向時，根據邊界條件式（5）—（7）、（9）和（11），并利用式（18）的正交性可得式（12）的解為：

式中，Ko1′ (ξ,-q)為函數Ko1(ξ,-q)的一階導數。

地震作用沿長軸方向時，橢圓截面柱體表面單位高度上的動水力為：

將式（21）代入式（23）整理得：

地震作用沿短軸方向時，橢圓截面柱體表面單位高度上的動水力為：

將式（22）代入式（25）整理得：

1.4 方法驗證

首先，通過圓形橋墩動水力的解析解（Li等，2013）驗證本文提出的橢圓形截面柱體動水壓力的解析解。假定結構為剛性，將單位高度動水力沿高度積分，可得到橢圓形柱體的均布附加質量系數為：

式中，m0為單位高度水體的附加質量，沿長軸 方向時m0=ρπb2，沿短軸方向時m0=ρπa2。 圖3為本文圓柱體附加質量系數的解與圓柱解析解的對比。由圖3可以看出，本文解與圓柱解析解很好地吻合。

進一步通過Wang等（2019）提出的數值方法驗證本文提出的橢圓形截面柱體動水壓力的解析解，該方法實質上是用垂向特征函數展開，將三維輻射問題簡化為二維問題；然后用外域特征函數展開解與結構截面附近內域有限元聯合求解。圖4為本文橢圓柱體附加質量系數的解與數值解的對比，橢圓柱的長軸和短軸尺寸分別為a=20m、b=10m。由圖4可以看出，本文的解與橢圓柱數值解很好地吻合。

圖3 圓柱解析解與本文解的對比 Fig.3 The proposed solution compared with the analytical solution for circular cylinders and this paper

圖4 橢圓柱數值解與本文解的對比 Fig.4 The proposed solution compared with the numerical solution for elliptical cylinders and this paper

2 水中橢圓柱體結構的地震反應求解

采用有限元方法將柔性柱體結構離散為梁單元，則地震作用下水中橢圓柱體的動力方程為（王勖成，2003）：

式中，M、C和K分別為結構質量、阻尼和剛度矩陣；us為結構柔性位移列向量；F為地震 動水力列向量。將地震動水力式（24）或（26）進行有限元離散，則F可以表示為

式中，Mp為附加質量矩陣，N為形函數列向量；Sj如式（24b）或（26b）。

將式（29）代入式（28）整理得：

式（32）可通過數值積分方法Newmark-β方法求解（劉晶波等，2005）。需要指出的是，附加質量矩陣Mp是滿陣的，難以在商業有限元中實現。本文將附加質量矩陣Mp的每一行元素進行集中化形成1個集中的附加質量矩陣Mg，即：

式中，mi,j=mj,i，mi=mi,1+mi,2+ …… +mi,L-1+mi,L，L表示柱體結構水下的節點數目。

通過柔性橢圓柱體驗證提出的集中附加質量矩陣的精度。圖5為地面運動位移，脈沖持時0.2s，頻譜覆蓋了地震作用頻段。橢圓柱體尺寸為H=80m、a=40m、b=20m、a1=30m、b1=15m；密度和彈性模量分別為2500kg/m3和30GPa；水深h=80m；梁單元長8m；不考慮阻尼作用。圖6為地面運動沿長軸方向時，采用精確附加質量模型和集中附加質量模型計算得到的結構頂部位移時程。由圖可見，采用集中附加質量方法計算的位移反應的周期偏大，即結構柔性運動引起的水體附加質量偏大。

圖5 地面位移時程和傅里葉譜 Fig.5 The time history of displacement of the ground motion and its Fourier spectrum

圖6 參考解和集中附加質量模型位移時程的對比 Fig.6 The comparison of time history of displacement of the cylinder obtained by reference solution and lumped added mass method

因此，本文進一步提出修正的集中附加質量方法，即將式（32）修正為如下形式：

式中，α≤1為集中附加質量矩陣修正系數。圖7為集中附加質量矩陣修正系數隨無量綱參數寬深比（l）和長短軸比（δ）的變化；無量綱參數寬深比（l）和長短軸比（δ）的定義如公式（36）和（37）所示：

式中，地震方向沿x軸時D= 2a，地震方向沿y軸時D= 2b。

圖7 集中附加質量修正系數 Fig.7 Correction factors for the lumped added mass matrix

圖8為地面運動沿長軸方向時，采用滿陣附加質量模型和修正集中附加質量模型計算得到的結構頂部位移時程，滿陣附加質量模型解為參考解。由圖可見，采用修正集中附加質量方法計算的位移反應與參考解很好地吻合。

圖8 參考解和修正集中附加質量模型的 位移時程比較 Fig.8 The comparison of time history of displacement of the cylinder obtained by reference solution and modified-lumped added mass method

3 橢圓形柱體均布附加質量簡化公式

式（34）所示的集中附加質量矩陣與通過剛性結構動水力解析公式計算得到的附加質量是一致的。需要指出的是，動水力附加質量解析公式的數學表達復雜，難以在工程中進行應用。因此，提出附加質量的簡化計算公式。假定附加質量沿高度均勻分布，相應的附加質量系數定義如式（27）所示。將式（24a）或（26a）代入式（27）整理得：

式中，Sj的定義見式（24b）或（26b）。需要指出的是，簡化公式是在0.2 ≤l≤ 2和0.2 ≤δ≤5的范圍內擬合得到的。

通過曲線擬合，圓形柱體均布附加質量系數的簡化公式為：

當地震動沿長軸方向時，通過曲線擬合得到橢圓形柱體沿長軸方向的均布附加質量系數的簡化公式為：

當地震動沿短軸方向時，通過曲線擬合，得到橢圓形柱體沿短軸方向的均布附加質量系數的簡化公式為：

圖9為橢圓柱體附加質量系數的解析解與簡化公式的對比，圖10則為簡化公式的誤差，可以看出簡化公式與解析解吻合較好。

圖9 附加質量系數解析解和簡化公式的對比 Fig.9 The comparison of coefficient of the added mass obtained by analytical solution and simplified formula

圖10 附加質量系數簡化公式的誤差 Fig.10 The error of the simplified formula for the added mass coefficient

4 結語

（1）基于橢圓坐標系，采用分離變量法將拉普拉斯方程轉換為馬蒂厄方程。通過求解馬蒂厄方程，提出了橢圓柱體結構地震動水壓力的解析解。

（2）建立了地震作用下橢圓柱體結構與水體相互作用的動力有限元方程，結構的動水力通過附加質量矩陣施加，該矩陣是滿陣的。

（3）為便于橢圓柱體結構附加質量矩陣在商業有限元中實現，提出了集中附加質量矩陣的方法，該方法中柔性運動引起的附加質量為集中附加質量矩陣和修正系數的乘積。

（4）基于剛性柱體結構動水力的解析解，通過曲線擬合的方法建立了橢圓柱體結構動水力的均布附加質量簡化公式，公式中的系數僅與無量綱參數寬深比和長短軸比相關。