武海浪，陳徐均，黃亞新，沈海鵬，苗玉基

（解放軍理工大學 野戰(zhàn)工程學院，南京210007）

0 引 言

帶支腿浮式結構是指在浮式平臺的基礎上設計加裝一定數(shù)量的可升降支腿而形成的海洋工程結構。結構在提樁和放樁過程中，樁腿并不觸及海底，主體平臺處于帶樁腿漂浮的狀態(tài)。對于一個確定的浮式結構，按照船舶水彈性理論對其進行水動力分析時，一般需要進行濕面網(wǎng)格定義、格林函數(shù)分析和水動力響應計算幾個步驟[1-5]。要研究水下樁腿長度變化以后的浮式結構的水彈性性質，就要針對變化的參數(shù)，重新進行浮式結構的干模態(tài)參數(shù)計算、重新定義新的濕面網(wǎng)格，然后進行水動力計算。文獻[6-8]就是采用了根據(jù)水下樁腿長度分工況的方法進行分別計算，認為水下樁腿的長度對浮式結構的水動力系數(shù)、波浪激勵力以及水動力響應都有顯著的影響。

三維線性水彈性分析的核心是格林函數(shù)分析，是比較消耗CPU時間的，尤其是當參與計算的面元較多的時候，計算效率會明顯地下降。前人對三維水彈性分析算法效率給予了關注，致力于提高計算效率：戴愚志等[9-10]用GMRES（Generalized Minimal RESidual algorithm）方法求解了大型離岸結構水彈性分析所獲得的復系數(shù)線性方程組，并按照向后誤差分析方法給出了它的終止準則。數(shù)值試驗表明GMRES方法可提高效率，減少計算時間，優(yōu)于直接方法。Wang等[11]引入兩項技術到三維水彈性計算當中，以提高超大型浮式結構水彈性計算的效率。其一是使用新的截斷標準來降低格林函數(shù)及其導數(shù)的分析時間，其二是使用了交互式稀疏求解方法代替高斯方法求解水彈性方程，結果顯示有效地提高了計算效率。Wu[12]為了提高水彈性程序的計算效率在其程序中使用了對稱復合勢方法，這種分析方法也一直沿用至今，所以在后面的分析中，本文也只建立了一半的濕面單元。

浮式結構提樁放樁時只是部分參數(shù)在確定的狀態(tài)下進行了微調，并未發(fā)生本質的變化。前面的研究方案盡管直觀、易行，但這種參數(shù)微調就不得不重新建模、反復分析和重新計算，這似乎并不夠高效，每一次新的計算包含了大量的重復工作。若采用這樣的計算方案，在這類浮式結構的設計工作中，或許將降低工作效率。本文在研究者所熟知的三維線性水彈性理論體系下，利用格林函數(shù)矩陣的性質，針對帶支腿浮式結構變參數(shù)的水動力系數(shù)計算問題提出了一種相對較高效的分析方法，主要是在格林函數(shù)分析和線性方程組求解方面作了相關的研究。

1 算法原理

1.1 變物面條件問題

假設彈性浮體周圍為均勻不可壓縮、無粘的理想流體，流場的運動是無旋的，自由表面波為微幅且?guī)е雀∈浇Y構在波浪中平衡位置附近作微幅運動，則流場的運動可以在歐拉坐標系中采用三維勢流理論來描述。

根據(jù)對流場的基本假定，在空間固定坐標系中帶支腿浮式結構周圍，如圖1所示，邊界SL+SF+S∞+SB所包圍的流場內輻射勢定解問題可以描述為：

其中：▽2為拉普拉斯算子，φr表示輻射波速度勢的空間分量，當浮體在平衡位置附近以入射波圓頻率ω，以第r階干模態(tài)振型作單位主坐標幅值振蕩時所誘導的流場周圍的流體運動速度勢。SL為物面邊界條件，表示帶支腿浮式結構的濕面積，主要由浮式結構平臺的底部、側舷和樁腿三個部分組成。由于樁腿是可以上下移動的，物面條件是可變的，物面條件的改變將導致方程解的改變。SF表示自由面邊界條件，S∞為無窮遠邊界條件，SB為海底邊界條件。

1.2 計算方法

假定：（1）帶支腿浮式結構中樁腿是平直的柱狀結構；（2）忽略因為樁腿的上升或者下放而導致的主體平臺水線面的變化；（3）樁腿底板的影響很小。選取帶支承腿浮式結構水下樁腿最長L1時為初始狀態(tài)，即將樁腿盡可能向下放置但不觸及海底的狀態(tài)，物面條件記為SL1。采用格林函數(shù)法對結構的初始狀態(tài)進行水動力分析。

對初始狀態(tài)進行水動力建模。圖2為任意一帶支腿浮式結構，此時樁腿處于水下伸長狀態(tài)，濕面網(wǎng)格分為船身主體部分和樁腿部分。整個結構一共劃分的單元數(shù)記為N，主體平臺單元編號記為1～Nhull，后面Nhull+1～N個編號給樁腿。對樁腿濕面元編號原則參照表1，樁腿濕面元編號的方法總體上遵循了以水深為標尺，編號大小遵循“由上到下”的原則。即離水線面向下越遠的面元，編號越靠后，豎直方向上，處于同一高度的面元編號越為接近。

