一種新開發的ANCF纜索單元

鄭茂盛, 陳建平, 童明波

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

線纜是一類具備大撓度變形特征的柔性結構,可用于控制安裝、設備連接、電力傳輸、光信號傳輸等。由于重量輕、阻尼小、抗拉強度高、柔性好等特點使其在工程實踐中被廣泛應用[1-2]。在現代工藝中幾乎可以制造出滿足各類工程需要的線纜。與此同時,復雜的層級結構、多樣的繩股纏繞方式、不同的層間材料屬性都極大地影響著線纜的動力學行為。就系繩操作的水下系統而言,系留線纜所帶來的不必要的阻力和非線性內力很大程度上影響著水下機器人運動。而對于空中系浮系統而言,無論是在升空階段、系泊階段還是回收階段,系留線纜始終起到承力的關鍵作用[3]。因此,對于該類纜索結構的動力學分析至關重要,精確有效的建模方案也引起了很多研究人員的關注。

多層圓形截面的纜索在材料屬性和幾何變形中都表現出了典型的非線性特征。對于這類具有幾何大變形特征的柔性結構進行動力學建模主要方法有以下幾種方法。

方法一是利用集中質量法對纜索進行大規模離散,相鄰質點之間利用彈簧-扭簧進行連接,并取各質點在慣性空間中的直角坐標為廣義坐標,利用 Lagrange 方程建立系統的運動微分方程,Buckham等[4]對此進行了充分的研究。這種方法簡便易懂,但其求解精度較低,且當節點劃分過密時,求解效率將大大降低。

方法二是由Kamman等[5]提出的利用基于偽剛體理論的有限段離散化方法進行動力學建模,但該方法的缺點同樣在于隨著有限段數目的增加,動力學方程的求解效率會大大降低。

方法三是可以利用連續介質力學對其進行動力學建模。將纜索處理為具有大撓度變形的細長彈性梁,再利用牛頓定律或Hamilton 原理推導出運動的偏微分方程,最后利用差分格式進行離散求解。Koh等[6]對此進行了詳細的敘述,但其難點在于非線性偏微分方程的數值求解。

方法四是由Shabana[7]于20世紀90年代提出的絕對節點坐標方法(absolute nodal coordinate formulation,ANCF)。其基本思想是將有限元中的單元節點用慣性坐標系中的絕對坐標進行描述,并用斜率矢量代替傳統的節點轉角坐標以得到系統在全局坐標系下的運動方程。該方法具有質量矩陣為常值矩陣、無科氏力和離心力等優點。其能夠很好地反映出柔性纜在拉伸、扭轉、彎曲等變形下的復雜動力學行為[8]。

與Gerstmayr等[9]提出的高階ANCF梁單元不同,現利用Shen等[10]的全參數高階ANCF梁單元的形狀函數對纜索進行運動學描述。原因在于后者在抑制泊松鎖定方面具有極大優勢。此外,在對單元子域數值積分的研究工作中,Orzechowski[11]選用高斯積分點對圓形截面單元進行了數值分析。隨后,Orzechowski等[12]又開發了一種能夠分析矩形橫截面的層合梁模型。此后,Lan等[13]對多層圓形截面單元子域的積分提出了合理、有效的數值策略。

針對于一些具有絕緣層的工程纜索而言,護套層通常選用橡膠、聚酯氨等材料以用于絕緣、抗磨和防腐化等。而這些材料的泊松比接近于0.5,能夠表現出不可壓縮和超彈性的力學特性。對于該類不可壓縮材料的研究工作在文獻[14]中得到概括總結,其中包括Neo-hookean模型、Mooney-Rivlin模型以及Yeoh模型[14]等。而由于應變能密度函數的展開式僅保留了有限的項數,Neo-hookean模型、Mooney-Rivlin 模型均無法準確地反映大變形狀態下材料的應力-應變關系,特別是無法反映“S”形應力-應變曲線,因此其只適用于應變水平較低的情況[14]。相反,雖然Yeoh 模型[14]在預測大變形簡單應變問題上具有一定優勢,但其在預測等雙軸拉伸問題上將會出現“偏軟”的較大誤差,無法精確地處理大變形復雜應變問題。

