李二明，華志宏，薛文良，程隆棣

(東華大學(xué) a． 理學(xué)院； b． 紡織面料技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室， 上海 201620)

纖維彈性細(xì)桿模型幾何大變形靜力學(xué)分析的有限元方法

李二明a，華志宏a，薛文良b，程隆棣b

(東華大學(xué) a． 理學(xué)院； b． 紡織面料技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室， 上海 201620)

以紡織工程中的氣流輔助加捻紡紗為背景，以彈性細(xì)桿為纖維模型，對其幾何位移大變形的非線性問題進(jìn)行分析和討論．采用空間桿單元有限元分析方法，導(dǎo)出了考慮軸力對彎曲影響的空間桿單元?jiǎng)偠染仃?，空間方位用歐拉角表示，幾何非線性采用分步加載逐次逼近方法處理，建立了彈性細(xì)桿靜力平衡方程組．對纖維彈性細(xì)桿有限元總體分析，運(yùn)用Matlab模擬氣流輔助加捻紡紗過程中纖維包纏的平衡狀態(tài)．

纖維；彈性細(xì)桿；有限元；大變形；分步加載

紡織工程各種紡紗工藝的成紗機(jī)理研究中，迫切需要了解纖維在各種復(fù)雜受力狀態(tài)下的運(yùn)動(dòng)行為表現(xiàn)．纖維是一種具有一定彈性的具有超大長徑比的柔性連續(xù)體材料．在以往的研究中為便于計(jì)算，通常把纖維簡化為較簡單的各種力學(xué)模型，例如把纖維簡化為是由許多剛性短桿相互鉸鏈串接的桿鏈模型[1 - 2]，這些模型與纖維的實(shí)際力學(xué)行為之間都存在著較大的失真．本文不作簡化直接把纖維看作連續(xù)的彈性細(xì)桿進(jìn)行計(jì)算研究．文獻(xiàn)[3]系統(tǒng)地介紹了研究彈性細(xì)桿的非線性力學(xué)方法，目前都以多參數(shù)表達(dá)的聯(lián)立偏微分方程組來表述彈性細(xì)桿的力學(xué)行為，較難用于一般空間三維問題的實(shí)際數(shù)值計(jì)算．本文采用空間桿單元有限元方法計(jì)算彈性細(xì)桿的空間靜力學(xué)問題．通常的桿單元有限元法只能處理系統(tǒng)小變形的線性問題，而紡織工程中的纖維受力變形屬于幾何大變形非線性問題．本文推導(dǎo)出了同時(shí)考慮桿的拉壓、彎曲和扭轉(zhuǎn)并計(jì)入拉壓軸力對彎曲影響的空間桿單元?jiǎng)偠染仃?，對于幾何大變形非線性問題采用分步加載逐次逼近的方法進(jìn)行計(jì)算，并以此方法計(jì)算了氣流輔助加捻紡紗過程中纖維包纏狀態(tài)分析．

1 空間桿單元?jiǎng)偠染仃?/h2>

桿單元在三維空間的變形可分解為軸向拉壓變形、兩個(gè)主平面內(nèi)的彎曲變形和扭轉(zhuǎn)變形的組合．當(dāng)彈性細(xì)桿發(fā)生幾何大變形或軸力較大時(shí)，必須考慮軸力對彎曲變形的影響．

1.1 單元的桿端位移和桿端力

將整個(gè)空間彈性細(xì)桿劃分成若干微小桿單元．為整個(gè)彈性細(xì)桿建立一個(gè)總體坐標(biāo)系O -xyz，再為每個(gè)桿單元建立一個(gè)單元局部坐標(biāo)系O′-uvw，如圖1所示．其u軸為桿單元的軸線方向，uO′v和uO′w為桿單元的兩個(gè)彎曲主平面方位．

圖1 空間桿單元及其坐標(biāo)系Fig.1 The spatial bar element and coordinate system

桿單元e節(jié)點(diǎn)i和節(jié)點(diǎn)j的桿端位移可以表示為

Δuj,Δvj,Δwj,Δθuj,Δθvj,Δθwj}T

(1)