表1 樁腿狀態(tài)編號及濕面元編號Tab.1 The numbering schedule of the wetted panels

樁腿處于初始狀態(tài)時，經(jīng)開拓后場內任一點P的輻射勢可按源分布表示為：

其中：G（ P， ）Q 為格林函數(shù)，應由初始物面條件SL1決定。使用文獻[12]給出的有限水深無航速問題的三維脈動源格林函數(shù)，其級數(shù)表達形式為：

積分方程離散為一組線性代數(shù)方程組，表示為：

上式可寫成矩陣形式，有：

其中：

當樁腿同時緩慢、均勻地由初始狀態(tài)長度L1恰好抬升到Li時，從實際效果上看是樁腿的水下長度變短。若重新進行水動力建模以及新的格林函數(shù)分析，可得到新的格林函數(shù)矩陣Gnew和新的影響矩陣Cnew。

其中：g′和c′分別代表樁腿底部的面元的影響。基于假定（3），可以認為樁腿底板的作用很小，忽略矩陣中g′和 c′在矩陣中的影響，即：

這樣，當樁腿由初始狀態(tài)連續(xù)地作抬升動作時，可以視為格林函數(shù)矩陣規(guī)模不斷縮小的過程，即連續(xù)地有矩陣中相應的行列退出計算。如圖4（b）所示，根據(jù)編號法則，位于樁腿底部面元的行列編號靠后，當樁腿由初始狀態(tài)作抬升動作時，等效于樁腿底部面元退出計算，相應的就是將格林函數(shù)矩陣中靠右的部分一層一層地剝離矩陣。直到樁腿全部抬升起來，等同于編號為Nhull+1～Nn-1的網(wǎng)格退出計算。這樣有規(guī)律的按編號進行矩陣縮減可以得到一系列格林函數(shù)矩陣，實際上是給出了流體域內速度勢定解問題中物面由SL1逐步變化到SLn的方程組的近似格林函數(shù)解。可用公式表示為：

格林函數(shù)及其影響矩陣的變化過程表示為

在濕面上進行源強積分計算，可得到輻射勢：

以上便為帶支腿浮式結構水動力系數(shù)快速計算的原理。根據(jù)算法原理，編寫程序，成為帶支腿浮式結構水動力系數(shù)快速計算子模塊。

2 算例介紹

選取某型帶支腿浮式結構為研究對象，結構長度約為130 m，寬38 m，型深8 m，吃水4 m。將主體平臺簡化為一個浮箱，6個升降腿分布于全船兩側，每個樁腿高68 m，截面為4 m×4 m方形。使用格林函數(shù)分析方法，分兩種方案計算水下樁腿長度為10 m，20 m，30 m，40 m和50 m時該型帶支腿浮式結構的水動力系數(shù)。

第一種方案采用分工況法計算。按樁腿長度建立5種工況分別進行水動力建模。如圖5所示，最短為工況1，最長為工況5（工況2～4圖略）。所研究的浮式結構是對稱結構，只需要定義一半濕面單元，即可進行水彈性分析。這樣，工況1就有265個面元；工況2時有289個面元；工況3為313個面元；工況4為337個面元；工況5為361個面元。使用水彈性程序對五種工況分別進行數(shù)值計算。

第二種方案采用以矩陣縮減為核心的快速算法計算。以浮式結構的樁腿向水下伸出50 m，即工況5為初始狀態(tài)，快速計算當水下樁腿長度為40 m，30 m，20 m，10 m時的水彈性響應。

表2 樁腿狀態(tài)編號及濕面元編號Tab.2 The numbering schedule of the wet panels

整個浮式結構處于初始狀態(tài)時共有361個濕面元，對已經(jīng)劃分好的濕面元網(wǎng)格進行編號，浮式結構主體部分編號為1～239，對樁腿部分面元主要按照由淺及深，由大及小的原則。表2說明了算例中對樁腿編號的規(guī)律。處于樁腿50 m處的濕面元組的編號為355～361，其狀態(tài)定為1；處于樁腿10 m處的面元組編號為239～262，其狀態(tài)編號為5。

3 結果檢驗

3.1 計算精度分析

如圖6～11，在坐標平面上，有一系列的散點，每一個點的橫坐標表示快速算法結果，縱坐標表示普通算法結果，通過衡量圖中點與直線y=x的位置來判斷快速算法結果和一般算法結果的相符程度，其中（a）圖為附加質量，（b）圖為附加阻尼。每一幅圖中各有四種顏色的點，分別代表帶支腿浮式結構的樁腿處于四種長度，在圖的右下方的方框內進行了說明。每一幅圖左上方標有一個4×3的數(shù)據(jù)矩陣。其中一行代表一個固定的樁腿水下長度；一列代表兩組計算結果定量的誤差分析結果，共有三列，分別為平均相對誤差（EAVE），均方相對誤差（RMVE）和相關系數(shù)（CC）。其計算方法為：