對于利用絕對節點坐標方法對上述工程纜索進行動力學建模的研究而言,近些年學者多將超彈性材料的Mooney-Rivlin模型[15-17]、Yeoh模型[14]應用到ANCF單元中。徐齊平等[18]通過引入罰函數,推導Yeoh模型[14]的彈性力矩陣,并采用絕對節點坐標方法建立大變形硅膠板的動力學方程。郭嘉輝等[19]借助平面二維ANCF梁單元分別構建了基于Neo-hookean模型、Mooney-Rivlin[15-17]模型以及Yeoh模型[14]的不可壓縮超彈性橡膠梁的動力學模型。

但是對于廣泛應用于浮空器系泊當中的系留纜索而言,應用ANCF方法對其精確動力學建模的相關研究較少。此外,為了進一步克服Mooney-Rivlin模型[15-17]難以準確地反映大變形狀態下超彈性材料的應力-應變關系以及應變適用水平低、范圍小等問題,首次將改進后的Yeoh模型[14]應用于系留纜索的動力學建模當中,并利用三維高階ANCF梁單元來精確刻畫纜索的運動及變形[20],開發一種新型的基于改進Yeoh超彈性模型的ANCF纜索單元。同時,通過比較幾種應用不同超彈性本構模型的ANCF單元,對其動力學響應的特征與差異進行對比分析,為后續對于徑向多層圓截面纜索的動力學建模研究提供一定的指導。

1 材料的力學特性

圖1為車載式系留監測系統的示意圖。針對于系留系統中的系留線纜而言,其徑向層級結構由內到外依次為內導體、絕緣層、承力層、外護套。銅、鋁等高導電的金屬材料多用于內導體;絕緣層多選用絕緣效果不錯的聚乙烯(PE)、交聯聚乙烯(XLPE)以及氟塑料等;芳綸、聚對苯撐苯并二噁唑纖維(PBO)等新型高強度材料往往作為抗拉件應用于承力層;橡膠、熱可塑性聚氨酯(TPU)材料則應用于護套層當中[3]。線纜的橫截面如圖2所示。

圖1 車載式系留檢測系統

圖2 線纜的橫截面

1.1 黏彈性材料的Kelvin-Voigt模型

García-Vallejo等[16]基于瑞利函數提出ANCF單元的內部阻尼模型。其在考慮泊松效應的同時可以很好地解釋多軸應力狀態下的線性黏彈性關系。對于均勻各向同性黏彈性材料而言,其本構方程的表達式為

(1)

(2)

D=

(3)

對黏彈性材料而言,有

(4)

體積模量K可表示為

(5)

式中：λ和μ分別為第一拉梅常數和第二拉梅常數;E為楊氏模量;v為泊松比;γd和γs分別為偏應力耗散因子和膨脹應力耗散因子。

1.2 改進后的Yeoh模型

引入右Cauchy-Green變形張量C=FTF,其中F為變形梯度。于是得到右Cauchy-Green變形張量的3個不變量I1、I2、I3,可分別表示為

(6)

式(6)中：tr為矩陣的跡。

可以將改進后的Yeoh模型[14]中的應變能密度函數寫為

Wp=C10(IC-3)+C20(IC-3)2+
C30(IC-3)3+C01(IIC-3)

(7)

式(7)中：IC和IIC分別為變形張量C的第一、第二不變量;C10、C20、C30、C01為4個模型常數。

第二Piola-Kirchhoff張量可表示為

(8)

式(8)中：ε為Green-Lagrange應變張量。

2 動力學建模

2.1 運動學分析

圖3為三維全參數高階ANCF梁單元的示意圖。

ri為節點i處的位置矢量;ri,x,ri,y,ri,z為節點i處的3個梯度矢量;ri,yz,ri,yy,ri,zz在幾何上解釋為曲率矢量,其對泊松鎖定的抑制具有關鍵作用

