佟 瑩， 夏 健

(南京航空航天大學 航空學院，南京 210016)

1 引 言

管道提升式深海采礦系統(tǒng)以其采集效率和技術(shù)可行性等優(yōu)勢成為最具發(fā)展前景的深海采礦模式。為適應復雜的深海采礦環(huán)境，需采取穩(wěn)定性好和環(huán)境適應性強的管道提升系統(tǒng)。全柔性立管系統(tǒng)[1,2]無需多級泵組推送，結(jié)構(gòu)精簡，系統(tǒng)靈活性高，大長度海洋柔性管道對復雜海底環(huán)境適應性強，采礦作業(yè)范圍廣，且柔性管路材料的耐腐蝕性和柔韌性有助于提升采礦系統(tǒng)使用壽命，具有廣闊的工程應用前景[3]。但全柔性管道提升式深海采礦系統(tǒng)發(fā)展時間較短，其理論研究還不成熟。

目前，對于結(jié)構(gòu)大變形及其非線性特性的研究廣泛采用有限單元方法。肖林京等[4]采用T.L.列式法描述深海揚礦鋼管系統(tǒng)的力系平衡，建立幾何非線性有限元基本方程，通過采礦船加速、勻速和減速拖航時揚礦管的偏移構(gòu)型、橫向速度、加速度及彎曲應力等參數(shù)表征其動力學響應和非線性動態(tài)特性。郭小剛等[5]基于懸鏈線理論對深海揚礦鋼管系統(tǒng)底部軟管段進行了非線性靜平衡分析。Garcia-Palacios等[6]采用U.L.列式法描述力系平衡并建立非線性有限元基本方程，基于ANSYS軟件框架編寫結(jié)構(gòu)模型的數(shù)值計算程序，對聚乙烯管的布放入水過程進行了幾何非線性動力學分析。馬今偉等[7]采用U.L.列式法基于非線性無額外自由度廣義有限元法對彈塑性大變形靜力學問題進行了非線性分析。樊玉新等[8]使用張量體系描述的非線性有限元模型分析柔性體大變形過程。與傳統(tǒng)有限元方法相比，基于連續(xù)介質(zhì)力學思想，張量體系描述的非線性有限元模型的模化和計算過程更加精簡，更便于代碼開發(fā)。張量體系描述的隨體坐標與等參單元的自然坐標之間的一致性在連續(xù)介質(zhì)和有限單元模型之間建立了光滑過渡。此外，克服了結(jié)構(gòu)運動小變形假設關于位移和變形的約束，適應性更廣。

基于三維固體有限變形理論[9]，本文采用張量體系描述結(jié)構(gòu)運動和力系平衡，建立全柔性立管系統(tǒng)的數(shù)學模型，推導了系統(tǒng)運動控制方程，針對幾何非線性和非保守載荷等關鍵問題提出了有效處理方案。采用非線性有限元法處理計算模型，通過工程算例的模擬結(jié)果驗證理論方法和計算模型的有效性和實用性。

2 物理模型

柔性立管系統(tǒng)將海面上靜止的能量供應母站與海底自推進集礦系統(tǒng)直接聯(lián)接，通過單級泵經(jīng)柔性立管將礦漿泵送至水面，依靠靈活柔性管道適應海底復雜作業(yè)環(huán)境，原理如圖1所示。管道大變形屬幾何非線性問題，海洋環(huán)境載荷隨管道構(gòu)型變化而變具有非保守力屬性，耦合雙重非線性特征成為發(fā)展數(shù)學模型的關鍵問題。

3 數(shù)學模型

本文從連續(xù)介質(zhì)力學有限變形理論出發(fā)，建立張量體系描述的非線性有限元模型。將柔性管道簡化為索單元，管路位形由其軸線位置描述，采用T.L.列式法，取t=0時刻無應變和位形已知的初始構(gòu)型C0為參考構(gòu)型，結(jié)構(gòu)運動和力系虛功(慣性力、內(nèi)力和外力)均相對于C0定義。取能量共軛的第二類Piola-Kirchhoff應力張量和Green應變張量度量結(jié)構(gòu)應力和應變，由達朗貝爾原理[10]建立連續(xù)系統(tǒng)動力學平衡方程，

(1)

式中L0和A0分別為C0索單元的初始長度和截面面積?；谑?1)，柔性管道運動的位移和應變不受小變形假設約束。

圖1 柔性管道提升式采礦系統(tǒng)原理

采用隨體坐標ξi描述管道大變形，初始時刻長度為ds0的線元dr=gidξi，當前構(gòu)型C變形后線元dR=Gidξi的長度為ds。gi和Gi為材料點ξi在C0和C的協(xié)變基矢量，為空間點位的函數(shù)。張量體系中Green應變張量E可由線元的長度變化定義，

(2)

式中g(shù)11和G11為C0和C的度量張量的協(xié)變分量，g11=g1·g1，G11=G1·G1。E11為位移的非線性函數(shù)。

第二類Piola-Kirchhoff應力張量S的逆變分量為

S11=Eg11g11E11

(3)

