求解大柔度梁/繩索的ANCF單元及隱式迭代格式

羅 鑫，魏泳濤

(1.四川大學(xué) 建筑與環(huán)境學(xué)院，成都 610065; 2.四川中銳信息技術(shù)有限公司，成都 610094)

1 引 言

繩索結(jié)構(gòu)應(yīng)用廣泛，如太空中的系留衛(wèi)星系統(tǒng)(Tethered Satellites System)[1,2]、海洋船舶工業(yè)、核工業(yè)和機(jī)場(chǎng)等領(lǐng)域內(nèi)的橋式起重機(jī)(Overhead Cranes)[3,4]以及高速列車(chē)系統(tǒng)中的懸鏈線(xiàn)(Pantograph/Cantenary)[5]等。彈簧-阻尼-集中質(zhì)量模型[6-8]是數(shù)值模擬繩索的常用方法，但該法的繩索處理成分段折線(xiàn)而不再保持光滑曲線(xiàn)的形態(tài)，因而精度較低。有限元法在分析繩索時(shí)[9,10]，通常采用具有極低抗彎剛度或大柔度的梁?jiǎn)卧员ＷC其形態(tài)的光滑性。在繩索運(yùn)動(dòng)中，往往存在大范圍剛體運(yùn)動(dòng)和大變形的耦合作用。但傳統(tǒng)梁?jiǎn)卧墓?jié)點(diǎn)自由度僅反映變形的撓度和轉(zhuǎn)角，即單元的剛體運(yùn)動(dòng)和變形是解耦的。浮動(dòng)坐標(biāo)法[11]將繩索運(yùn)動(dòng)考慮成隨浮動(dòng)坐標(biāo)系的大范圍運(yùn)動(dòng)與相對(duì)于浮動(dòng)坐標(biāo)系的變形的疊加，但由此導(dǎo)出的科氏力和向心力等慣性力表達(dá)式十分復(fù)雜。

為克服浮動(dòng)坐標(biāo)法中的難點(diǎn)，Shabana[12]提出了以節(jié)點(diǎn)坐標(biāo)和斜率作為節(jié)點(diǎn)自由度的絕對(duì)節(jié)點(diǎn)坐標(biāo)方法ANCF(Absolute nodal coordinate formulation)。該方法的基本思想為在連續(xù)介質(zhì)力學(xué)中對(duì)即時(shí)構(gòu)形進(jìn)行物質(zhì)描述時(shí)，運(yùn)動(dòng)和變形本身就是耦合的。因此，若將梁/繩索的節(jié)點(diǎn)坐標(biāo)和斜率而非節(jié)點(diǎn)撓度和轉(zhuǎn)角設(shè)置為節(jié)點(diǎn)自由度，則可自動(dòng)耦合大范圍運(yùn)動(dòng)和大變形，也無(wú)需單獨(dú)考慮科氏力和向心力等慣性力，且質(zhì)量矩陣是常數(shù)矩陣。 ANCF 已有關(guān)于繩索的研究，Gerstmayr等[13]將標(biāo)準(zhǔn)梁?jiǎn)卧淖杂啥葴p少并引入附加的形函數(shù)，構(gòu)造了低階的空間兩節(jié)點(diǎn)繩索單元；Seo等[14]使用該單元模擬了高速列車(chē)系統(tǒng)接觸網(wǎng)的三維大變形問(wèn)題；Shan等[15]使用無(wú)抗彎剛度的該單元分析了太空中的系留網(wǎng)。

ANCF中單元的節(jié)點(diǎn)等效力計(jì)算遠(yuǎn)比傳統(tǒng)有限元復(fù)雜。劉鋮等[16]基于第一類(lèi)Piola-Kirchhoff應(yīng)力張量導(dǎo)出了節(jié)點(diǎn)等效力及雅可比矩陣的解析表達(dá)式。范紀(jì)華等[17]在節(jié)點(diǎn)等效力求解中引入了彈性線(xiàn)方法。Gerstmayr等[13]逐一計(jì)算了節(jié)點(diǎn)等效力各元素，特別是對(duì)應(yīng)于彎曲變形的情況。Zemljaric等[18]對(duì)節(jié)點(diǎn)等效力計(jì)算方法進(jìn)行了一定簡(jiǎn)化。Berzeri等[19]推導(dǎo)了小變形下歐拉-伯努利梁?jiǎn)卧墓?jié)點(diǎn)等效力。這也是目前ANCF分析繩索時(shí)采用顯式積分的主要原因。顯式積分不需要迭代，既避免了不收斂的問(wèn)題，也免去了由節(jié)點(diǎn)等效力導(dǎo)出切線(xiàn)剛度矩陣的繁瑣過(guò)程，但必須按算法的穩(wěn)定原則確定足夠小的時(shí)間步長(zhǎng)。而隱式積分在各時(shí)間步上需迭代求解以滿(mǎn)足平衡方程，在采用較大時(shí)間步長(zhǎng)時(shí)也能保證精度要求。