利用如式(9)所示的多項式來近似表示單元任意點于慣性坐標系下的位移。

r=[a0+a1x+a2y+a3z+a4xy+a5xz+a6x2+a7x3+a8yz+a9y2+a10z2+a11xyz+a12xy2+a13xz2b0+b1x+b2y+b3z+b4xy+b5xz+b6x2+b7x3+b8yz+b9y2+b10z2+b11xyz+b12xy2+b13xz2c0+c1x+c2y+c3z+c4xy+c5xz+c6x2+c7x3+c8yz+c9y2+c10z2+c11xyz+c12xy2+c13xz2]

(9)

此時利用節點的全局位置矢量、3個方向的梯度矢量以及曲率矢量作為單元的廣義坐標,可將式(9)變換為

(10)

式(10)中：I為單位矩陣;q為單元的廣義坐標列陣。

s1=1-3ξ2+2ξ3

s2=l(ξ-2ξ2+ξ3)

s3=lη(1-ξ)

s4=lζ(1-ξ)

s5=l2ηζ(1-ξ)

s8=3ξ2-2ξ3

s9=l(-ξ2+ξ3)

s10=lξη

s11=lξζ

s12=l2ξηζ

(11)

單元內任意點的速度與加速度坐標列陣同樣可表示為

(12)

2.2 動力學分析

2.2.1 慣性力的虛功

基于2.1節中對ANCF單元運動學的分析,其慣性力的虛功可表示為

(13)

式(13)中：ρ1、ρ2分別為內、外層的材料密度;A1、A2分別為內、外層的截面面積;VI、VII分別為內、外層的體積積分域。

結合式(10)、式(12),可將式(13)變換為

(14)

式(14)中：δW1為慣性力的虛功;M為纜索單元的質量矩陣,是一個42×42階的對稱正定常值矩陣;δq為廣義坐標的變分。

2.2.2 黏彈性力的虛功

基于連續介質力學的相關理論,6個應變分量于全局坐標系下分量值可寫為

(15)

式(15)中：rx、ry、rz為位置矢量r相對于x、y、z坐標方向的3個梯度矢量。

基于1.1節中對黏彈性材料本構關系的分析,此時黏彈性力的虛功可表示為

δW2=-δqTQLayerI

(16)

式(16)中：δW2為黏彈性力的虛功;QLayerI為廣義非線性黏彈性力列陣,可以記為

(17)

式(17)中：K為42×1階的廣義彈性力列陣,其具備高度的非線性特征;C為42×42階的廣義阻尼系數矩陣。

K與C也可以記為

(18)

2.2.3 不可壓縮材料彈性力的虛功

基于1.2節中對護套層力學性能的基本分析,可以得到護套層非線性彈性力的虛功表達式為

δW3=-δqTQLayerII

(19)

式(19)中：δW3為護套層彈性力的虛功;QLayerII為其非線性彈性力列陣,QLayerII的表達式為

(20)

而改進后的Yeoh模型常數取文獻[14]中的擬合結果為：C10=1.7×105Pa,C20=-1.55×103Pa,C30=46.1 Pa,C01=5.24×103Pa。

2.2.4 外力的虛功

若纜索的內、外層單元分別受到均布體力f1、f2和作用于a位置的集中力p影響時,其由外力項所引起的虛功可表示為

(21)

式(21)中：ra為a位置于全局坐標系下的位置矢量;δW4為外力的虛功。將式(21)可簡化為

δW4=δqTQ

(22)

Q為單元的廣義外力列陣,可表示為

(23)

式(23)中：Sa為位置a處的形函數矩陣S。

2.2.5 動力學方程

基于虛功原理,綜合上述對相關虛功的具體分析,建立ANCF纜索單元的動力學方程。

根據虛功原理有

δW1+δW2+δW3+δW4=0

(24)