對比分析通過快速算法得到的前六階剛體模態(tài)的水動力系數(shù)。如圖6，圖7和圖11所示，縱蕩、橫蕩和艏搖的水動力系數(shù)散點與直線y=x重合程度極高，這說明對這三種剛體模態(tài)來說，快速算法結果和普通算法結果相符程度極高。計算得到的平均相對誤差、均方相對誤差和相關系數(shù)上也證實了這樣的結論，這三種模態(tài)的平均相對誤差和均方相對誤差都在0.5%以下，而兩者之間的相關系數(shù)都在99%以上，可以認為采用快速算法得到的結果與普通算法結果一致；

和前面這三個剛體模態(tài)相比，能夠看到垂蕩、橫搖和縱搖三個模態(tài)相對應圖中的散點相對直線y=x有明顯的偏移，如圖8，圖9和圖10所示，這表明快速算法結果和普通算法結果有一定的誤差存在，從數(shù)據(jù)上看，這幾個模態(tài)的平均相對誤差、均方相對誤差幾乎都在1%以上，樁腿處于10 m時的橫搖附加阻尼系數(shù)誤差達到了7%，如圖9（b）所示，誤差明顯增高；注意到圖中的散點位置高于直線y=x，如圖10（a）所示，這說明，和普通算法相比較，快速算法在一定程度上低估了橫搖、縱搖和垂蕩的附加質量和附加阻尼。

同時，這三幅圖中誤差統(tǒng)計量表明，快速算法結果的誤差隨著樁腿長度變短而增大，如圖9（b）中的數(shù)據(jù)，當樁腿為40 m時，平均相對誤差為1.4597%和0.3421%，當樁腿長度處于10 m狀態(tài)時，這兩個數(shù)據(jù)增加到7.1231%和1.4454%。原因在于普通算法中針對樁腿不同長度時的浮式結構分別進行水動力建模，樁腿的底部加了底板，使結構物封閉。而在快速算法中，在進行矩陣拆解時，這部分濕面元沒有參與計算，不包含樁腿底板的信息。當樁腿長度較長時，浮式結構的整體濕面積相對較大，底板的影響偏??；而當樁腿長度為10 m時，盡管樁腿底板的面積并沒有變化，但是浮式結構的整體濕面積變小了，樁腿底板在濕面積中占有的比例變大，因此底板影響就顯現(xiàn)出來了。

3.2 計算效率分析

表3給出了分工況方法和快速算法耗費的CPU時間統(tǒng)計。在算例中一共分析了200個頻率，20個模態(tài)（包含14個彈性體模態(tài)）。

表3 CPU計算時間統(tǒng)計表Tab.3 Statistics of CPU time

當使用分工況法進行分析時，計算工況5耗時67 min；從工況1到工況4的計算時間和為208 min；當使用快速算法進行分析時，工況5被選為初始狀態(tài)，分析時間與傳統(tǒng)方法一致，但工況4到工況1的結果為一次性算出，消耗CPU時間為16 min，和前一種分析方法相比較，一共節(jié)約192 min，效率提高了92.3%。這還只是統(tǒng)計的在計算機上的絕對計算時間，尚不包括網(wǎng)格劃分、參數(shù)定義等前期手工準備時間。計算效率有顯著的提升。效率提高的主要原因在于后續(xù)的分析主要是在進行矩陣的縮減，實際上是在對電腦中已有的數(shù)據(jù)進行行變換和列變換，而之前方法的計算主要是進行格林函數(shù)分析和格林函數(shù)偏導數(shù)的分析，這樣節(jié)約出了大量的機時。若濕面網(wǎng)格劃分更多，格林函數(shù)矩陣的規(guī)模更大，計算效率或許將進一步提高，算法的優(yōu)勢將進一步體現(xiàn)。

4 結 論

本文在三維勢流理論體系下，利用格林函數(shù)矩陣的性質，針對帶支腿浮式結構變參數(shù)的水動力系數(shù)計算問題提出了一種相對較高效的方法。將求解格林函數(shù)的過程視為格林函數(shù)矩陣G和影響矩陣C矩陣規(guī)模減小的過程，即水下樁腿任意長度為Li時的格林函數(shù)矩陣和影響矩陣中包含著L＜Li時的信息。這樣在進行下一步計算的時候，可以通過上述已經(jīng)確定下來的編號找到退出計算的網(wǎng)格，對其進行消除，而將已有的信息保留下來，而不必重復計算，從而達到簡化計算的目的。本文還給出了通過快速算法計算結果的誤差和快速算法的效率。

結果表明：快速算法所計算的剛體的縱蕩、橫蕩和艏搖三種水動力系數(shù)與普通算法結果一致，可完全反應這些模態(tài)的水動力特性；快速算法所計算的垂蕩、橫搖和縱搖三種剛體模態(tài)的結果與普通算法有較小的誤差，但基本可以反應這些模態(tài)的水動力特性；當使用格林函數(shù)矩陣縮減法進行分析以后，算例的計算效率可以提高92%以上，非?？捎^。