式中E為彈性模量，g11為C0度量張量的逆變分量，g11g11=1。同理，S11為位移的非線性函數(shù)。

作用于大變形柔性管道的海洋環(huán)境載荷具有非保守力屬性。將C中隨時間變化的線元AdL矢量映射至C0得與時間無關的edL0，取T為流體面力矢量，則有

(4)

將有限元列式的單元平衡方程匹配至全局[12]得矩陣形式的系統(tǒng)運動微分方程，

(5)

4 計算模型

4.1 非線性靜力學方程求解

非線性靜力學方程為內(nèi)力與外力間的平衡，外力作用下張緊管線獲取抵抗外力的結(jié)構(gòu)剛度。內(nèi)力向量R為節(jié)點位移向量U的非線性函數(shù)[8]，對R進行Taylor序列展開至線性項，定義非線性的切線剛度矩陣K，

Rn + 1=Rn+[Kn] ΔUn + 1

(6)

由逐步加載與Newton-Raphson迭代混合的數(shù)值求解方法抑制非線性方程的數(shù)值解漂移及迭代計算不收斂等問題。混合法保證每次增量加載引起的位形變化滿足特定誤差精度，且能有效抑制純增量加載的誤差累計。采用粘性松弛技術(shù)[9]引入人工阻尼D，以保證初始無應變狀態(tài)下剛度矩陣K非奇異。第k+1個增量載荷作用下，第n+1次迭代，系統(tǒng)的靜力學平衡方程為

(7)

為加速收斂，人工阻尼矩陣D隨增量推進滿足如下關系，

(8)

[D]=μk + 1[D]

(9)

式中μ為阻尼因子，γ為衰減因子。阻尼因子μ隨著載荷增量推進而變化。衰減因子作為粘性松弛技術(shù)的控制參數(shù)，可視求解問題而定。本文靜力學分析算例中，初始阻尼因子為0.01，衰減因子為0.1，初始阻尼元素初值為 1.0e +4。

4.2 非線性動力學方程求解

為有效求解系統(tǒng)非線性動力學響應，本文采用廣義α隱式時間推進格式[10]，綜合了Newmark法、HHT-α法和WBZ-α法的特征，通過調(diào)整數(shù)值參數(shù)保證高頻段具備高速衰減的數(shù)值阻尼過濾虛假高階振型，而在低頻段控制較小的數(shù)值阻尼保證數(shù)值解的精度?；趶V義α格式，t+Δt時刻節(jié)點的速度和位移滿足如下關系,

(10)

式中數(shù)值參數(shù)γ和β滿足關系(11，12)時，算法無條件穩(wěn)定且具備二階精度,

(11,12)

αm和αf為插值系數(shù)，經(jīng)t和t+Δt時刻參量插值，t+Δt時刻節(jié)點的運動學參量及外載修正如下,

Ut + Δ t=(1-αf)Ut + Δ t+αfUt

Ft + Δ t=(1-αf)Ft + Δ t+αfFt

(13)

直接積分法通過構(gòu)造t+Δt時刻位移與t時刻已知變量的遞推表達式完成時間推進。由式(13)可推導得系統(tǒng)等效線性方程組

[Keff] ΔU=Feff

(14)

式中

本文采用Rayleigh阻尼矩陣C將實際結(jié)構(gòu)阻尼矩陣簡化為質(zhì)量矩陣M和剛度矩陣K的線性組合，系數(shù)cM和cK由實驗給出,

C=cMM+cKK

(15)

本文取cM和cK均為5.0e -4。

5 非線性有限元程序驗證

基于上述計算模型，采用FORTRAN語言編寫非線性有限元數(shù)值計算程序。目前，認為Verification & Validation[13]是評估非線性問題數(shù)值計算模擬可靠性的最客觀方式，其驗證過程如圖2所示。

圖2 模型驗證和確認的工作范圍

為了驗證本文方法處理幾何非線性問題的有效性，采用網(wǎng)格收斂指數(shù)(GCI)方法[14]客觀評定漸近解的離散誤差范圍。初始長度L=0.35 m，線密度ρA=0.314 kg/m，重力加速度g=2.0 m/s2，彈性模量E=1.4e+6 Pa的均質(zhì)鏈條自0.1π位置釋放的自重驅(qū)動懸掛鏈條擺動，時間步長為 5.0e -3 s。圖3展示了基于三組網(wǎng)格密度(10 elems，20 elems，40 elems)，不同時刻滿足收斂標準的數(shù)值模擬結(jié)果及自由端位移和速度的GCI指數(shù)。數(shù)值結(jié)果表明，40 elems網(wǎng)格密度下，本文計算模型對幾何非線性問題具有穩(wěn)定的預測能力。速度場漸近解的離散誤差范圍大于位移場，但收斂指數(shù)仍在置信范圍內(nèi)，說明計算模型能夠有效處理數(shù)學模型的非線性動力學特征。

圖3 計算模型驗證