由節(jié)點(diǎn)等效力導(dǎo)出隱式積分所需的的切線(xiàn)剛度矩陣，目前尚無(wú)此方面完整研究的報(bào)道。此外，基于ANCF方法分析繩索時(shí)，多針對(duì)單根繩索，而對(duì)工程中常見(jiàn)的多根繩索系于一點(diǎn)的情況，還未有考慮。本文對(duì)ANCF大柔度梁/繩索單元，導(dǎo)出了用于平衡迭代的切線(xiàn)剛度矩陣的全部公式；其關(guān)鍵環(huán)節(jié)在于將兩矢量的叉乘表示為反對(duì)稱(chēng)矩陣與矢量的點(diǎn)乘，以此改造ANCF大柔度梁/繩索單元的節(jié)點(diǎn)等效力公式；由此導(dǎo)出的節(jié)點(diǎn)等效力和切線(xiàn)剛度矩陣均為矩陣形式，也非常方便代碼開(kāi)發(fā)；對(duì)不同ANCF單元的鉸接和剛性連接，引入罰方法的處理；最后結(jié)合HHT(Hilber-Hughes-Taylor)隱式格式，對(duì)弦微幅和小幅振動(dòng)、雙桿擺和T字形桿以及柔繩進(jìn)行了數(shù)值模擬。

2 基于ANCF的大柔度梁/繩索單元

研究細(xì)長(zhǎng)的等截面大柔度梁/繩索，在運(yùn)動(dòng)/變形中始終維持光滑曲線(xiàn)，且橫截面是正多邊形，即橫截面對(duì)通過(guò)形心的任何軸都具有相同的慣性矩。本文基于文獻(xiàn)[13]提出的2節(jié)點(diǎn)12自由度的ANCF梁/繩索單元，該單元服從平截面假設(shè)，不考慮剪切變形和扭轉(zhuǎn)變形。當(dāng)單元具有極低或零抗彎剛度時(shí)，即可用于繩索模擬。

r=[S1IS2IS3IS4I]e=Se

(1)

式中I為3×3的單位陣,S1～S4為埃爾米特插值基函數(shù),e為將節(jié)點(diǎn)i和j處的位置矢量和切線(xiàn)矢量排列而得的單元自由度向量。

(2)

Sx,Sx x,rx和rx x為相應(yīng)的一階和二階導(dǎo)數(shù)。

單元慣性力在虛位移δr上的虛功為

(3)

(4)

圖1 單元內(nèi)任一點(diǎn)P的位置的插值

由歐拉-伯努利梁理論，單元內(nèi)力的虛功為

(5)

為簡(jiǎn)化書(shū)寫(xiě)，記矢量c=rx×rx x。將矢量叉乘運(yùn)算改寫(xiě)為矩陣乘積： -rx x×rx=Grx,rx×rx x=Qrx x,c×rx=Rrx;G,Q和R為反對(duì)稱(chēng)矩陣，對(duì)偶矢量分別為-rx x,rx和c。

(6)

(7)

(8)

3 隱式動(dòng)力學(xué)迭代格式

為將隱式積分格式應(yīng)用于前述建立的ANCF大柔度梁/繩索單元，需建立平衡迭代格式，為此應(yīng)從節(jié)點(diǎn)等效力導(dǎo)出相應(yīng)的切線(xiàn)剛度矩陣。

采用HHT隱式時(shí)間積分方法建立繩索的隱式動(dòng)力學(xué)迭代格式，根據(jù)HHT方法，

(1+α)Ft + Δt-αFt-Pt + (1 + α)Δt=0

