基于絕對節點坐標法的鋼索結構動力學模型與特性

錢彥懿　王　慧　余海東

1.上海交通大學上海市復雜薄板結構數字化制造重點實驗室，上海，200240 2.航空工業第一飛機設計研究院，西安，710089

0　引言

鋼索傳動具有結構簡單、重量輕、抗拉強度高的特點，被廣泛應用于航空飛機的軟式操縱系統。隨著大型及超大型飛機的出現，對鋼索在長距離傳動過程中的變形及響應遲滯時間提出了更高的要求。由于鋼索具有大柔性、大變形、大位移等顯著的幾何非線性特點，加之材料本身的物理非線性特征，所以鋼索在運動變形中會發生剛度變化；且鋼索通常由多股鋼絲捻制而成，加載后由于股內鋼絲間的相互作用也會產生截面內的扭轉等變形。傳統的有限元模型對繩索結構力學行為的精確描述存在一定的局限，因此需提出新的方法來描述鋼索受力過程中的剛度變化及截面變形，實現鋼索結構變形的精確預測。

對于柔性繩索的建模，通常分為解析法［1-2］和有限元法兩類。隨著計算機技術的發展，有限元法在繩索結構非線性分析領域的應用日漸廣泛。目前較為成熟的索單元模型有兩節點直桿單元、兩節點曲線索單元和多節點等參元。直桿單元［3］形式簡單但精度較差，僅適用于位移、變形很大的空間繩索。曲線索單元不適用于張緊結構，且達到較高求解精度所需單元數多，計算量大［4］。多節點等參元能較好地模擬繩索真實構型且精度較高［5］，但由于模型節點多，自由度增加顯著，對大跨多索結構難以處理，在小垂度拉索結構中不能明顯提高計算精度，反而大大增加了計算量。以上模型均未考慮繩索運動過程中變剛度及截面變形對整體結構變形的影響，且無法兼顧求解效率和計算精度問題。

SHABANA［6］提出絕對節點坐標法（absolute nodal coordinate formulation，ANCF）來描述大柔性體的運動變形，其理論基礎是連續介質力學和有限元法［7］。絕對節點坐標法突破了傳統有限元方法的截面不變形假設，可描述單元體的剪切、扭轉、彎曲及拉伸等各類變形，對求解大范圍轉動、大變形問題具有明顯優勢［8-9］，且適用于具有變剛度特性的柔性繩索類物體建模。

本文基于絕對節點坐標法，考慮鋼索的大柔性特點及其在加載過程中出現的截面內旋轉變形，引入高階位移插值函數描述鋼索截面變形和旋轉，推導了鋼索單元的質量矩陣、廣義外力矩陣和剛度矩陣，并對不同載荷條件、結構參數以及材料特性的鋼索進行了靜力學變形和動力學響應遲滯仿真分析，仿真結果能夠更真實地反映鋼索實際變形狀況。

1　基于絕對節點坐標法的鋼索單元建模

1.1　鋼索單元變形描述

以三維鋼索單元為對象，描述其在空間內的拉伸、彎曲和剪切變形。圖1為變形后的基于絕對節點坐標法鋼索單元（下稱ANCF鋼索單元）模型示意圖，其中，每個鋼索單元包含I、J兩個節點，分別位于端點處。鋼索單元內任意一點P的位置坐標列陣r及單元節點坐標列陣e均在全局坐標系XYZ中描述。依據絕對節點坐標法，三維鋼索單元內任意一點P在全局坐標系下的位置坐標列陣r可通過對節點I、J的插值來求解：

式中，x為變形前點P在單元坐標系xyz中的位置坐標矢量；S為單元形函數矩陣；r1、r2、r3分別為點P在全局坐標系下沿X軸、Y軸和Z軸方向的投影分量。

圖1　變形后的ANCF鋼索單元模型Fig.1　Geometric definition of the deformed ANCF cable element model

對于三維變形鋼索單元，其位移場可用單元坐標x、y、z的三次插值多項式表達：

