纖維彈性細桿模型動力學分析的有限元方法及應用

劉忠樂, 華志宏, 薛文良

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

纖維彈性細桿模型動力學分析的有限元方法及應用

劉忠樂a, 華志宏a, 薛文良b

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

基于紡織工程中的氣流輔助加捻紡紗, 以彈性細桿作為纖維模型, 討論纖維空間大變形的動力學問題．采用空間桿單元有限元分析方法, 對彈性細桿動力學分析的有限元模型進行簡化, 導出了彈性桿單元剛度方程和剛體質量單元動力學方程, 然后綜合出彈性細桿整體的動力學運動微分方程．節點坐標姿態采用歐拉四元素表示, 導出了相關的坐標變換矩陣．以此方法對氣流輔助加捻紡紗動力學模型進行實例計算．

纖維模型；彈性細桿；有限元；動力學分析

隨著氣流技術在紡紗加工過程中的廣泛應用, 在其工藝機理研究中迫切需要通過數值模擬來了解纖維在各種復雜受力狀態下的運動行為．纖維是一種具有一定彈性的且超大長徑比的柔性連續體材料, 近年來一般都將纖維簡化成珠鏈模型[1-2]進行分析計算．由于珠鏈模型不考慮纖維各截面的角位移, 故不能準確表達纖維的彎曲與扭轉特性．本文同時考慮纖維各截面的線位移和角位移, 把纖維看作連續的彈性細桿進行研究．雖然目前關于彈性細桿幾何大變形動力學問題的研究已經有了較完整的理論[3], 但現有的理論和方法都是以多參數表達的聯立偏微分方程組來表述彈性細桿的力學行為, 較難用于一般空間三維問題的實際數值計算, 其實際應用的相關研究極少．本文采用有限元方法計算彈性細桿三維空間幾何大變形動力學問題, 導出了彈性桿單元剛度方程和剛體質量單元動力學方程, 然后綜合出彈性細桿整體的動力學運動微分方程, 并以此方法進行數值模擬并計算了某氣流輔助加捻紡紗工藝中纖維在氣流驅動下隨時間的運動演變過程．

1 彈性細桿動力學分析的有限元模型 簡化

采用有限元方法對彈性細桿進行動力學分析．將整個彈性細桿分割成n+1個相互獨立的微段, 在每個微段中點的橫截面中心設置一節點, 每個節點橫截面固結一連體動坐標系, 即節點坐標系．將固定于地面的慣性參考系作為有限元分析的總體坐標系．每個節點坐標系相對于固定總體坐標系有6個自由度, 設置6個坐標參數確定其相對位置, 其中3個描述節點的移動線位移, 其余3個描述節點橫截面轉動角位移．考慮彈性細桿的質量時, 將每個微段的質量看作隨其節點坐標系作空間一般運動的剛體(見圖1)．考慮彈性細桿的彈性時, 將每相鄰兩節點之間看作為只計彈性不計質量的空間彈性桿單元的相連接(見圖2)．作用于彈性細桿上的外力都向其節點簡化．

圖1 作空間一般運動的微段剛體Fig.1 Rigid body unit for general motion in space

圖2 彈性細桿的彈性桿單元有限元模型Fig.2 Finite element model of elastic bar elements of elastic thin rod

2 微段剛體空間一般運動的動力學方程

彈性細桿動力學分析有限元模型中隨某節點i作空間一般運動的微段剛體如圖1所示．圖中Oxyz為固定總體坐標系,O′uvw為固結于剛體的連體動坐標系, 即節點坐標系, 其坐標原點位于剛體質心, 坐標軸方向為該剛體質心慣性主軸方向,u軸為微段剛體的軸線方向,v和w軸位于微段剛體過質心的橫截面上．

剛體在空間的任意一般運動可分解為隨其質心的平移和相對于質心的轉動．根據質心運動定理知, 剛體隨質心O′平移的動力學方程為

ma=F

(1)

式中:m為剛體質量；a為質心相對于固定總體坐標系Oxyz的絕對加速度；F為剛體上外力的主矢．

根據動量矩定理[4], 可導出剛體繞質心O′轉動的動力學方程為

Jα+ω×(J·ω)=M

(2)

式中:J為剛體對質心的慣量；M為剛體上外力對質心的主矩;ω為剛體的角速度；α為ω相對動坐標系O′uvw對時間的導數.可以證明α也等于ω相對定坐標系Oxyz對時間的導數.

合并式(1)和(2), 可將節點i微段剛體的動力學方程寫成矩陣形式為

(3)

或簡寫為

(4)

相對于定坐標系Oxyz, 式中

(5)

(6)

(7)

其質量矩陣為

(8)

設微段剛體相對質心慣性主軸O′uvw的慣量為

(9)

則

(10)

式中O′uvw相對于Oxyz的方向余弦矩陣為

(11)

3 以節點坐標位置表示彈性桿單元的 桿端力

圖2中的任意某相鄰的兩個節點及其之間的彈性桿單元如圖3所示, 其中i和j為彈性細桿在任意時間t的節點位置．取節點i的節點坐標系作為該單元的局部坐標系, 坐標面O′vw位于節點i的橫截面上,i-j0位置為桿單元無受力變形的參考位置．

圖3 彈性桿單元和局部坐標系Fig.3 Elastic bar element and local coordinate system

彈性桿單元的桿端變形位移與桿端力相對于總體坐標系Oxyz的關系可表示為

(12)

式中單元相對Oxyz的桿端力和桿端變形位移為

而單元剛度矩陣

(13)

現因取節點i的節點坐標系O′uvw作為該單元的局部坐標系(見圖3), 故有

(14)

而

(15)

記

則由式(12)可將桿單元e的剛度方程寫成

(16)