(9)

(10)

(11)

式中Ψ為不平衡力向量，C為阻尼矩陣，P為等效外載荷。

當(dāng)積分參數(shù)滿(mǎn)足-0.3≤α≤0，β=0.25(1-α)2，γ=0.5-α?xí)r，HHT方法是無(wú)條件穩(wěn)定的。其中α=0時(shí)，HHT退化成β=0.25，γ=0.5的Newmark法，該方法的主要缺點(diǎn)是無(wú)法濾掉虛假的高階分量。

動(dòng)力學(xué)切線(xiàn)剛度矩陣(即迭代中的雅可比矩陣)為

(12)

(13)

(14)

(15)

(16)

此處略去切線(xiàn)剛度矩陣繁瑣的推導(dǎo)過(guò)程。由此在t+Δt時(shí)刻的Newton-Raphson迭代格式為

(17)

需要注意的是，HHT法在單元過(guò)程中要單獨(dú)計(jì)算上一時(shí)刻應(yīng)力的節(jié)點(diǎn)等效力，且外載荷是t+(1+α)Δt時(shí)刻的。

4 鉸接和剛性連接的處理

ANCF大柔度梁/繩索單元使用切線(xiàn)的斜率為自由度，在剛性連接中，連接處的節(jié)點(diǎn)具有不同的切線(xiàn)方向，因此應(yīng)作特別處理。

如圖2所示，節(jié)點(diǎn)i和節(jié)點(diǎn)j分別鉸接或剛性連接，約束方程為

(18)

圖2 鉸鏈和剛性連接

(19)

(20)

式中0為3×1的零矢量或3×3的零矩陣，C1和C2是罰因子，即相應(yīng)的彈簧剛度。本文算例表明，可分別取為拉伸和彎曲切線(xiàn)剛度矩陣最大元素的1000倍。切線(xiàn)剛度矩陣第三項(xiàng)的推導(dǎo)中考慮矢量模的倒數(shù)為常數(shù)，從而簡(jiǎn)化計(jì)算過(guò)程。

5 數(shù)值算例

5.1 弦的微幅和小幅振動(dòng)

弦在張緊狀態(tài)下作橫向微幅振動(dòng)時(shí)，內(nèi)部張力可視為常數(shù)。圖3所示的弦，長(zhǎng)為1 m，直徑為0.2 mm，彈性模量為200 GPa，密度為7800 kg/m3，兩端簡(jiǎn)支，弦的張力T=5 N。初始時(shí)受外部約束而水平靜置，突然撤去外部約束后，弦在重力作用下發(fā)生微幅振動(dòng)，其理論解為

(21)

圖3 振動(dòng)的弦

使用20個(gè)ANCF繩索單元模擬，HHT中積分參數(shù)α取為-0.005，時(shí)間步長(zhǎng)取為10-4s。圖4是弦中點(diǎn)橫向位移的理論解(21)與ANCF數(shù)值解的對(duì)比(每30個(gè)時(shí)間步輸出一個(gè)結(jié)果)。弦中點(diǎn)的最大振幅為0.1201 mm，符合微幅振動(dòng)的假設(shè)。模擬結(jié)果與理論解符合得很好。

圖4 中點(diǎn)位移的對(duì)比

進(jìn)一步考察在弦中點(diǎn)受沖擊載荷作用下的響應(yīng)，載荷為0 s ～0.02 s內(nèi)一幅值Q為1 N或10 N的三角形力波。由于未有該情況下的解析解，故將ANCF的數(shù)值解與ANSYS中采用三次插值模式的Beam188單元的結(jié)果進(jìn)行對(duì)比。將弦均分為20個(gè)單元，ANCF采用α為-0.005得HHT格式，時(shí)間步長(zhǎng)為10-4s；ANSYS采用Newmark格式，時(shí)間步長(zhǎng)為10-5s。

圖5列出了在沖擊載荷幅值分別是1 N和 10 N 時(shí)，弦中點(diǎn)的位移時(shí)程，此時(shí)最大振幅分別超過(guò)了30 mm和65 mm，已不再符合微幅振動(dòng)假設(shè)。ANCF數(shù)值結(jié)果與ANSYS結(jié)果符合得較好。

5.2 雙桿擺和T字形桿擺