由于廣義坐標變分的獨立性,將式(14)、式(16)、式(19)、式(22)代入式(24)中,可以得到單元的動力學方程為

(25)

對于n-1個ANCF單元而言,其共有n個不同的特征節點。結合式(25)可以得到整體的動力學方程為

(26)

式(26)中：Ci為第i個單元的阻尼系數矩陣;Ki為內層的線彈性力矩陣;QLayerIIi為外層的彈性力矩陣;Qi為外力矩陣;qi為廣義坐標矩陣;i=1,2,…,n。

在此基礎上,利用布爾邏輯型矩陣B可以得到如式(27)所示的轉換關系。

(27)

式(27)中：w為各個節點位置處的廣義坐標所組成的列陣,其不包含重復節點。

于是式(26)可進一步改寫為

(28)

(29)

3 數值方法

3.1 子域積分的數值方法

對于第2節中所涉及的體積分問題實質上很難得到一個準確的解析解。在Matlab程序調試當中利用Syms符號函數進行積分求解會極大地降低工作效率。為此,采用Gauss-Legendre數值積分策略,先利用坐標變換將空間圓柱型積分域變換成嚴格的高斯型積分域,再取用5個×5個×5個高斯積分點精確的求解相應的空間體積分問題。5個×5個高斯積分點在平面上的位置排布如圖4所示。

圖4 5×5高斯積分點在平面上的分布

內芯層體積分可以改寫成式(30)所示的形式。

(30)

式(30)中：F為以笛卡爾坐標為自變量的待積函數;r1為內芯層圓截面的半徑;r為纜索單元的半徑;Φ為積分的轉角變量;L為纜索單元的長度;Ai、Bj、Ck為對應序號的高斯積分點;wi、wj、wk為對應積分點的權重系數。

同理,護套層體積分可以改寫成

(31)

式(31)中：r2為纜索單元截面的半徑。

3.2 動力學求解的數值方法

為保證微分方程數值求解的穩定性,避免出現由于方程剛性較強而產生的失穩現象,利用隱式方法進行數值求解,這是因為相較于顯示方法而言,隱式法的穩定域通常更大,不必選用過分小的步長來滿足苛刻的穩定性條件。但是任何隱式方法均涉及非線性方程組解的數值迭代,其實現過程復雜、計算效率低下。為此,設計一種具有隱式效果的差分格式,以便對上述動力學微分方程組進行高效準確的數值求解。

將式(28)進行等價性變換得

(32)

式(32)中：υ為各個節點位置處的廣義速度所組成的列陣,其不包含重復節點。

一階隱式歐拉法的格式可以寫為

(33)

式(33)中：wi、vi分別為ti時刻下的位移變量與速度變量;h為固定的時間步長。

(34)

式(34)保留了部分隱式方法的基本特征,但卻是一種顯示差分格式,其在擴大穩定域的同時改善了原來隱式方法求解效率低的特征。

3.3 靜力學求解的數值方法

為提高非線性方程組數值迭代的收斂速度,因此采用Newton-Raphson方法進行數值迭代。將式(28)進行等價性變換得

(35)

其迭代格式可以寫為

(36)

4 數值算例

4.1 靜力學數值仿真

為了將相關理論解同數值解進行對比,選用靜態小變形問題進行說明。將水平力與垂直力分別施加在懸臂纜的自由端,并分析自由端位移的變化情況。其中,水平力Fh與垂直力Fv的大小均為10 N。受力情況如圖5所示。模型的基本參數如表1所示。

表1 仿真參數Table 1 Parameters of simulation

圖5 兩個小變形問題

在水平力Fh的作用下,末端伸長量Δl在小變形靜力學問題下的理論解為

(37)

式(37)中：E1、E2分別為內、外層材料的楊氏模量;Δl=1.586 6×10-7m。

利用Newton-Raphson迭代的結果如表2所示。

表2 水平力下的數值解Table 2 Numerical solution under horizontal force

在垂直力Fv的作用下,撓度w在小變形靜力學問題下的理論解為