式中，ai、bi、ci為插值多項式的待定系數。

式（2）中包含24個待定系數，對于兩節點鋼索單元（I，J），每個節點應設置12個廣義坐標分量（即邊界條件），以確定式中的待定系數。每個單元節點坐標列陣e由其節點位置坐標列陣和對x、y、z分別求微分的斜率矢量來描述，共24個廣義坐標形成單元節點坐標列陣：

式中，eI、eJ分別為節點 I( 0 ,0)和J( l ,0)的廣義坐標，l為單元初始長度；rI、rJ分別為節點I和節點J在全局坐標系下的位置坐標；?rI/?x、?rJ/?x分別為節點I和節點J的斜率坐標，表示鋼索截面在全局坐標系中的方向。

將單元的兩個節點位置坐標eI和eJ代入式（2），可求得式（2）中的24個待定系數，它們均為單元廣義坐標的函數，經整理即可得單元的形函數：

式中，I為3×3的單位矩陣。

1.2　鋼索單元質量矩陣

將式（1）對時間求導，可得到鋼索軸向任意一點的絕對速度

鋼索結構的動能

式中，ρ為材料密度；V為鋼索單元體積；M為鋼索單元質量矩陣。

由式（6）可知，質量矩陣為常量矩陣，對于等截面鋼索，其單元質量矩陣可由沿x方向積分得到：

式中，A為鋼索橫截面積。

1.3　鋼索單元廣義外力矩陣

鋼索系統所受廣義外力包括集中力、力矩及分布力（如重力等），可依據虛功原理進行求解。設F=[FxFyFz]T為鋼索單元上任意一點所受的外力，則此外力所做的虛功

對于大跨度鋼索結構，應考慮鋼索自重對其變形的影響，將重力考慮為分布力，其虛功可通過對單元體積求積分得到，進而求得單元廣義重力矩陣：

式中,g為重力加速度。

1.4　考慮阻尼影響的鋼索單元剛度矩陣

鋼索單元的廣義彈性力可由連續介質力學中非線性應變位移關系推導得到，鋼索單元的應變張量可用鋼索單元變形梯度J描述。對于鋼索單元模型，考慮鋼索在預緊力下運動及變形特點，其變形梯度

式中，r、r0分別為變形前后任意一點在全局坐標系中的坐標；J0為常量矩陣，與單元的初始構型有關。

當選取單元初始構型為直線型，且單元坐標系平行于全局坐標系時，J0為單位矩陣，因此式（11）可簡化為

考慮鋼索結構的大變形問題，利用單元變形梯度，由彈性力引起的鋼索單元應變可由格林拉格朗日應變張量定義：

εa為對稱張量，可將其寫為由6個獨立應變分量構成的列陣ε：

對于本文軸向加載的鋼索模型，考慮到鋼索截面尺寸遠小于其軸向尺寸，因此沿鋼索軸向的拉伸應力為其廣義彈性力主要組成部分。同時由于傳動鋼索在加載之前處于張緊狀態，因此初始張力引起的結構廣義彈性力也不容忽視。故包含鋼索初始張力及系統阻尼的單元軸向應變ε?1可表示為

式中，ε1、ε?1分別為廣義彈性力引起的軸向應變和應變率；F0為鋼索初始張力；E為材料彈性模量；C為材料阻尼系數。

根據材料本構模型，得到單元應力

根據虛功原理，可得由廣義彈性力所做的虛功

2　鋼索變形仿真計算分析

鋼索受力模型見圖2。鋼索水平放置，一端固定，另一端受軸向載荷。初始時刻鋼索在初始張力的作用下保持張緊狀態，在末端軸向載荷的作用下產生變形。由于重力方向垂直于鋼索傳動方向，故鋼索自重產生的撓度對其運動變形的影響不可忽略，為確保仿真精度，應考慮鋼索自重的影響，因此，鋼索結構在廣義彈性力Qe、廣義重力Qg及由末端載荷T引起的廣義外力Qf的作用下保持平衡。由式（9）、式（10）及式（16）可得鋼索的受力平衡方程