為考察本文算法對(duì)于大剛度梁的模擬和約束處理，模擬圖6的直徑為2 cm的圓截面鋼梁，彈性模量為200 GPa，密度為7800 kg/m3，其中圖6(a)為兩根長(zhǎng)為0.5 m的桿鉸接而形成的雙擺； 圖6(b)為T(mén)字形擺，各段長(zhǎng)度均為0.5 m。均由水平位置無(wú)初速度釋放，各段等分為兩個(gè)單元，采用α為-0.005的HHT方法，時(shí)間步長(zhǎng)均為10-3s。

圖5 沖擊載荷下中點(diǎn)位移的對(duì)比

圖7為鉸接的兩個(gè)節(jié)點(diǎn)距離的時(shí)間歷程，距離越小表明約束效果越好，圖7列出了不同罰因子C1水平的模擬結(jié)果。罰因子0.01EA/L(L=1 m)時(shí)，鉸接節(jié)點(diǎn)的最大距離為94.07 μm，約束效果好。增大罰因子能達(dá)到更好的約束效果，但罰因子增大到1013EA/L會(huì)導(dǎo)致不收斂。

T字?jǐn)[的模擬中，C1設(shè)置為103EA/L(L=1 m)，C2參照抗彎剛度EI/L進(jìn)行設(shè)置。圖8為節(jié)點(diǎn)切線(xiàn)夾角θ0與直角之差的時(shí)間歷程。罰因子0.1EI/L能達(dá)到較好的約束效果，θ0與直角相差不超過(guò)1.5°，進(jìn)一步增大罰因子至EI/L，θ0與直角相差不超過(guò)0.2°，但增大至100EI/L不收斂。

圖6 雙桿擺和T字形擺

圖7 鉸接節(jié)點(diǎn)距離的時(shí)間歷程

5.3 自由下落的柔繩

計(jì)算一根水平放置、左端鉸接于地面的柔繩，無(wú)初速度釋放，在重力作用下自由下落。繩長(zhǎng)為 1.2 m，圓截面面積為0.0018 m2，截面慣性矩為 1.215×10-8m4，密度為5540 kg/m3，彈性模量為 0.7 MPa。

將模擬結(jié)果與已有結(jié)果[18]進(jìn)行比較。節(jié)點(diǎn)等效力分別為小變形的Model I[19]、小/大變形的Model V[18]、符號(hào)運(yùn)算的Model S以及本文的公式(6～8)。此前的算例采用四階龍格-庫(kù)塔法的顯式格式，本文算例采用HHT法的隱式格式，時(shí)間步長(zhǎng)設(shè)置為10-3s。

圖8 θ0和90°的對(duì)比

圖9 T字形桿擺動(dòng)的ANCF結(jié)果和解析解

圖10是自由端鉛垂位置的時(shí)程曲線(xiàn)。本文算法的模擬結(jié)果與已有結(jié)果符合得較好，證明了本文公式及隱式迭代格式的正確性?？紤]繩索水平釋放，在0.7 s時(shí)解除約束的情況，圖11為柔繩各時(shí)刻的形狀。

圖10 自由端的鉛垂位置

圖11 繩索各時(shí)刻的形狀

6 結(jié) 論

(1) 為有效模擬大柔度梁/繩索結(jié)構(gòu)的大范圍運(yùn)動(dòng)和變形，針對(duì)ANCF大柔度梁/繩索單元建立了基于HHT方法的繩索隱式動(dòng)力學(xué)迭代格式。

(2) 將矢量叉乘改寫(xiě)為矩陣相乘，得到了 ANCF 大柔度梁/繩索單元的節(jié)點(diǎn)等效力公式，進(jìn)一步由節(jié)點(diǎn)等效力導(dǎo)出了切線(xiàn)剛度矩陣的全部公式。

(3) 本文算法的模擬結(jié)果與理論解及已有結(jié)果符合得較好，驗(yàn)證了節(jié)點(diǎn)等效力公式、切線(xiàn)剛度矩陣公式、及所建立的隱迭代格式的正確性。

(4) 采用罰方法處理ANCF的約束方程，實(shí)現(xiàn)了梁/繩索間的鉸接和剛性連接。數(shù)值模擬結(jié)果表明罰方法能夠有效實(shí)施此兩種約束。