(38)

其具體的數值結果為w=2.100 7×10-3m。利用Newton-Raphson迭代的結果如表3所示。

表3 垂直力下的數值解Table 3 Numerical solution under vertical force

由表2、表3可知,ANCF單元數目的增加使得仿真解更加貼近理論解。自由端位移的精度與收斂性得到了一定的保證。

ANCF單元的數目為5個時,撓度迭代的基本情況如圖6所示。而對于應用Mooney-Rivlin模型與原始的Yeoh模型的ANCF單元而言,迭代過程中數值解的變化規律如圖7所示。

圖6 垂直力下撓度解的迭代

圖7 兩種模型下數值解的迭代

由圖6、圖7可知,應用Mooney-Rivlin的單元迭代收斂的速度較慢,而對于應用另外兩種模型的單元而言,數值迭代均能夠快速、高精度地收斂到正確解。

4.2 動力學數值仿真

4.2.1 物理擺

水平放置的柔性擺在重力的作用下下落,研究系統的動力學響應特征以及能量的變化規律。柔性擺的基本參數如表4所示。偏應力耗散因子與膨脹應力耗散因子分別取0.01與0.1。

表4 纜索的基本參數Table 4 Parameters of cable

圖8為擺的自由端沿著x、y方向的位移變化曲線。圖9為擺落過程中動能、重力勢能、應變能以及耗散能隨時間的變化規律。

圖8 自由端的位移

圖9 能量的變化情況

由圖9可知,能量的變化規律符合預期的結果。動能與重力勢能呈現出反相的基本特征,即系統的動能取極大值時,重力勢能取極小值。這是由于柔擺下落過程中勢能主要向動能轉化,僅有一小部分經過阻尼力耗散后的能量轉化成了單元內部的應變能?？傮w來說,系統的總能量在整個模擬時間內保持不變,單元具備能量守恒的基本特征。這也證明了動力學數值方法的可行性與單元動力學響應的可靠性。

4.2.2 懸臂纜的擺落

主要探究應用幾種不同的不可壓縮材料模型對單元數值響應的影響,這是因為不同的超彈性解釋模型在大變形、強非線性特征下的力學特性各有差異。表4給出了纜索的基本參數信息。

圖10為應用3種不可壓縮材料模型的ANCF單元于不同時刻的位形曲線?？梢钥闯?應用改進后Yeoh模型的單元于不同時刻下的位形曲線與原始Yeoh模型結果高度一致。這意味著原始Yeoh模型和改進后Yeoh模型對單元動力學行為的影響效果基本一致。而對于應用Mooney-Rivlin模型的ANCF單元而言,其自由端的位移隨著時間的不斷推移明顯超前于應用另外兩個模型的纜索單元。其勢能轉化為動能的速度更快。在總能量守恒的前提下,Mooney-Rivlin模型所產生的應變能低于原始Yeoh模型和改進后Yeoh模型,導致其無法準確地反映柔性纜在大尺度變形下的運動情況。

圖10 應用3種不可壓縮模型的單元的位置構型

5 結論

提出一種適用于線纜動力學建模的ANCF纜索單元。其中利用Kelvin-Voigt模型來解釋內層結構的黏彈性性質,而將改進的Yeoh模型應用于護套層單元當中。在此基礎上,通過一系列靜力學與動力學數值仿真得到如下結論。

(1)單元能夠高精度快速地收斂到理論解,并且對比了應用3種不可壓縮模型的單元的收斂特性,證實了本單元具備較好的收斂性能。

(2)通過柔性擺的數值模擬證實了單元具備能量守恒的基本特征,證明了數值方法的可行性與單元動力學行為的可靠性。

(3)應用改進的Yeoh模型的單元能夠很好地反映柔性纜在大尺度變形下的運動。而相較于應用Mooney-Rivlin模型的ANCF單元而言,其無法準確地反映柔性纜在大變形條件下的運動特征,進一步證實了所提出的ANCF纜索單元的合理性與優越性。