圖2　鋼索受力模型Fig.2　The cable with an axial force

仿真中將鋼索劃分為10個單元，整個系統包含11個節點，每個節點包含12個節點坐標，固定端點處存在移動約束，因此系統共有129個節點坐標。

鋼索材料密度ρ=6 140 kg/m3，彈性模量E=45 GPa，泊松比ν=0.3，重力加速度g=9.81 m s2，鋼索截面為圓形，其長度通過改變單元長度來控制。對鋼索的變形進行仿真計算，將MATLAB仿真結果與鋼索拉伸實驗結果進行比較，驗證模型對鋼索變形求解的正確性，再進一步研究初始張力和鋼索結構參數對其變形的影響。

2.1　仿真與實驗結果比較

以原長600 mm、直徑3.2 mm的國產碳鋼鋼索為研究對象，初始張力F0=200 N，軸向載荷T從50 N增至1 800 N，用MATLAB計算得到加載過程中各時刻的末端軸向位移，以此作為鋼索的變形量u，并與對應時刻拉伸實驗的測量結果作比較。結果表明：仿真計算與實驗測量結果間的誤差相對于鋼索原長的百分比在0.001 2%～0.021 4%之間，仿真結果與實驗結果非常接近。圖3為鋼索的載荷變形曲線，可以看出，隨著軸向載荷的增加，基于ANCF鋼索單元計算得到的鋼索變形量的變化趨勢顯現出一定程度的非線性，這主要是因為鋼索的旋轉運動消耗了小部分應變能；數值計算結果和實際拉伸實驗測量結果基本吻合，從而驗證了基于ANCF鋼索單元的鋼索靜力學變形模型具備較高的求解精度。

圖3　仿真及實驗結果的載荷　變形曲線Fig.3　Load-displacement curves for numerical and experimental results

2.2　初始張力對鋼索變形的影響

保持其余參數不變，初始張力F0由50 N增至1 500 N，分析初始張力變化對鋼索變形的影響。其中，鋼索直徑d=3.2 mm，軸向載荷T=500 N，鋼索長度L為1 m、3 m、6 m及10 m，分別計算得到鋼索變形量u。考慮到鋼索的變形量與其自身長度線性相關，進行量綱一化，采用單位長度變形量（即變形率α）作為衡量不同長度鋼索變形程度的指標。不同長度的鋼索在不同初始張力下的變形率見表1。由表1可知，相同長度鋼索的變形率隨初始張力的增大呈減小趨勢，但減小的幅度極小，幾乎可忽略不計，其主要原因是初始張力的增大在一定程度上導致鋼索結構剛度的增大，即鋼索結構抵抗變形的能力增強，從而在相同載荷下變形量減小，變形率減小；但本文所討論的是鋼索張緊后，再對其加載產生的變形與初始張力的關系，由于加載前鋼索已處于張緊狀態，且實際上使其張緊所需初始張力很小，因而在本文所取初始張力范圍內的鋼索結構剛度變化較小，故其變形率的減小幅度較小。通過對比相同初始張力、不同長度鋼索的變形率可知，變形率與鋼索長度無關。

表1　不同初始張力及鋼索長度下的變形率Tab.1　Deformation rate of cables with various preloads and lengths

圖4所示為4種長度的鋼索在不同初始張力下的變形量，可以看出，4條曲線均呈微幅線性下降趨勢，即鋼索變形量隨初始張力的增大而緩慢減小。

圖4　不同初始張力及鋼索長度下的鋼索變形量Fig.4　Deformation of cables with various preloads and lengths

2.3　鋼索幾何參數對其變形的影響

2.3.1不同直徑鋼索的變形

為得到鋼索直徑對其變形的影響規律，保持其余參數不變，分別計算直徑為1.6 mm、3.2 mm及6.3 mm的3種鋼索在不同軸向載荷下的變形。其中，鋼索長度均為1 m，初始張力F0=500 N，軸向載荷T由50 N增至2 000 N。