式中：Δui, Δvi, Δwi和Δuj, Δvj, Δwj為節(jié)點(diǎn)線位移；Δθui, Δθvi, Δθwi和Δθuj, Δθvj, Δθwj為橫截面角位移．

與桿端位移對應(yīng)的桿端力可表示為

Fu j, Fv j, Fw j, Mu j, Mv j, Mw j}T

(2)

式中：Fu i, Fv i, Fw i和Fu j, Fv j, Fw j為節(jié)點(diǎn)力；Mu i, Mv i, Mw i和Mu j, Mv j, Mw j為節(jié)點(diǎn)力偶．

1.2 桿單元的彎曲變形分析

分析圖2所示uO′v平面內(nèi)的彎曲變形．考慮桿的軸力對彎矩的影響，桿單元任意截面位置u的彎矩為

M=Fv iu+Fn(v-Δvi)-Mw i

(3)

式中：v為桿的撓度；Fn為桿的軸力，F(xiàn)n=-Fui=Fuj．假設(shè)桿單元的剛性位移遠(yuǎn)大于其彈性變形，式(3)中的v近似取

(4)

式中：l為桿單元長度．式(4)代入式(3)得

(5)

圖2 主平面內(nèi)的彎曲變形Fig.2 Bending deformation in the principal plane

假設(shè)相對于局部坐標(biāo)系O′-uv，桿單元的位移為小量，其撓曲線的一階導(dǎo)數(shù)很小，則桿的彎曲撓曲線近似微分方程[4]為

(6)

式中：Iw為截面慣性矩；E為材料彈性模量．式(5)代入式(6)得

(7)

通過對u積分，并根據(jù)桿單元節(jié)點(diǎn)邊界條件及其桿單元力平衡方程可解得

(8)

同理，在uO′w平面內(nèi)的彎曲變形也可導(dǎo)出類似的桿端力與桿端位移關(guān)系．

1.3 桿單元的剛度矩陣

對于桿的軸向拉壓變形和扭轉(zhuǎn)變形，可導(dǎo)出下列關(guān)系

(9)

(10)

式中：A為橫截面面積；G為材料剪切彈性模量；Ip為截面極慣性矩．

(11)

而空間桿單元的桿端力和桿端位移關(guān)系可簡寫為

(12)

2 坐標(biāo)變換

2.1 剛度矩陣的坐標(biāo)變換

在圖1所示總體坐標(biāo)系中，單個(gè)節(jié)點(diǎn)的位移為沿x、y、z軸的線位移和繞x、y、z軸的角位移，桿單元的桿端位移和桿端力為

{Δ}e={ Δxi, Δyi, Δzi, Δθxi, Δθyi, Δθzi,

Δxj, Δyj, Δzj, Δθxj, Δθyj, Δθzj}T

(13)

Fxj, Fyj, Fzj, Mxj, Myj, Mzj}T

(14)

在總體坐標(biāo)系O-xyz和局部坐標(biāo)系O′-uvw之間存在如下坐標(biāo)變換

(15)

其中

(16)

為桿單元e的單元坐標(biāo)變換矩陣，而

(17)

為單元局部坐標(biāo)系O′-uvw相對于總體坐標(biāo)系O-xyz的方向余弦矩陣．

(18)

其中

(19)

即為桿單元相對于總體坐標(biāo)系的剛度矩陣．

2.2 坐標(biāo)系方位的歐拉角坐標(biāo)表示

如圖3所示，將單元局部坐標(biāo)系O′-uvw的方位看作是坐標(biāo)系O-xyz依次繞圖中z軸、x1軸和z2軸的3個(gè)轉(zhuǎn)動(dòng)角ψ1、 ψ2和ψ3實(shí)現(xiàn)的．(ψ1, ψ2, ψ3)通常被稱為歐拉坐標(biāo)．以歐拉角表示的坐標(biāo)方向余弦矩陣[5]如下

(20)

圖3 歐拉角Fig.3 Euler’s angles

(21)

(22)

3 非線性問題的分步加載計(jì)算法

3.1 分步加載逐次逼近概念

