金玉珍， 李 俊， 林培鋒， 吳震宇

(浙江理工大學 機械與自動控制學院， 浙江 杭州 310018)

?

纖維與氣流耦合的數值模擬

金玉珍， 李 俊， 林培鋒， 吳震宇

(浙江理工大學 機械與自動控制學院， 浙江 杭州 310018)

結合纖維的本構方程與氣流控制方程，建立二維的纖維與氣流耦合方程。采用相容時間積分與迭代耦合算法，對纖維在噴氣織機主噴嘴中拉伸彎曲運動特性、纖維受力情況以及與壁面接觸時的微觀運動特性進行研究。結果表明:纖維在氣流的作用下會產生拉伸彎曲變形，呈現正弦或者余弦形式波動；纖維在自由端的彎曲變形效果比在約束端的明顯；纖維受到氣流作用時, 最大應力值往往發(fā)生在纖維彎曲部位以及與壁面相接觸的部位。

流固耦合； 數值模擬； 纖維； 自適應網格

柔性體與流體耦合的數值模擬一直是流體動力型紡織機械基礎研究領域的一大難題，國內外不少學者采用各種方法對纖維與氣流兩相耦合運動特點進行研究分析。Smith等[1]研究了纖維在二維加速層流中的運動狀態(tài)。構造能夠描述纖維慣性、空氣拖曳力和彎曲剛度的非線性方程，數值模擬帶有彎鉤形狀的纖維在氣流的作用下被伸直。Kong等[2]構建了一個新的空氣/纖維二維模型來模擬纖維在約束管道中的運動。通過模擬得到了沿著氣流方向的纖維會產生彎曲，但是這個模型不能體現纖維的伸長、扭轉等特性；Yamamoto等[3]提出了一種能夠在流場中模擬剛性和柔性纖維運動的方法。纖維假設為N個相同半徑的球相互連接而成，通過改變球與球之間的距離、調度、扭轉角來獲得纖維的拉伸、彎曲和扭轉特性。通過求解流體動力和扭轉力對單個球體的位移和旋轉方程來獲得纖維模型的運動。但在具體應用上，他們忽略了纖維的扭轉，并且提出的珠鏈模型的計算成本很高。Lindstrom等[4]提出了一個柔性纖維在黏性流體中的運動模型，纖維被看成由纖維段連接而成，通過與氣流的黏性與動態(tài)拖曳力來描述。Takemura等[5]采用浸沒邊界法模擬單根纖維在低雷諾數的剪切流的運動變化。數值模擬的結果顯示了流體黏度的增加會使纖維運動變得更加復雜。Tornberg等[6]采用非局部細長體理論將流體對細絲的作用力與細絲的速度聯(lián)系起來，以此建立纖維彎曲與拉力的數學模型。采用有限差分法模擬了單根柔性細絲以及相互作用的多根細絲在剪切流中的運動，結果顯示流體剪應力會引起細絲的彎曲和松弛；Dong等[7]分別采用隨機模擬和大渦模擬這2種方法模擬了柔性纖維在三維湍流中的運動變化，考慮了纖維與壁面之間的相互作用。結果顯示纖維濃度在接近壁面附近呈線性變化，而其在遠離壁面保持恒定。裴澤光等[8]采用N-S方程加ALE方程作為氣流控制方程，并采用自適應網格的方法來改善網格畸變，研究纖維的運動變形。其研究結果能體現纖維柔性變形的特點，但對纖維的本構方程構建、纖維的變形特性未做詳細的剖析。徐存強等[9]通過數值模擬與實驗測試相結合的方法，反推出在不同相對速度下氣流與緯紗的摩擦因素，得到了氣流引緯主噴嘴內氣流速度與緯紗間摩擦因素的關系。

已有國內外不少學者對柔性體與流體耦合的數值模擬做了很多探索性的工作，但迄今為止還沒有形成公認可靠的、計算穩(wěn)定的、快速的數值模擬方法。

針對纖維-氣流耦合特性問題的研究，在原有Elastic-orthotropic材料模型基礎上，結合纖維材料的特性修正了Elastic-orthotropic材料模型的本構方程以及相關參數，基于有限體積方法求解氣流場，再對材料結構中的每個質點進行拉格朗日追蹤，流體方程和材料本構方程按順序相互迭代求解，直到兩相耦合系統(tǒng)的解達到收斂，繼而通過數據分析得到纖維材料在高速氣流中的運動特性。本文研究為流體動力型紡織機械機制特性的進一步研究提供了一定的參考。

1 數學模型

1.1 纖維本構方程

纖維具有特殊的物理特征:長徑比大，有彈性，當受到流場的作用力時會產生拉伸、彎曲等變形特性。假設纖維材料是線性彈性的，可以采用Elastic-orthotropic材料模型。二維模型中，纖維的本構方程[10]可以表示為

