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

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

(東華大學 a． 理學院； b． 紡織面料技術教育部重點實驗室， 上海 201620)

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

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

(東華大學 a． 理學院； b． 紡織面料技術教育部重點實驗室， 上海 201620)

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

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

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

1 空間桿單元剛度矩陣

桿單元在三維空間的變形可分解為軸向拉壓變形、兩個主平面內的彎曲變形和扭轉變形的組合．當彈性細桿發生幾何大變形或軸力較大時，必須考慮軸力對彎曲變形的影響．

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

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

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

桿單元e節點i和節點j的桿端位移可以表示為

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

(1)

式中：Δui, Δvi, Δwi和Δuj, Δvj, Δwj為節點線位移；Δθui, Δθvi, Δθwi和Δθuj, Δθvj, Δθwj為橫截面角位移．

與桿端位移對應的桿端力可表示為

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為節點力；Mu i, Mv i, Mw i和Mu j, Mv j, Mw j為節點力偶．

1.2 桿單元的彎曲變形分析

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

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

(3)

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

(4)

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

(5)

圖2 主平面內的彎曲變形Fig.2 Bending deformation in the principal plane

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

(6)

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

(7)

通過對u積分，并根據桿單元節點邊界條件及其桿單元力平衡方程可解得

(8)

同理，在uO′w平面內的彎曲變形也可導出類似的桿端力與桿端位移關系．

1.3 桿單元的剛度矩陣

對于桿的軸向拉壓變形和扭轉變形，可導出下列關系

(9)

(10)

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

(11)

而空間桿單元的桿端力和桿端位移關系可簡寫為

(12)

2 坐標變換

2.1 剛度矩陣的坐標變換

在圖1所示總體坐標系中，單個節點的位移為沿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)

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

(15)

其中

(16)

為桿單元e的單元坐標變換矩陣，而

(17)

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

(18)

其中

(19)

即為桿單元相對于總體坐標系的剛度矩陣．

2.2 坐標系方位的歐拉角坐標表示

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

(20)

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

(21)

(22)

3 非線性問題的分步加載計算法

3.1 分步加載逐次逼近概念

纖維彈性細桿模型超大位移變形問題為幾何非線性問題．由式(12)或式(18)表示的單元剛度方程只有在桿單元相對于單元局部坐標系的位移足夠小的情況下才成立．為此，必須保證桿單元的長度足夠的小，桿單元的局部坐標系的位置距離彈性細桿大位移變形后的精確位置足夠的近．

現采用載荷分步加載逐次逼近的方法求解．載荷從零開始由小到大分若干次逐漸加大進行分步計算，并保證每一步載荷的增量為一小量．在每步計算中，都以上一步的載荷作用下求得的位置精確值，作為求解下一步的載荷作用下的位置精確值的初始近似值，來求解下一步位置的精確值．對于載荷為零開始的第一加載步計算，將整個彈性細桿無載荷狀態的位置作為初始位置，在該位置為每個桿單元建立單元局部坐標，然后載荷增加一個小量，用相對于單元局部坐標小位移情況下成立的單元剛度方程求解在該載荷作用下的各桿單元的位移，從而得到整個彈性細桿變形后的位置；對于任意第N+1加載步計算，將第N步求得的彈性細桿位置作為初始位置，在該位置再為每個桿單元建立單元局部坐標，將載荷在第N步的基礎上再增大一個小量，繼續用相對于單元局部坐標小位移情況下成立的單元剛度方程式求解在該新的載荷作用下的各桿單元的位移，從而得到整個彈性細桿的一個新的位置．

3.2 逐次逼近計算式

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

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

由圖4可看出

(23)

其中，相對于總體坐標系O-xyz

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

(24)

將式(24)代入單元剛度方程式(18)，得每一加載步的單元位移增量方程如下

(25)

其中

4 有限元總體分析

4.1 整體位移增量方程

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

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

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

(26)

將圖5所示的全部單元組合在一起，其全部節點的靜力學平衡方程為

(27)

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

或簡寫成如下形式

(28)

4.2 算例

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

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

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

運用Matlab編程計算，代入具體數據，求得纖維包纏的靜平衡形態如圖7所示．

(a) 空間三維視圖

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

5 結 語

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

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

[1] 郭會芬．噴氣紡紗噴嘴內三維旋轉氣流場及柔性纖維運動的研究[D]．上海：東華大學紡織學院，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] 劉延柱．彈性細桿的非線性力學[M]．北京：清華大學出版社，2006：1452．

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

[5] 洪嘉振．計算多體系統動力學[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

上海市自然科學基金資助項目(13ZR1400900)

李二明(1986—)，男，河南駐馬店人，碩士研究生，研究方向為固體力學．E-mail: 2131410@mail.dhu.edu.cn 華志宏(聯系人)，男，副教授，E-mail: hzh@dhu.edu.cn

TS 101.2

A