劉忠樂, 華志宏, 薛文良
(東華大學(xué) a.理學(xué)院;b.紡織面料技術(shù)教育部重點實驗室, 上海 201620)
纖維彈性細桿模型動力學(xué)分析的有限元方法及應(yīng)用
劉忠樂a, 華志宏a, 薛文良b
(東華大學(xué) a.理學(xué)院;b.紡織面料技術(shù)教育部重點實驗室, 上海 201620)
基于紡織工程中的氣流輔助加捻紡紗, 以彈性細桿作為纖維模型, 討論纖維空間大變形的動力學(xué)問題.采用空間桿單元有限元分析方法, 對彈性細桿動力學(xué)分析的有限元模型進行簡化, 導(dǎo)出了彈性桿單元剛度方程和剛體質(zhì)量單元動力學(xué)方程, 然后綜合出彈性細桿整體的動力學(xué)運動微分方程.節(jié)點坐標姿態(tài)采用歐拉四元素表示, 導(dǎo)出了相關(guān)的坐標變換矩陣.以此方法對氣流輔助加捻紡紗動力學(xué)模型進行實例計算.
纖維模型;彈性細桿;有限元;動力學(xué)分析
隨著氣流技術(shù)在紡紗加工過程中的廣泛應(yīng)用, 在其工藝機理研究中迫切需要通過數(shù)值模擬來了解纖維在各種復(fù)雜受力狀態(tài)下的運動行為.纖維是一種具有一定彈性的且超大長徑比的柔性連續(xù)體材料, 近年來一般都將纖維簡化成珠鏈模型[1-2]進行分析計算.由于珠鏈模型不考慮纖維各截面的角位移, 故不能準確表達纖維的彎曲與扭轉(zhuǎn)特性.本文同時考慮纖維各截面的線位移和角位移, 把纖維看作連續(xù)的彈性細桿進行研究.雖然目前關(guān)于彈性細桿幾何大變形動力學(xué)問題的研究已經(jīng)有了較完整的理論[3], 但現(xiàn)有的理論和方法都是以多參數(shù)表達的聯(lián)立偏微分方程組來表述彈性細桿的力學(xué)行為, 較難用于一般空間三維問題的實際數(shù)值計算, 其實際應(yīng)用的相關(guān)研究極少.本文采用有限元方法計算彈性細桿三維空間幾何大變形動力學(xué)問題, 導(dǎo)出了彈性桿單元剛度方程和剛體質(zhì)量單元動力學(xué)方程, 然后綜合出彈性細桿整體的動力學(xué)運動微分方程, 并以此方法進行數(shù)值模擬并計算了某氣流輔助加捻紡紗工藝中纖維在氣流驅(qū)動下隨時間的運動演變過程.
采用有限元方法對彈性細桿進行動力學(xué)分析.將整個彈性細桿分割成n+1個相互獨立的微段, 在每個微段中點的橫截面中心設(shè)置一節(jié)點, 每個節(jié)點橫截面固結(jié)一連體動坐標系, 即節(jié)點坐標系.將固定于地面的慣性參考系作為有限元分析的總體坐標系.每個節(jié)點坐標系相對于固定總體坐標系有6個自由度, 設(shè)置6個坐標參數(shù)確定其相對位置, 其中3個描述節(jié)點的移動線位移, 其余3個描述節(jié)點橫截面轉(zhuǎn)動角位移.考慮彈性細桿的質(zhì)量時, 將每個微段的質(zhì)量看作隨其節(jié)點坐標系作空間一般運動的剛體(見圖1).考慮彈性細桿的彈性時, 將每相鄰兩節(jié)點之間看作為只計彈性不計質(zhì)量的空間彈性桿單元的相連接(見圖2).作用于彈性細桿上的外力都向其節(jié)點簡化.
圖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
彈性細桿動力學(xué)分析有限元模型中隨某節(jié)點i作空間一般運動的微段剛體如圖1所示.圖中Oxyz為固定總體坐標系,O′uvw為固結(jié)于剛體的連體動坐標系, 即節(jié)點坐標系, 其坐標原點位于剛體質(zhì)心, 坐標軸方向為該剛體質(zhì)心慣性主軸方向,u軸為微段剛體的軸線方向,v和w軸位于微段剛體過質(zhì)心的橫截面上.
剛體在空間的任意一般運動可分解為隨其質(zhì)心的平移和相對于質(zhì)心的轉(zhuǎn)動.根據(jù)質(zhì)心運動定理知, 剛體隨質(zhì)心O′平移的動力學(xué)方程為
ma=F
(1)
式中:m為剛體質(zhì)量;a為質(zhì)心相對于固定總體坐標系Oxyz的絕對加速度;F為剛體上外力的主矢.
根據(jù)動量矩定理[4], 可導(dǎo)出剛體繞質(zhì)心O′轉(zhuǎn)動的動力學(xué)方程為
Jα+ω×(J·ω)=M
(2)
式中:J為剛體對質(zhì)心的慣量;M為剛體上外力對質(zhì)心的主矩;ω為剛體的角速度;α為ω相對動坐標系O′uvw對時間的導(dǎo)數(shù).可以證明α也等于ω相對定坐標系Oxyz對時間的導(dǎo)數(shù).
合并式(1)和(2), 可將節(jié)點i微段剛體的動力學(xué)方程寫成矩陣形式為
(3)
或簡寫為
(4)
相對于定坐標系Oxyz, 式中
(5)
(6)
(7)
其質(zhì)量矩陣為
(8)
設(shè)微段剛體相對質(zhì)心慣性主軸O′uvw的慣量為
(9)
則
(10)
式中O′uvw相對于Oxyz的方向余弦矩陣為
(11)
圖2中的任意某相鄰的兩個節(jié)點及其之間的彈性桿單元如圖3所示, 其中i和j為彈性細桿在任意時間t的節(jié)點位置.取節(jié)點i的節(jié)點坐標系作為該單元的局部坐標系, 坐標面O′vw位于節(jié)點i的橫截面上,i-j0位置為桿單元無受力變形的參考位置.
圖3 彈性桿單元和局部坐標系Fig.3 Elastic bar element and local coordinate system
彈性桿單元的桿端變形位移與桿端力相對于總體坐標系Oxyz的關(guān)系可表示為
(12)
式中單元相對Oxyz的桿端力和桿端變形位移為
而單元剛度矩陣
(13)
現(xiàn)因取節(jié)點i的節(jié)點坐標系O′uvw作為該單元的局部坐標系(見圖3), 故有
(14)
而
(15)
記
則由式(12)可將桿單元e的剛度方程寫成
(16)
由此即可根據(jù)各節(jié)點坐標位置得出各彈性桿單元的桿端力.
對于固結(jié)于每個節(jié)點i的微段剛體, 式(4)為對應(yīng)于每個節(jié)點的6個自由度表示的動力學(xué)方程.現(xiàn)將其節(jié)點自由度編號轉(zhuǎn)換為總體節(jié)點自由度編號, 將式(4)改寫成
(17)
對于每個彈性桿單元e, 將式(16)表示的每個彈性桿單元的節(jié)點力按對應(yīng)于總體節(jié)點自由度編號進行轉(zhuǎn)換,將其改寫成
(18)
(19)
將式(17)和(18)代入式(19), 得
(20)
式中:
(21)
(22)
(23)
式(20)即為彈性細桿有限元總體動力學(xué)分析的運動微分方程.
根據(jù)剛體有限轉(zhuǎn)動的歐拉定理, 圖1中任意節(jié)點坐標系O′uvw相對于總體坐標系Oxyz的姿態(tài)方位, 可看作是過O′點的xyz方位繞過O′的某單位矢量p轉(zhuǎn)動某θ角而實現(xiàn).設(shè)單位矢量p在Oxyz坐標系的坐標為(p1,p2,p3), 可定義以下4個變量[7]
(24)
作為歐拉四元素坐標, 用于表示節(jié)點坐標系的姿態(tài)方位.
直接驗算可證實歐拉四元素之間存在如下關(guān)系
(25)
因此, 歐拉四元素中只有3個獨立變量.
可以導(dǎo)出以歐拉四元素表示的節(jié)點坐標系O′uvw相對于總體坐標系Oxyz的角速度為
(26)
圖3和式(15)中節(jié)點j的方位相對i的微小角位移為
(27)
式(11)所示的方向余弦矩陣為
(28)
由此對有限元總體動力學(xué)分析的運動微分方程式(20)進行數(shù)值積分, 就可求得整個彈性細桿各節(jié)點在不同時刻的速度和位置.
紡織工程中的一種氣流輔助加捻過程簡化模型如圖4所示.圖4中大圓柱為固定的空心腔體, 內(nèi)部中心圓柱體CD代表中心須條, 被看作相對固定的剛體, 偏心小圓柱AB為一根做包纏運動的纖維, 被看作彈性細桿, 其A端固定,B端自由, 假設(shè)在圖示初始位置的速度為零.在腔體內(nèi)存在有做高速螺旋運動的旋轉(zhuǎn)氣流, 其氣體繞中心x軸高速旋轉(zhuǎn), 同時沿著x軸方向低速移動.現(xiàn)分析包纏纖維在高速氣流的驅(qū)動下逐漸纏繞到中心須條上的動力學(xué)過程.
圖4 氣流輔助加捻紡紗的動力學(xué)模型Fig.4 Dynamic model of airflow assisted twisting yarn
采用本文的有限元方法進行動力學(xué)分析, 將整個包纏纖維均分成n+1個微段, 對于氣流作用于包纏纖維的分布荷載向節(jié)點簡化.包纏纖維與剛性中心圓柱及腔體之間的接觸約束采用罰函數(shù)法處理, 即在兩者之間恰當(dāng)?shù)貥?gòu)造一個排斥力, 該力在兩者之間的距離較接近時快速增大, 而存在一定距離后就快速衰減到零.運用Matlab編程計算, 代入具體數(shù)據(jù), 求得包纏纖維在一系列不同時刻的包纏形態(tài)位置如圖5所示.
圖5 包纏纖維隨時間變化的包纏過程Fig.5 The wrapping process of wrapping fiber with time
為研究纖維空間運動動力學(xué)過程, 本文以彈性細桿作為纖維模型, 采用有限元方法建立了彈性細桿總體動力學(xué)分析的運動微分方程, 并以此方法模擬了一種氣流輔助加捻過程中纖維在氣流驅(qū)動下隨時間的動力學(xué)運動演變過程, 可為其紡紗設(shè)備與工藝的設(shè)計與改進提供理論基礎(chǔ).本文為紡織工程中數(shù)值模擬纖維在各種不同氣流驅(qū)動下的動力學(xué)運動過程的建模提供了理論方法.
[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] 張勇, 曾泳春, 王云俠, 等.基于珠-桿模型的噴氣渦流紡噴嘴氣流場中的纖維運動規(guī)律[J].東華大學(xué)學(xué)報(自然科學(xué)版), 2013, 39(5): 583-589.
[3] 劉延柱.彈性細桿的非線性力學(xué)[M].北京: 清華大學(xué)出版社, 2006: 124-150.
[4] 劉延柱.高等動力學(xué)[M].北京: 高等教育出版社, 2001: 105- 118.
[5] 趙經(jīng)文, 王宏鈺.結(jié)構(gòu)有限元分析[M].北京: 科學(xué)出版社, 2001: 17-20.
[6] 李二明, 華志宏, 薛文良, 等.纖維彈性細桿模型幾何大變形靜力學(xué)分析的有限元方法[J].東華大學(xué)學(xué)報(自然科學(xué)版), 2016, 42(6):916-921.
[7] 洪嘉振.計算多體系統(tǒng)動力學(xué)[M].北京: 高等教育出版社, 1999: 44-46.
(責(zé)任編輯: 楊 靜)
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—),男,湖北石首人,碩士研究生,研究方向為固體力學(xué).E-mail: 2141386@mail.dhu.edu.cn 華志宏(聯(lián)系人),男,副教授,E-mail: hzh@dhu.edu.cn
TS 101.2
A