羅 鑫，魏泳濤

(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ī)場等領(lǐng)域內(nèi)的橋式起重機(jī)(Overhead Cranes)[3,4]以及高速列車系統(tǒng)中的懸鏈線(Pantograph/Cantenary)[5]等。彈簧-阻尼-集中質(zhì)量模型[6-8]是數(shù)值模擬繩索的常用方法，但該法的繩索處理成分段折線而不再保持光滑曲線的形態(tài)，因而精度較低。有限元法在分析繩索時(shí)[9,10]，通常采用具有極低抗彎剛度或大柔度的梁單元以保證其形態(tài)的光滑性。在繩索運(yùn)動中，往往存在大范圍剛體運(yùn)動和大變形的耦合作用。但傳統(tǒng)梁單元的節(jié)點(diǎn)自由度僅反映變形的撓度和轉(zhuǎn)角，即單元的剛體運(yùn)動和變形是解耦的。浮動坐標(biāo)法[11]將繩索運(yùn)動考慮成隨浮動坐標(biāo)系的大范圍運(yùn)動與相對于浮動坐標(biāo)系的變形的疊加，但由此導(dǎo)出的科氏力和向心力等慣性力表達(dá)式十分復(fù)雜。

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

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

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

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

研究細(xì)長的等截面大柔度梁/繩索，在運(yùn)動/變形中始終維持光滑曲線，且橫截面是正多邊形，即橫截面對通過形心的任何軸都具有相同的慣性矩。本文基于文獻(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處的位置矢量和切線矢量排列而得的單元自由度向量。

(2)

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

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

(3)

(4)

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

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

(5)

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

(6)

(7)

(8)

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

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

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

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

(9)

(10)

(11)

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

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

動力學(xué)切線剛度矩陣(即迭代中的雅可比矩陣)為

(12)

(13)

(14)

(15)

(16)

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

(17)

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

4 鉸接和剛性連接的處理

ANCF大柔度梁/繩索單元使用切線的斜率為自由度，在剛性連接中，連接處的節(jié)點(diǎ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)的彈簧剛度。本文算例表明，可分別取為拉伸和彎曲切線剛度矩陣最大元素的1000倍。切線剛度矩陣第三項(xiàng)的推導(dǎo)中考慮矢量模的倒數(shù)為常數(shù)，從而簡化計(jì)算過程。

5 數(shù)值算例

5.1 弦的微幅和小幅振動

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

(21)

圖3 振動的弦

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

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

進(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)行對比。將弦均分為20個(gè)單元，ANCF采用α為-0.005得HHT格式，時(shí)間步長為10-4s；ANSYS采用Newmark格式，時(shí)間步長為10-5s。

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

5.2 雙桿擺和T字形桿擺

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

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

圖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會導(dǎo)致不收斂。

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

圖6 雙桿擺和T字形擺

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

5.3 自由下落的柔繩

計(jì)算一根水平放置、左端鉸接于地面的柔繩，無初速度釋放，在重力作用下自由下落。繩長為 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]、符號運(yùn)算的Model S以及本文的公式(6～8)。此前的算例采用四階龍格-庫塔法的顯式格式，本文算例采用HHT法的隱式格式，時(shí)間步長設(shè)置為10-3s。

圖8 θ0和90°的對比

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

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

圖10 自由端的鉛垂位置

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

6 結(jié) 論

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

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

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

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