(1)

式中：σ11、σ22表示纖維在二維坐標軸上的2個主應力；相應的ε11、ε22為正應變；σ12代表纖維在二維坐標軸上的切向應力；相應的γ12為切向應變；kij為材料的剛度。材料的剛度可以通過材料的泊松比以及彈性模量等參數來調整。通過該本構方程就能夠得到纖維應力與應變對應的函數關系。

1.2 纖維與氣流控制方程

在二維直角坐標系下，不考慮密度變化以及熱交換可以得到流體運動控制方程[11]：

(2)

(3)

(4)

式中：v和w分別是y軸和z軸上的速度；ρ和p分別是氣流密度和壓強。根據流體控制方程給出的變量，得出了流體模型的所需的入口條件以及初始條件，本文采用壓力入口條件。

纖維的運動方程為

(5)

(6)

由于纖維被氣流場包圍，纖維與氣流相互作用的部分發(fā)生在2個區(qū)域的交界面上，要滿足以下2個條件。

運動學[12]條件為

(7)

動力學條件為

(8)

式中：df和ds分別表示纖維和氣流在流固耦合面上的位移；σf和σs分別表示纖維和氣流在流固交界面上的應力；n表示界面上的法向向量。

1.3 流固耦合相容時間積分

流體和結構模型中的時間積分必須是相容的。雖然流體模型采用的是Euler坐標系，結構模型采用的是Lagrange坐標系，但在流固耦合界面上都采用的Lagrange坐標系；因此，流體和結構模型中交界面上的位移、速度、加速度都是相同的。流體的速度和加速度[13]分別為

(9)

(10)

式中，α為時間復合積分系數。t+△t時刻的速度和加速度可以用位移未知量來表示。

(11)

(12)

式中，ξt、ηt分別為t時刻有關速度和位移的變化量。

將式(9)～(12)應用到耦合系統(tǒng)中。最終時間積分格式可以表示為

(13)

(14)

2 幾何模型及邊界條件設置

2.1 幾何模型及相關參數

圖1示出二維噴氣織機主噴嘴與纖維的幾何模型。模型中黑色的矩形代表了纖維所占據的區(qū)域，并將纖維的初始位置固定于噴嘴的中心軸線處，纖維的左端受到人為的約束作用，自由度為0；纖維的右端不受任何約束，自由度為y方向以及z方向的2個自由度。圖1還示出噴嘴模型的壓力入口、壓力出口以及流固耦合壁面位置。氣流場與纖維主要相關參數如表1、2所示。

圖1 幾何模型及其邊界條件Fig.1 Geometric model and boundary conditions

入口條件/MPa出口條件/MPa單元大小/mm時間步長網格數量0.20.12×10-40.0011.1736×105

表2 纖維主要相關參數

2.2 纖維/壁面接觸函數定義

纖維在高速氣流中會產生大變形和自由運動，且在運動過程中有可能與壁面發(fā)生接觸。纖維/壁面接觸的動力方程[14]為

圖2所示為纖維與壁面的接觸，暫不考慮纖維與壁面之間的摩擦力。

設SXB為纖維與壁面接觸的作用面，SBX為壁面與纖維接觸的作用面，Sc=SXB∩SBX為t時刻纖維與壁面的接觸面。定義接觸函數為

(16)

將f表示交界面上的法向接觸壓力，則上述接觸作用應滿足的法向條件：

圖2 纖維與壁面接觸系統(tǒng)Fig.2 System of fiber and wall in contact

g≥0,f≥0,gf=0

(17)

如果g>0，則表示纖維與壁面沒有接觸，則f=0；如果g=0，則表示纖維與壁面之間存在接觸，則f>0。

在模擬纖維與氣流耦合的時候，由于接觸函數的存在，當纖維接觸壁面的時候，不會因為纖維與壁面的接觸而停止計算。由于接觸時f>0，壁面對纖維的作用力會使得纖維往反方向運動。

3 數值結果與分析

3.1 纖維的運動變形特性

對于纖維的彎曲變形特性，從宏觀物理學上分析，即為纖維受到外力作用而產生體積或形狀的改變；從微觀定量上分析，即為纖維邊界點偏離水平軸線的位移大小。研究中選取了6個時刻的纖維彎曲變形特點，如圖3所示。

圖3 纖維在6個時刻的彎曲變形情況Fig.3 Bending deformation of fiber at six moments