圖5為不同直徑鋼索的加載變形曲線，可以看出，不同直徑鋼索的變形量均隨載荷的增大呈線性增大趨勢，且鋼索直徑越小，鋼索變形量的增大速率也越快（即圖中曲線的斜率越大）。這表明小直徑鋼索對載荷變化的敏感程度要高于大直徑鋼索對載荷變化的敏感程度。為了更直觀地反映鋼索直徑對其變形難易程度的影響，保持軸向載荷不變（均為500 N），比較不同直徑的鋼索變形量。由圖6可以看出，當加載條件相同時，鋼索的變形量隨直徑的增大，呈非線性減小趨勢。這是由于鋼索直徑的增大導致其結構剛度增大，因此在相同加載條件下，鋼索變形量減小。

圖5　不同直徑鋼索的載荷　變形曲線Fig.5　Load-deformation curves of cables with various diameters

圖6　不同直徑鋼索的變形量Fig.6　Deformation of cables with different diameters

2.3.2不同長度鋼索的變形

為研究鋼索長度對其變形的影響，保持其余參數不變，分別計算原長1～15 m的鋼索在相同加載條件下的變形量。其中，鋼索直徑均為3.2 mm，初始張力F0=1 000 N，軸向載荷T=1 000 N。圖7所示為不同長度鋼索的變形量，可以看出，在相同加載條件下，鋼索的變形量隨著鋼索長度的增大呈線性增大趨勢，而變形率保持不變。這表明鋼索的變形難易程度與其自身長度無關，不同長度的鋼索在同等條件下的變形量與其原長成正比，但其變形量相對于原長的百分比不變。

圖7　不同長度鋼索的變形量Fig.7　Deformation of cables with various lengths

3　鋼索動力學仿真計算分析

本文對鋼索動力學特性的研究主要著眼于鋼索傳動過程中在階躍載荷作用下的響應遲滯效應。初始時刻鋼索在200 N的初始張力作用下保持水平張緊，在末端軸向激勵的作用下產生運動變形。在整個仿真過程的前0.5 ms內，鋼索一端受到F=500 N的階躍載荷，另一端被固定，鋼索各個部分從加載的一端開始依次張緊，并獲得位移和速度。對于大跨度柔性鋼索結構，由于鋼索自重產生的撓度對其動力學特性的影響不可忽略，為確保仿真精度，應考慮鋼索自重的影響。由前面推導的M、Qe、Qg、Qf建立基于牛頓方程的鋼索系統動力學方程：式（18）為一組二階常微分方程，其求解可通過四階龍格庫塔方法實現。仿真中鋼索分為10個單元，整個系統包含11個節點，每個節點包含12個節點坐標，固定端點處存在移動約束，因此系統共有129個節點坐標。

圖8　節點軸向位移Fig.8　Displacement of node in X direction

在仿真計算中，采用位移從加載端節點傳遞到前端第2個節點的時間td來衡量遲滯效應的程度。圖8為仿真過程中該節點軸向位移隨時間的變化曲線，可以看出，由于遲滯效應的存在，最前端節點在初始的一段時間（約0.2 ms）的位移幾乎為零，之后位移急劇增大至接近其靜力拉伸變形量后保持不變，因此，遲滯時間td可定義為從初始時刻至前端第2個節點位移開始大幅增大時刻所經歷的時長。

依據應力波的傳播理論，應變和質點速度可看作以波的形式沿繩索傳播，其傳播速度v可依據鋼索物質特性計算求得：

從而可計算得到理論遲滯時間

仿真過程中，彈性模量E=45 GPa，泊松比ν=0.3，重力加速度g=9.81 m s2，鋼索截面為圓形，其長度通過改變單元長度來控制。對于鋼索的動力學仿真，主要研究鋼索結構參數及材料特性對其響應遲滯的影響，并將MATLAB計算得到的遲滯時間td與解析解t0d進行比較，以驗證模型對鋼索動態遲滯預測的合理性。

3.1　鋼索結構參數對其遲滯的影響