纖維彈性細(xì)桿模型超大位移變形問題為幾何非線性問題．由式(12)或式(18)表示的單元?jiǎng)偠确匠讨挥性跅U單元相對于單元局部坐標(biāo)系的位移足夠小的情況下才成立．為此，必須保證桿單元的長度足夠的小，桿單元的局部坐標(biāo)系的位置距離彈性細(xì)桿大位移變形后的精確位置足夠的近．

現(xiàn)采用載荷分步加載逐次逼近的方法求解．載荷從零開始由小到大分若干次逐漸加大進(jìn)行分步計(jì)算，并保證每一步載荷的增量為一小量．在每步計(jì)算中，都以上一步的載荷作用下求得的位置精確值，作為求解下一步的載荷作用下的位置精確值的初始近似值，來求解下一步位置的精確值．對于載荷為零開始的第一加載步計(jì)算，將整個(gè)彈性細(xì)桿無載荷狀態(tài)的位置作為初始位置，在該位置為每個(gè)桿單元建立單元局部坐標(biāo)，然后載荷增加一個(gè)小量，用相對于單元局部坐標(biāo)小位移情況下成立的單元?jiǎng)偠确匠糖蠼庠谠撦d荷作用下的各桿單元的位移，從而得到整個(gè)彈性細(xì)桿變形后的位置；對于任意第N+1加載步計(jì)算，將第N步求得的彈性細(xì)桿位置作為初始位置，在該位置再為每個(gè)桿單元建立單元局部坐標(biāo)，將載荷在第N步的基礎(chǔ)上再增大一個(gè)小量，繼續(xù)用相對于單元局部坐標(biāo)小位移情況下成立的單元?jiǎng)偠确匠淌角蠼庠谠撔碌妮d荷作用下的各桿單元的位移，從而得到整個(gè)彈性細(xì)桿的一個(gè)新的位置．

3.2 逐次逼近計(jì)算式

考察圖4所示的一個(gè)加載步中的任意某個(gè)單元．圖4中i -j位置為該加載步初始位置，i′-j′位置為載荷增大一小量后單元的新位置，i -j0位置為桿單元無受力變形的位置，坐標(biāo)面vO′w位于單元i端的橫截面上， Δri和Δrj為該單元的位移增量．

圖4 單元分步加載分析Fig.4 Step loading analysis of the unit

由圖4可看出

(23)

其中，相對于總體坐標(biāo)系O-xyz

上述每個(gè)位移量中均包含線位移分量和角位移分量．而單元桿端位移

(24)

將式(24)代入單元?jiǎng)偠确匠淌?18)，得每一加載步的單元位移增量方程如下

(25)

其中

4 有限元總體分析

4.1 整體位移增量方程

如圖5所示，一空間彈性細(xì)桿被劃分為n個(gè)桿單元，共有n+1個(gè)節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)有6個(gè)自由度，全部節(jié)點(diǎn)共有6(n+1)個(gè)自由度，即共有6(n+1)個(gè)位移分量．

圖5 桿的單元分解Fig.5 The element dividing of the rod

對于由式(25)表示的每個(gè)單元e的單元位移增量方程，相對于總體節(jié)點(diǎn)自由度，可改寫成如下維數(shù)矩陣

(26)

將圖5所示的全部單元組合在一起，其全部節(jié)點(diǎn)的靜力學(xué)平衡方程為

(27)

將式(26)代入式(27)，得

或簡寫成如下形式

(28)

4.2 算例

將紡織工程中氣流輔助加捻紡紗過程簡化成圖6所示力學(xué)模型．圖6中大圓柱體CD代表中心須條，看作相對固定的剛體．彈性細(xì)桿AB代表一根包纏纖維，在由繞中心須條旋轉(zhuǎn)的氣流場產(chǎn)生的空氣阻力q作用下包纏到中心須條上．現(xiàn)把彈性細(xì)桿A端的約束看作固定約束，把空氣阻力看作彈性細(xì)桿所受的載荷，分析計(jì)算纖維在給定載荷作用下包纏在中心須條上的靜平衡形態(tài)．

圖6 氣流輔助加捻紡紗的力學(xué)模型Fig.6 Mechanical model of airflow assisted twisting yarn