由圖3可以看到，纖維在初始段內，因為沒有受到氣流的作用，保持著水平狀態(tài)。但是，隨著時間的增加，纖維開始出現彎曲變形。纖維首先從約束部分開始變形，然后延伸到纖維的中間部位，然后自由端開始接觸噴嘴導紗管的下壁面，隨著氣壓不斷變化，纖維的自由端碰到下壁面后有一個向上變形的趨勢，最后觸碰到導紗管上壁面。纖維在t=1.582 ms時，纖維基本恢復到與初始相一致的狀態(tài)，完成1個周期的運動。

為了定量分析纖維的波動變形情況，在纖維的邊界上均勻地取17個坐標點，并且在以上6個時刻中取了3個時刻來分析這17個坐標點在z坐標軸上的位移，通過這17個點相對噴嘴中心軸線的距離來定量地分析纖維的彎曲變形情況。彎曲變形特性如圖4所示。

圖4 纖維在3個時刻的波動情況Fig.4 Volatility of fiber at three moments

圖4中每條曲線代表了不同時間點上纖維相對噴嘴中心軸線的距離。在最初的時間點上，纖維并沒有明顯離開噴嘴中心軸線的趨勢，圖上正方體模型曲線就代表了纖維最初的波動情況。

纖維左端被約束的區(qū)域附近并沒有發(fā)生非常明顯的彎曲變形。變形最大的區(qū)域是接近自由端的區(qū)域，從圖4可看出，在這3個時刻纖維每個部位基本都存在彎曲變形，并且呈現三角正弦或者余弦波形式的波動。

纖維不僅在高速氣流中會產生彎曲變形，也會產生拉伸和壓縮變形。對纖維進行離散分析時將纖維表面用相同大小的等長單元尺寸來進行離散。在纖維約束端截取了一小段長度，一端受到約束條件的限制無法運動，另一個點可以在Y-Z平面上任意運動，截取的纖維段如圖5所示。

圖5 纖維段上的約束點與自由點Fig.5 Constrained point and free point of fiber

分析自由點在Y軸上的運動情況，并且以約束點作為參考點，約束點任何時刻在Y方向上的位移都為0，而自由點不受約束，它在Y軸上的位移會隨著時間的變化而發(fā)生變化。圖6示出自由點和約束點在Y方向上隨時間變化的情況。

圖6 自由點和約束點在Y方向上的位移變化量Fig.6 Displacement variation in Y-direction of free point and constraint point over time

圖6中虛線為約束點在Y軸上的位移變化量，實線為自由點在Y軸上的位移變化量。由于受到約束作用，約束點在任何時間點上都不發(fā)生任何位移變化，在圖中表現為一條水平直線，即位移變化量為0；而對于自由點，在剪切力的作用下就會在Y軸上發(fā)生位移變化，正值代表伸長量，負值代表壓縮量。顯然相對圖4纖維在Z方向上的位移變化量就可以看出，纖維的拉伸壓縮量非常小。

3.2 纖維的受力特性

纖維受到的彎曲影響主要與纖維表面受到的主應力有關，而纖維受到的拉伸變形作用與纖維表面受到的切應力有關；因此，分析纖維在3個時刻上受到的主應力與切應力的分布。結果如表3、4所示。

表3 3個時刻主應力在纖維上的分布

表4 3個時刻切應力在纖維上的分布

Tab.4 Distribution of shear stress in fiber at three moments

纖維長度/mm切應力/MPat=0.625mst=0.79mst=1.012ms00.202.351.453.20.283.402.096.40.035.532.279.60.254.042.9412.60.252.033.5015.80.146.653.7718.80.089.553.2122.00.030.433.43

從表3、4可以看出，在t=0.625 ms的時間上，受到的最大應力值都集中在纖維的中間部位，這與圖3中t=0.625 ms時纖維的彎曲集中在中間位置符合一致，說明纖維在彎曲時受到的應力是最大的。在t=0.79 ms時，纖維受到2個應力的最大值均位于18.8 mm附近，而在圖3中纖維在該區(qū)域正與壁面接觸，說明纖維在與壁面接觸時受到的應力也比較大，并且由于壁面對纖維的反作用力會使得纖維往反方向運動。在t=1.012 ms時，與圖3比較可以分析纖維前半部分并沒有發(fā)生明顯的彎曲變形，該部分受到的應力相對穩(wěn)定。

4 結 論

1) 修正纖維的本構方程可以得到纖維應力與應變對應的函數關系；在Arbitrary-Lagrangian-Eulerian坐標系下構建的纖維與氣流耦合控制方程可以數值模擬纖維的彎曲變形。

2) 纖維在氣流的作用下會產生彎曲、拉伸變形，約束端的纖維在氣流中的彎曲變形不太明顯；而自由端的纖維會產生明顯的彎曲變形特性，并且呈現正弦或者余弦形式的波動。