鋼索結構參數（包括鋼索直徑及鋼索長度）會對其動態響應遲滯產生影響。保持其余參數不變，分別對3種直徑（1.6 mm、3.2 mm和6.3 mm）及8種不同長度（0.6～15 m）的鋼索進行動力學仿真計算，由動態特性曲線得到鋼索響應遲滯時間。其中，鋼索材料密度ρ=6 140 kg/m3，初始張力F0=200 N，軸向階躍載荷為500 N。

圖9和圖10分別為不同直徑及不同長度鋼索的動力學仿真遲滯時間與理論遲滯時間曲線。由圖9可以看出，鋼索直徑對其遲滯效應幾乎沒有影響，不同直徑鋼索的遲滯時間相同。仿真計算結果變化趨勢與理論值相符，誤差百分比約4.6%。由圖10可以看出，鋼索長度對遲滯時間的影響極大，遲滯時間隨鋼索長度的增大，呈線性增大趨勢，符合理論變化曲線，且誤差百分比不超過3.8%。由應力波理論可知，在鋼索材料屬性不變而長度增大的情況下，由于應力波的傳播速度不變而傳播距離增加，因而傳播時間增加，導致遲滯時間增加。

圖9　遲滯時間隨鋼索直徑變化的規律Fig.9　Delay time of cables with various diameters

圖10　遲滯時間隨鋼索長度變化的規律Fig.10　Delay time of cables with various lengths

圖11　遲滯時間隨鋼索密度變化的規律Fig.11　Delay time of cables with various mass densities

圖12　遲滯時間隨彈性模量變化的規律Fig.12　Delay time of cables with various elasticity modulus

3.2　鋼索材料特性對其遲滯的影響

鋼索材料屬性影響鋼索結構剛度，因而對其動態特性存在影響。作為初步研究，鋼索材料密度分別取 2 660 kg/m3、4 550 kg/m3、6 140 kg/m3；彈性模量分別取 0.45 GPa、4.5 GPa、45 GPa。仿真中鋼索直徑均為3.2 mm，鋼索長度為1 m，初始張力為200 N，軸向階躍載荷為500 N。圖11和圖12分別為不同材料密度及不同彈性模量鋼索的遲滯時間變化曲線。由圖11和圖12可以看出，仿真得到的不同材料密度及彈性模量的鋼索動態響應遲滯時間曲線均與其對應的理論遲滯時間曲線相吻合，最大誤差百分比分別為4.6%和10%。鋼索密度及彈性模量均對其動態響應有較大影響，遲滯時間隨著鋼索密度的增大，呈近似線性增大趨勢，密度越大遲滯時間增加越緩慢；彈性模量對遲滯時間的影響趨勢剛好相反，隨著彈性模量的增大，遲滯時間呈非線性減少趨勢，并且減少速率逐漸趨緩。依據應力波傳播理論和式（19）可知，應力波的傳播速度由材料密度及彈性模量決定，彈性模量增大或材料密度減小導致傳播速度增大，因此在傳播距離不變的情況下，傳播時間減少，故遲滯時間減少。

4　結論

（1）基于絕對節點坐標法建立了鋼索單元動力學模型，引入高階位移插值函數描述鋼索截面內變形，同時單元剛度矩陣隨結構運動變形而變化，從而實現對鋼索變形的精確描述。

（2）采用準靜態實驗驗證了ANCF鋼索單元求解的精確性，分析加載條件、結構參數和材料特性對鋼索變形和動態遲滯的影響，通過解析解驗證了模型動態遲滯預測的合理性。

（3）結果表明：在本文所討論的加載條件下，鋼索結構剛度主要受其直徑影響，初始張力對鋼索變形的影響有限，而鋼索長度與其變形程度無關；鋼索的動態響應遲滯時間由其長度及材料特性決定，與直徑無關。在保證其他條件相同的情況下，采用大直徑鋼索并適當增大初始張力，可在一定程度上增大鋼索結構剛度，減小傳動過程中的末端變形，改善長距離工況下鋼索的傳動性能。