取該算例的實驗與解析解[15]及現(xiàn)有數(shù)值模擬結(jié)果[16]為參考，驗證本文方法計算結(jié)果的準確性。圖4為自釋放起2.4 s內(nèi)鏈條構(gòu)型隨時間變化歷程，右側(cè)自由端的踢腿現(xiàn)象與實驗觀察結(jié)果一致，說明張量體系描述的非線性有限元法成功反映了張緊力為結(jié)構(gòu)運動提供剛度的真實物理現(xiàn)象。此外，相比于常規(guī)有限元法忽略剛度引起抖動的踢腿，本文數(shù)值結(jié)果更加光滑，說明本文計算模型能夠有效處理速度梯度過大引發(fā)數(shù)值計算不穩(wěn)定的問題。由10次往復擺動歷程數(shù)值結(jié)果計算平均周期(T=2.196 s)和特征參數(shù)gT2/L=27.556(由解析解計算其理論值為27.3[14])。數(shù)值結(jié)果與理論解的一致性驗證了本文計算模型的準確性和有效性。

圖4 自重驅(qū)動鎖鏈構(gòu)型變化時間歷程

6 柔性管道非線性特性分析

以完成410 m淺水水試的全柔性立管系統(tǒng)[1]為例，采用本文數(shù)值計算程序分兩階段對柔性管路非線性響應特性進行分析。柔性管路長500 m，內(nèi)外徑為0.4445 m和0.508 m，彈性模量為1.4 GPa，單位長度管線濕重為3336.9 N/m，揚礦礦漿混合液密度為1130 kg/m3，海水密度為1030 kg/m3。自尾端10 m位置起裝有80 m范圍的浮力裝置，結(jié)構(gòu)如圖1 所示。離散有限單元模型由100個等參索單元構(gòu)成，水面母站聯(lián)接Node 1，集礦機聯(lián)接Node 101。

6.1 柔性管道非線性靜力學分析

第一階段為柔性立管系統(tǒng)穩(wěn)定作業(yè)狀態(tài)的靜力學分析。圖5為水面母站靜止，集礦機與母站之間水平距離X=40 m，80 m，120 m，160 m，200 m和240 m時，柔性立管系統(tǒng)的靜平衡構(gòu)型。X較小時，浮力塊作用下海底管路呈S型剖面，該模擬結(jié)果與水試實驗結(jié)果一致[1]。隨X增大，S型剖面逐漸拉直，富余管線的收縮特性使得海洋柔性管道具備靈活適應復雜海底環(huán)境的優(yōu)勢。

圖5 集礦機不同位置對應的柔性立管系統(tǒng)靜平衡構(gòu)型

圖6為設置浮力塊的局部管線在靜平衡構(gòu)型中的軸向應力。無浮力設備時，管線軸向應力自海面向下線性降低，海底富余長度的管線由于較小的軸向應力而出現(xiàn)局部松弛。浮力設備的施加調(diào)整了相應位置的軸向應力分布，有效改善了海底富余管路的松弛現(xiàn)象。緊湊的S型剖面的軸向應力呈V型分布，存在最小應力的拐點。隨著S型剖面逐漸拉直，管路軸向應力呈線性增加，過渡拐點消失。

圖6 不同靜平衡構(gòu)型中浮力施加節(jié)點處的軸向應力

6.2 柔性管道入水布放過程非線性動力學分析

第二階段為柔性立管系統(tǒng)入水布放過程的動力學分析。柔性管道一端固定于水面靜止的母站，另一端聯(lián)接濕重8.5 t的集礦系統(tǒng)。釋放至集礦系統(tǒng)著床的柔性管道系統(tǒng)入水過程中，構(gòu)型隨時間變化歷程如圖7所示。集礦系統(tǒng)自水面釋放經(jīng)歷8.75 s接觸海床面，假定集礦機接觸海床面以后位置不再變化，落地位置與水面母站的水平位移為226.87 m。數(shù)值結(jié)果表明，浮力設備及海洋環(huán)境載荷對集礦系統(tǒng)的著床位置有顯著影響。圖8為集礦系統(tǒng)接觸海床面后，柔性立管四個局部位置的節(jié)點坐標隨時間變化歷程。在海洋環(huán)境載荷的非保守力作用下，柔性立管動力學響應呈類周期性。

圖7 柔性立管系統(tǒng)構(gòu)型隨時間變化歷程

圖8 局部節(jié)點位置隨時間變化歷程

7 結(jié) 論

本文基于三維固體有限變形理論建立全柔性立管系統(tǒng)的數(shù)學模型，采用非線性有限元方法得到了非保守載荷作用下幾何非線性柔性立管系統(tǒng)的數(shù)值計算模型，開發(fā)相應計算軟件應用于柔性立管系統(tǒng)靜動力學問題的數(shù)值模擬。有限變形理論對結(jié)構(gòu)位移和應變大小不做任何限制，嚴格考慮幾何非線性因素并代入公式推導，廣義α隱式時間推進格式在數(shù)值計算過程中有效處理高頻段和低頻段的數(shù)值阻尼，提升了數(shù)值模擬軟件對非線性動力學問題的處理能力和適應性。本文軟件有望處理柔性立管系統(tǒng)工程實際中更加復雜的非線性動力學問題。