3) 通過對纖維的受力分析可知，纖維受到應力最大的部分往往集中在纖維的彎曲部位和纖維與壁面接觸部位。

4) 通過研究纖維的彎曲變形特性以及受力情況將有助于進一步研究流體動力型紡織機械的機制性。

[1] SMITH A C, ROBERTS W W. Straightening of crimped and hooked fibers in converging transport ducts: computational modeling[J]. Journal of Textile Research, 1994, 64(6): 335-344.

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

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

[4] LINDSTRTIM SB, UESAKA T. Simulation of the motion of flexible fibers in viscous fluid flow[J]. Physics of Fluids, 2007, 19(11): 113307-1-16.

[5] TAKEMURA M, CHIBA K, NAKAMURA K. Motion of flexible fibers in a newtonian flow: part 2:evolution of the configuration of a single fiber in a simple shear flow[J]. Journal of Textile Engineering, 2001, 47(3): 77-91.

[6] TORNBERG AK, SHELLEY MJ. Simulating the dynamics and interactions of flexible fibers in Stokes flows[J]. Journal of Computational Physics, 2004,196(1): 8- 40.

[7] DONG S, FENG X, SALCUDEAN M, et al. Concentration of pulp fibers in 3D turbulent channel flow[J]. International Journal of Multiphase Flow, 2003, 29: 1-21.

[8] PEI Zeguang, YU Chongwen. Numerical study on the effect of nozzle pressure and yarn delivery speed on the fiber motion in the nozzle of Murata vortex spinning[J]. Journal of Fluids and Structures, 2011,27:121-133.

[9] 徐存強，馮志華，董騰中,等.氣流引緯主噴嘴內氣流速度與緯紗間摩擦因數的關系[J].紡織學報， 2013,34(11):147-150. XU Cunqiang, FENG Zhihua, DONG Tengzhong, et al. The relationship between friction coefficients of the weft yarns and velocity of the airflow in main nozzle during weft insertion[J]. Journal of Textile Research, 2013, 34(11): 147-150.

[10] 李翠玉,張義同,徐家福.機織物懸垂屈曲的數值模擬[J]. 天津大學學報，2007,40(1):46-50. LI Cuiyu, ZHANG Yitong, XU Jiafu. Numerical simulation of draping and buckling of woven fabric[J]. Journal of Tianjin University, 2007, 40(1):46-49.

[11] 汪靜，金玉珍，胡旭東,等.基于Win CE的噴氣織機人機交互系統(tǒng)設計[J]. 紡織學報，2010,31(7):122-126. WANG Jing, JIN Yuzhen, HU Xudong, et al. Win CE based design of human-machine interaction system for air-jet loom[J]. Journal of Textile Research, 2010, 31(7): 122-126.

[12] ZHANG Hou, KLAUS-JURGEN Bathe. Direct and iterative computing of fluid flows fully coupled with structures[J]. Computational Fluid and Solid Mechanics, 2001(1): 1440-1443.

[13] BATHE Klaus-Jurgen, ZHANG Hou. Finite element developments for general fluid flows with structural interactions[J]. International Journal for Numerical Methods in Engineering, 2004, 60:213-232.

[14] PANTUSO D, BATHE K J, BOUZINOV P A. A finite element procedure for the analysis of thermo-mechanical solids in contact[J]. Computers and Structures, 2000, 75: 551-573.

Numerical simulation on coupling of fiber and air flow

JIN Yuzhen, LI Jun, LIN Peifeng, WU Zhenyu

(FacultyofMechanicalEngineeringandAutomation,ZhejiangSci-TechUniversity,Hangzhou,Zhejiang310018,China)

A two-dimensional model of fluid-structure coupling was set up based on the equation of the fiber and the fluid control. The consistent time integration and iterative coupling algorithm were used to solve the equation. The interaction of the fiber and the airflow in the main nozzle, such as the fiber motional, deformational characteristics and the force of the fiber, was analyzed. The results show that the fiber under airflow will come into stretching and bending deformation and present in the form of sine or cosine. The deformation effect of the free end of the fiber is more obvious than that of the restrained end. The maximum stress often occurred in the bending position as well as the part where the fiber contacted with the wall.

fluid-structure interaction; numerical simulation; fiber; adaptive mesh

10.13475/j.fzxb.201501015206

2013-12-02

2014-09-30

浙江省自然科學基金資助項目(LZ14E050004, LQ12A02002)；浙江理工大學流體工程技術創(chuàng)新團隊資助項目(11132932611309)；浙江理工大學研究生創(chuàng)新研究項目(YCX13023)

金玉珍(1979—)，女，副教授，博士。研究方向為紡織裝備流體傳動。E-mail: gracia1101@foxmail.com。

TS 101.2

A