基于多剛體離散元法的深海采礦揚(yáng)礦管線計(jì)算分析

基于多剛體離散元法的深海采礦揚(yáng)礦管線計(jì)算分析

典型的深度開采系統(tǒng)主要由采礦船、主采板管、中間艙、插管和采礦機(jī)組成。深海采礦整體聯(lián)動(dòng)系統(tǒng)的動(dòng)力學(xué)分析是研究采礦系統(tǒng)在開采作業(yè)的聯(lián)動(dòng)狀態(tài)下采礦船和集礦機(jī)的各單體運(yùn)動(dòng)和聯(lián)合運(yùn)動(dòng)狀態(tài)、揚(yáng)礦硬管的管線形態(tài)和情況、揚(yáng)礦軟管的空間形態(tài)和情況,采礦整體系統(tǒng)在開采作業(yè)中采礦船、硬管、中間艙、軟管、集礦機(jī)的同步運(yùn)動(dòng)性能和穩(wěn)定性,對(duì)可能出現(xiàn)的問題及應(yīng)對(duì)措施和處理效果進(jìn)行分析和模擬,為系統(tǒng)的研制和開采(試驗(yàn))作業(yè)提供理論依據(jù)和技術(shù)參考。深海采礦系統(tǒng)的揚(yáng)礦管線長達(dá)數(shù)千米,在變形過程中呈現(xiàn)出大位移、小變形的幾何非線性特性。國內(nèi)外深海采礦管線計(jì)算分析大多是基于有限元的思想,其分析方法分為2種:更新拉格朗日法(UL)和隨體旋轉(zhuǎn)法(CR)。二者均采用增量加載平衡迭代法求解。UL法因?yàn)榻Y(jié)構(gòu)運(yùn)動(dòng)中包含剛體旋轉(zhuǎn)成分,在大位移、大轉(zhuǎn)動(dòng)中它們?yōu)橛邢蘖?因此,在采用增量法求解時(shí),其加載步必須很小,小到可以按矢量操作,否則剛體運(yùn)動(dòng)的誤差積累,造成發(fā)散;CR法的坐標(biāo)系是粘附于結(jié)構(gòu)上,隨后一起平移和旋轉(zhuǎn)運(yùn)動(dòng)的剛性坐標(biāo)系。因此,梁單元相對(duì)于該坐標(biāo)系的運(yùn)動(dòng)為純粹的變形運(yùn)動(dòng),從而能將剛體運(yùn)動(dòng)從總位移中消除,使其余變形運(yùn)動(dòng)相對(duì)于局部剛性坐標(biāo)系始終較小。故可以運(yùn)用小變形梁理論進(jìn)行有限元建模,不會(huì)因剛體運(yùn)動(dòng)的誤差積累而造成發(fā)散,并大大減少增量迭代求解的程序計(jì)算量。但是,隨體剛性坐標(biāo)系的確定較困難,即使采用CR法,其增量步也必須很小。由于深海采礦系統(tǒng)非常復(fù)雜,采用有限元法進(jìn)行深海采礦整體聯(lián)動(dòng)系統(tǒng)動(dòng)力學(xué)分析時(shí),增量迭代計(jì)算的截荷步不能取得太大,否則將不收斂,計(jì)算工作量大,占用大量內(nèi)存和機(jī)時(shí)。1992年,G.G.W.Mustoe等提出采用基于多剛體動(dòng)力學(xué)的離散元方法計(jì)算分析揚(yáng)礦硬管,并對(duì)二維深海采礦揚(yáng)礦硬管進(jìn)行了動(dòng)態(tài)分析;隨后,CHENGBao-rong等采用該方法對(duì)三維深海采礦揚(yáng)礦硬管進(jìn)行結(jié)構(gòu)模態(tài)分析。在此,作者采用離散元方法對(duì)整體采礦系統(tǒng)進(jìn)行三維動(dòng)力學(xué)計(jì)算與分析。1動(dòng)力學(xué)分析方法深海采礦系統(tǒng)的揚(yáng)礦硬管由若干根硬管用管接頭連接而成,每一段硬管不能變形,可以視為剛性構(gòu)件,但長達(dá)數(shù)千米的管線,已經(jīng)不能視為一個(gè)剛性構(gòu)件,必須考慮其大位移、小變形的幾何非線性特性。對(duì)于柔性軟管,更應(yīng)采用基于大位移和旋轉(zhuǎn)的非線性分析技術(shù)進(jìn)行分析、計(jì)算?；诙鄤傮w動(dòng)力學(xué)理論,采用離散元方法(discreteelementmethod,DEM)對(duì)上述管線系統(tǒng)(包括硬管、中間倉和軟管)進(jìn)行分析。其基本思想是:將管結(jié)構(gòu)離散成有限剛段,以其表示物體的慣量特性,而其剛度向接點(diǎn)等效移植,外力向段的質(zhì)心移植。如圖1所示,將管結(jié)構(gòu)離散為n個(gè)等長度且不發(fā)生彈性變形的剛性單元(剛體元),相鄰的剛體元之間用無質(zhì)量、無長度的柔性聯(lián)接元聯(lián)接。該柔性聯(lián)接元等效于如下6個(gè)彈簧:a.1個(gè)縱向彈簧,kx=EA/L;b.2個(gè)橫向彈簧,ky=12EIz/L3,kz=12EIy/L3;c.1個(gè)扭轉(zhuǎn)彈簧,krx=GIp/L;d.2個(gè)彎曲彈簧,kry=EIy/L,krz=EIz/L。式中:kx,ky,krx,kry,krz分別為相應(yīng)彈簧的剛度;E為彈性模量;G為剪切模量;A為剛體元的截面積;L為剛體元的長度;Iy和Iz分別為截面對(duì)應(yīng)于y軸和z軸的慣性矩;Ip為極慣性矩。在動(dòng)力學(xué)分析中,剛體元只有動(dòng)能而無彈性勢(shì)能;而柔性聯(lián)接元只有彈性勢(shì)能而無動(dòng)能。管結(jié)構(gòu)的變形能僅儲(chǔ)存在柔性聯(lián)接元中,剛體元中任一點(diǎn)的位移可用剛體元質(zhì)心的剛體位移來描述。取剛體i的質(zhì)心在慣性參考系中的3個(gè)直角坐標(biāo)和反映剛體方位的3個(gè)歐拉角作為笛卡兒廣義坐標(biāo),即qi=[xyzψθφ]Ti。采用Lagrange乘子法建立系統(tǒng)動(dòng)力學(xué)方程為:ddt(?Τ?˙q)Τ-(?Τ?q)Τ+fΤqρ=Q。(1)完整約束方程為:f(q,t)=0。其中:T為系統(tǒng)動(dòng)能列陣;q為系統(tǒng)廣義坐標(biāo)列陣,q=[qΤ1qΤ2…qTn]T;Q為廣義力列陣;ρ為對(duì)應(yīng)于完整約束的拉氏乘子列陣;˙q為系統(tǒng)廣義速度列陣;f為完整約束方程列陣。將式(1)寫成更一般的形式:{F(q,u,˙u,λ,t)=0,G(u,˙q)=u-˙q=0,Φ(q,t)=0。(2)式中:u為系統(tǒng)廣義速度列陣;˙u為廣義加速度列陣;λ為約束反力及作用力列陣;F為系統(tǒng)動(dòng)力學(xué)微分方程列陣;G為描述廣義速度的代數(shù)方程列陣;Φ為描述約束的代數(shù)方程列陣。推導(dǎo)出上述二階常微分方程和代數(shù)方程后,應(yīng)用Gear剛性積分算法對(duì)式(2)進(jìn)行數(shù)值求解。首先,定義系統(tǒng)的狀態(tài)矢量y=[qTuTλT]T;然后,根據(jù)當(dāng)前時(shí)刻的系統(tǒng)狀態(tài)矢量值,用Taylor級(jí)數(shù)預(yù)估下一個(gè)時(shí)刻系統(tǒng)的狀態(tài)矢量值:yn+1=yn+?yn?th+12!?2yn?t2h2+?。(3)式中:時(shí)間步長h=tn+1-tn。采用這種預(yù)估算法得到新的時(shí)刻的系統(tǒng)狀態(tài)矢量值通常不準(zhǔn)確,可采用Gear法k+1階積分進(jìn)行校正。如果預(yù)估值與校正值的差值小于規(guī)定的積分誤差限,接受該解,進(jìn)行下一時(shí)刻的求解。否則拒絕該解,并減小積分步長,重新進(jìn)行預(yù)估-校正過程,直到求解時(shí)間達(dá)到規(guī)定的模擬時(shí)間為止。采用多剛體離散元法所得到的模型具有簡單性和力學(xué)直觀性的特點(diǎn),可以用多剛體動(dòng)力學(xué)理論描述單元的大位移和大轉(zhuǎn)動(dòng)。由于組裝在單元交界點(diǎn)進(jìn)行,因此,總剛度矩陣為一個(gè)高度稀疏的矩陣,只與單元編號(hào)排序有關(guān),而與結(jié)點(diǎn)編號(hào)無關(guān),將顯著減少計(jì)算量,節(jié)省計(jì)算機(jī)內(nèi)存。2計(jì)算與文本測(cè)試從靜態(tài)分析和動(dòng)態(tài)分析兩方面分別選擇典型算例,檢驗(yàn)上述多剛體離散元模型的正確性與可行性。2.1應(yīng)力分析結(jié)果如圖2所示圓形截面懸臂梁,彈性模量E=207GPa;長度L=0.61m;直徑d=0.025m。在端點(diǎn)處存在3種載荷:扭矩M=22.5N·m;豎向載荷Fv=-111N;軸向載荷Fa=445N。分以下3種情況對(duì)其進(jìn)行應(yīng)力分析:a.僅有扭矩M作用時(shí)固接處最大剪切應(yīng)力τ;b.當(dāng)豎向載荷Fv與軸向載荷Fa同時(shí)作用時(shí)固接處最大正應(yīng)力σ;c.扭矩M,豎向載荷Fv與軸向載荷Fa三者同時(shí)作用時(shí)固接處最大主應(yīng)力σ。表1所示為采用離散元法計(jì)算得到的固接處應(yīng)力與解析法的對(duì)比。可見,離散元算法的計(jì)算精度較高。2.2載荷作用下的共振設(shè)簡支梁的直徑為6mm,長度為1.216m,彈性模量E=200GPa,質(zhì)量密度ρ=7.8×103kg/m3。在簡支梁中點(diǎn)處分別施加2種載荷:正弦載荷F=2.5sin(10t)及正弦載荷F=2.5sin(51t),考察中點(diǎn)處的響應(yīng)。由于該簡支梁系統(tǒng)的第1階固有頻率為51rad/s,所以,在第2種載荷作用下將產(chǎn)生共振。圖3和圖4所示為采用解析法和離散元法所得結(jié)果?？梢?計(jì)算值與解析值較吻合。3計(jì)算機(jī)模擬結(jié)果應(yīng)用多剛體離散元法對(duì)中國大洋多金屬結(jié)核1km海試系統(tǒng)進(jìn)行整體系統(tǒng)聯(lián)動(dòng)動(dòng)力學(xué)分析。海試系統(tǒng)構(gòu)成與參數(shù)詳見文獻(xiàn)。應(yīng)用多剛體離散元法將長度為900m的揚(yáng)礦硬管離散為90個(gè)剛體元,相鄰剛體元用柔性聯(lián)接元相連;將長度為400m的軟管離散為200個(gè)剛體元,相鄰剛體元采用柔性聯(lián)接元相連。揚(yáng)礦硬管上端通過球鉸與采礦船相連,軟管下端通過十字鉸與集礦機(jī)相連。每段管線系統(tǒng)所受液動(dòng)力采用莫里森公式計(jì)算,并施加到每個(gè)剛體元的質(zhì)心。設(shè)海流沿-X方向或+X方向流動(dòng),海洋表面海流速度為0.772m/s,海底海流速度為0.15m/s。集礦機(jī)與采礦船朝+X方向聯(lián)動(dòng),且采礦船同時(shí)做升沉運(yùn)動(dòng),幅度為4m,周期為8s。針對(duì)如下2種典型工況進(jìn)行計(jì)算機(jī)模擬:a.集礦機(jī)與采礦船均以0.5m/s的速度沿+X軸方向順流運(yùn)行250m(海流朝+X方向流動(dòng));b.集礦機(jī)與采礦船均以0.5m/s的速度沿+X軸方向逆流運(yùn)行250m(海流朝-X方向流動(dòng))。圖5所示為中間艙相對(duì)采礦船的位置變化曲線。可見,在a.工況下,中間艙X方向的初始偏移為-9.1m,系統(tǒng)運(yùn)行250m后中間倉偏移為-22.5m;在b.工況下,中間艙沿X方向初始偏移為-13.8m,系統(tǒng)運(yùn)行250m后中間倉偏移為-51.3m。圖6所示為集礦機(jī)沿X方向所受軟管作用力變化曲線。可見,當(dāng)系統(tǒng)順流運(yùn)行250m后集礦機(jī)沿X方向所受軟管作用力穩(wěn)定在4.1kN左右;當(dāng)系統(tǒng)逆流運(yùn)行250m后集礦機(jī)沿X方向所受軟管作用力穩(wěn)定在-4.3kN左右。圖7所示為集礦機(jī)沿Z方向受軟管作用力變化曲線。可見,集礦機(jī)沿Z方向所受軟管作用力約在17.6kN附近變化,且順流時(shí)集礦機(jī)沿Z方向受力比逆流時(shí)集礦機(jī)沿Z方向受力稍大。為了驗(yàn)證離散元方法的效率,采用有限元法對(duì)上述2個(gè)工況進(jìn)行分析。表2所列為采用2種方法計(jì)算所得中間艙相對(duì)采礦船偏移量與計(jì)算耗時(shí)的對(duì)比。計(jì)算均在相同配置的PC機(jī)(PentiumIV1.8GHz,內(nèi)存512M)上完成。由表2可見,在計(jì)算誤差僅為1%～2%時(shí),離散元方法的計(jì)算速度較有限元方法的計(jì)算速度提高了近8倍。以上研究結(jié)果表明,對(duì)深海長管線這類形狀較為簡單的桿系結(jié)構(gòu),在不需要特別討論應(yīng)力和應(yīng)變等情況時(shí),采用離散元法進(jìn)行動(dòng)力學(xué)分析是