錢彥懿　王　慧　余海東

1.上海交通大學(xué)上海市復(fù)雜薄板結(jié)構(gòu)數(shù)字化制造重點實驗室，上海，200240 2.航空工業(yè)第一飛機設(shè)計研究院，西安，710089

0　引言

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

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

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

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

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

1.1　鋼索單元變形描述

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

式中，x為變形前點P在單元坐標系xyz中的位置坐標矢量；S為單元形函數(shù)矩陣；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為插值多項式的待定系數(shù)。

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

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

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

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

1.2　鋼索單元質(zhì)量矩陣

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

鋼索結(jié)構(gòu)的動能

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

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

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

1.3　鋼索單元廣義外力矩陣

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

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

式中,g為重力加速度。

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

鋼索單元的廣義彈性力可由連續(xù)介質(zhì)力學(xué)中非線性應(yīng)變位移關(guān)系推導(dǎo)得到，鋼索單元的應(yīng)變張量可用鋼索單元變形梯度J描述。對于鋼索單元模型，考慮鋼索在預(yù)緊力下運動及變形特點，其變形梯度

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

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

考慮鋼索結(jié)構(gòu)的大變形問題，利用單元變形梯度，由彈性力引起的鋼索單元應(yīng)變可由格林拉格朗日應(yīng)變張量定義：

εa為對稱張量，可將其寫為由6個獨立應(yīng)變分量構(gòu)成的列陣ε：

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

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

根據(jù)材料本構(gòu)模型，得到單元應(yīng)力

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

2　鋼索變形仿真計算分析

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

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

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

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

2.1　仿真與實驗結(jié)果比較

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

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

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

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

表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　鋼索幾何參數(shù)對其變形的影響

2.3.1不同直徑鋼索的變形

為得到鋼索直徑對其變形的影響規(guī)律，保持其余參數(shù)不變，分別計算直徑為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可以看出，當(dāng)加載條件相同時，鋼索的變形量隨直徑的增大，呈非線性減小趨勢。這是由于鋼索直徑的增大導(dǎo)致其結(jié)構(gòu)剛度增大，因此在相同加載條件下，鋼索變形量減小。

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

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

2.3.2不同長度鋼索的變形

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

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

3　鋼索動力學(xué)仿真計算分析

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

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

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

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

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

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

3.1　鋼索結(jié)構(gòu)參數(shù)對其遲滯的影響

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

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

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

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

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

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

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

鋼索材料屬性影響鋼索結(jié)構(gòu)剛度，因而對其動態(tài)特性存在影響。作為初步研究，鋼索材料密度分別取 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可以看出，仿真得到的不同材料密度及彈性模量的鋼索動態(tài)響應(yīng)遲滯時間曲線均與其對應(yīng)的理論遲滯時間曲線相吻合，最大誤差百分比分別為4.6%和10%。鋼索密度及彈性模量均對其動態(tài)響應(yīng)有較大影響，遲滯時間隨著鋼索密度的增大，呈近似線性增大趨勢，密度越大遲滯時間增加越緩慢；彈性模量對遲滯時間的影響趨勢剛好相反，隨著彈性模量的增大，遲滯時間呈非線性減少趨勢，并且減少速率逐漸趨緩。依據(jù)應(yīng)力波傳播理論和式（19）可知，應(yīng)力波的傳播速度由材料密度及彈性模量決定，彈性模量增大或材料密度減小導(dǎo)致傳播速度增大，因此在傳播距離不變的情況下，傳播時間減少，故遲滯時間減少。

4　結(jié)論

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

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

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