采用桿單元有限元法，將整個(gè)彈性細(xì)桿均分成n個(gè)單元．對于氣流作用于桿的分布載荷，在計(jì)算時(shí)將其向各節(jié)點(diǎn)簡化．關(guān)于彈性細(xì)桿與中心剛性圓柱之間的接觸約束，采用罰函數(shù)法處理，即在兩者之間構(gòu)造一個(gè)排斥力，該力在兩者之間的距離較接近時(shí)快速增大，而存在一定距離后就快速衰減到零，從而在分步加載逼近過程中阻止彈性細(xì)桿進(jìn)入中心剛性圓柱體內(nèi)．

運(yùn)用Matlab編程計(jì)算，代入具體數(shù)據(jù)，求得纖維包纏的靜平衡形態(tài)如圖7所示．

(a) 空間三維視圖

(b) yz平面視圖圖7 纖維包纏的靜平衡形態(tài)Fig.7 The static equilibrium configuration of the fiber wrapping

5 結(jié) 語

本文采用空間桿單元有限元方法研究了纖維彈性細(xì)桿幾何大變形靜力平衡狀態(tài)，直接對纖維彈性細(xì)桿幾何大變形進(jìn)行數(shù)值求解，得到纖維在氣流作用下包纏在中心剛性圓柱上的靜平衡位置形態(tài)， 解

決了桿單元有限元分析方法難以用于空間三維幾何大變形問題實(shí)際計(jì)算的問題．改變數(shù)值模擬原始條件，可以得出不同載荷作用下纖維彈性細(xì)桿包纏須條的狀態(tài)，為氣流輔助加捻紡紗設(shè)備的設(shè)計(jì)與改進(jìn)提供了理論基礎(chǔ)，同時(shí)也為纖維彈性細(xì)桿幾何大變形動(dòng)力學(xué)分析提供了基礎(chǔ)和方向．

[1] 郭會(huì)芬．噴氣紡紗噴嘴內(nèi)三維旋轉(zhuǎn)氣流場及柔性纖維運(yùn)動(dòng)的研究[D]．上海：東華大學(xué)紡織學(xué)院，2009：110111．

[2] KONG L X, PLATFOOT R A. Computational two-phase air/fiber flow within transfer channels of rotor spinning machines [J]. Textile Res J, 1997, 67(4): 269278.

[3] 劉延柱．彈性細(xì)桿的非線性力學(xué)[M]．北京：清華大學(xué)出版社，2006：1452．

[4] 劉鴻文．材料力學(xué)[M]．5版．北京：高等教育出版社，2011：174177．

[5] 洪嘉振．計(jì)算多體系統(tǒng)動(dòng)力學(xué)[M]．北京：高等教育出版社，1999：4647．

Finite-Element Method of Statics Analysis of the Geometric Large Deformation of Fiber Elastic Thin Rod Model

LI Er-minga，HUA Zhi-honga， XUE Wen-liangb， CHENG Long-dib

(a． College of Science； b． Key Laboratory of Textile Science & Technology， Ministry of Education，Donghua University， Shanghai 201620， China)

Based on the airflow assisted twisting of yarn spinning in the textile engineering, using the elastic thin rod as the fiber model, the geometric nonlinear problems for this model’s large deformation were analyzed and discussed. Using finite-element method and taking the bending problem influenced by the axial force into consideration, the spatial bar element stiffness matrixes were derived. The elastic thin rod static equilibrium equations were established using consecutive loading approach to deal with geometric nonlinearity problem and using Euler’s angle to denote spatial orientation. The unity of the fiber elastic thin rod was analyzed by the finite-element method,and the static equilibrium configurations of the fiber in the airflow assisted twisting yarn for textile engineering were simulated by Matlab.

fiber; elastic thin rod; finite-element; large deformation; consecutive loading

16710444 (2016)060916-06

20150914

上海市自然科學(xué)基金資助項(xiàng)目(13ZR1400900)

李二明(1986—)，男，河南駐馬店人，碩士研究生，研究方向?yàn)楣腆w力學(xué)．E-mail: 2131410@mail.dhu.edu.cn 華志宏(聯(lián)系人)，男，副教授，E-mail: hzh@dhu.edu.cn

TS 101.2

A