由此即可根據各節點坐標位置得出各彈性桿單元的桿端力．

4 彈性細桿有限元總體動力學分析

對于固結于每個節點i的微段剛體, 式(4)為對應于每個節點的6個自由度表示的動力學方程．現將其節點自由度編號轉換為總體節點自由度編號, 將式(4)改寫成

(17)

對于每個彈性桿單元e, 將式(16)表示的每個彈性桿單元的節點力按對應于總體節點自由度編號進行轉換,將其改寫成

(18)

(19)

將式(17)和(18)代入式(19), 得

(20)

式中：

(21)

(22)

(23)

式(20)即為彈性細桿有限元總體動力學分析的運動微分方程．

5 節點坐標系姿態的歐拉四元素坐標 表示

根據剛體有限轉動的歐拉定理, 圖1中任意節點坐標系O′uvw相對于總體坐標系Oxyz的姿態方位, 可看作是過O′點的xyz方位繞過O′的某單位矢量p轉動某θ角而實現．設單位矢量p在Oxyz坐標系的坐標為(p1,p2,p3), 可定義以下4個變量[7]

(24)

作為歐拉四元素坐標, 用于表示節點坐標系的姿態方位.

直接驗算可證實歐拉四元素之間存在如下關系

(25)

因此, 歐拉四元素中只有3個獨立變量.

可以導出以歐拉四元素表示的節點坐標系O′uvw相對于總體坐標系Oxyz的角速度為

(26)

圖3和式(15)中節點j的方位相對i的微小角位移為

(27)

式(11)所示的方向余弦矩陣為

(28)

由此對有限元總體動力學分析的運動微分方程式(20)進行數值積分, 就可求得整個彈性細桿各節點在不同時刻的速度和位置．

6 實 例

紡織工程中的一種氣流輔助加捻過程簡化模型如圖4所示．圖4中大圓柱為固定的空心腔體, 內部中心圓柱體CD代表中心須條, 被看作相對固定的剛體, 偏心小圓柱AB為一根做包纏運動的纖維, 被看作彈性細桿, 其A端固定,B端自由, 假設在圖示初始位置的速度為零．在腔體內存在有做高速螺旋運動的旋轉氣流, 其氣體繞中心x軸高速旋轉, 同時沿著x軸方向低速移動．現分析包纏纖維在高速氣流的驅動下逐漸纏繞到中心須條上的動力學過程．

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

采用本文的有限元方法進行動力學分析, 將整個包纏纖維均分成n+1個微段, 對于氣流作用于包纏纖維的分布荷載向節點簡化．包纏纖維與剛性中心圓柱及腔體之間的接觸約束采用罰函數法處理, 即在兩者之間恰當地構造一個排斥力, 該力在兩者之間的距離較接近時快速增大, 而存在一定距離后就快速衰減到零．運用Matlab編程計算, 代入具體數據, 求得包纏纖維在一系列不同時刻的包纏形態位置如圖5所示．

圖5 包纏纖維隨時間變化的包纏過程Fig.5 The wrapping process of wrapping fiber with time

7 結 語

為研究纖維空間運動動力學過程, 本文以彈性細桿作為纖維模型, 采用有限元方法建立了彈性細桿總體動力學分析的運動微分方程, 并以此方法模擬了一種氣流輔助加捻過程中纖維在氣流驅動下隨時間的動力學運動演變過程, 可為其紡紗設備與工藝的設計與改進提供理論基礎．本文為紡織工程中數值模擬纖維在各種不同氣流驅動下的動力學運動過程的建模提供了理論方法．

[1] YAMAMOTO S, MATSUOKA T. A method for dynamic simulation of rigid and flexible fiber in a flow field [J]. Journal of Chemical Physics, 1993, 98(1): 644-650.

[2] 張勇, 曾泳春, 王云俠, 等．基于珠-桿模型的噴氣渦流紡噴嘴氣流場中的纖維運動規律[J]．東華大學學報(自然科學版), 2013, 39(5): 583-589．

[3] 劉延柱．彈性細桿的非線性力學[M]．北京: 清華大學出版社, 2006: 124-150．

[4] 劉延柱．高等動力學[M]．北京: 高等教育出版社, 2001: 105- 118．

[5] 趙經文, 王宏鈺．結構有限元分析[M]．北京: 科學出版社, 2001: 17-20．

[6] 李二明, 華志宏, 薛文良, 等．纖維彈性細桿模型幾何大變形靜力學分析的有限元方法[J]．東華大學學報(自然科學版), 2016, 42(6)：916-921．

[7] 洪嘉振．計算多體系統動力學[M]．北京: 高等教育出版社, 1999: 44-46．

(責任編輯: 楊 靜)

Finite-Element Method and Application of Dynamic Analysis of Fiber Elastic Thin Rod Model

LIUZhonglea,HUAZhihonga,XUEWenliangb

(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 textile engineering, the dynamic problems of this model’s spatial large deformation were discussed using the elastic thin rod as the fiber model. Using spatial bar finite-element method for the dynamic analysis of elastic thin rod model, the spatial bar element stiffness equation and the rigid body mass unit dynamic equation were derived, and then the kinetic differential equations of whole elastic thin rod were synthesized. Also, the related coordinate transformation matrices were obtained using the Euler quaternion for the determination of node coordinate attitude. As an example, the method was used in calculation of the airflow assisted twisting spinning dynamics model.

fiber model; elastic thin rod; finite-element; dynamic analysis

1671-0444 (2017)03-0454-05

2016-06-13

劉忠樂(1993—)，男，湖北石首人，碩士研究生，研究方向為固體力學．E-mail: 2141386@mail.dhu.edu.cn 華志宏(聯系人)，男，副教授，E-mail: hzh@dhu.edu.cn

TS 101.2

A