纖維與氣流耦合的數(shù)值模擬

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

(浙江理工大學(xué) 機(jī)械與自動(dòng)控制學(xué)院， 浙江 杭州 310018)

?

纖維與氣流耦合的數(shù)值模擬

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

(浙江理工大學(xué) 機(jī)械與自動(dòng)控制學(xué)院， 浙江 杭州 310018)

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

流固耦合； 數(shù)值模擬； 纖維； 自適應(yīng)網(wǎng)格

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

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

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

1 數(shù)學(xué)模型

1.1 纖維本構(gòu)方程

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

(1)

式中：σ11、σ22表示纖維在二維坐標(biāo)軸上的2個(gè)主應(yīng)力；相應(yīng)的ε11、ε22為正應(yīng)變；σ12代表纖維在二維坐標(biāo)軸上的切向應(yīng)力；相應(yīng)的γ12為切向應(yīng)變；kij為材料的剛度。材料的剛度可以通過材料的泊松比以及彈性模量等參數(shù)來調(diào)整。通過該本構(gòu)方程就能夠得到纖維應(yīng)力與應(yīng)變對(duì)應(yīng)的函數(shù)關(guān)系。

1.2 纖維與氣流控制方程

在二維直角坐標(biāo)系下，不考慮密度變化以及熱交換可以得到流體運(yùn)動(dòng)控制方程[11]：

(2)

(3)

(4)

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

纖維的運(yùn)動(dòng)方程為

(5)

(6)

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

運(yùn)動(dòng)學(xué)[12]條件為

(7)

動(dòng)力學(xué)條件為

(8)

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

1.3 流固耦合相容時(shí)間積分

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

(9)

(10)

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

(11)

(12)

式中，ξt、ηt分別為t時(shí)刻有關(guān)速度和位移的變化量。

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

(13)

(14)

2 幾何模型及邊界條件設(shè)置

2.1 幾何模型及相關(guān)參數(shù)

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

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

入口條件/MPa出口條件/MPa單元大小/mm時(shí)間步長網(wǎng)格數(shù)量0.20.12×10-40.0011.1736×105

表2 纖維主要相關(guān)參數(shù)

2.2 纖維/壁面接觸函數(shù)定義

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

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

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

(16)

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

圖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。

在模擬纖維與氣流耦合的時(shí)候，由于接觸函數(shù)的存在，當(dāng)纖維接觸壁面的時(shí)候，不會(huì)因?yàn)槔w維與壁面的接觸而停止計(jì)算。由于接觸時(shí)f>0，壁面對(duì)纖維的作用力會(huì)使得纖維往反方向運(yùn)動(dòng)。

3 數(shù)值結(jié)果與分析

3.1 纖維的運(yùn)動(dòng)變形特性

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

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

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

為了定量分析纖維的波動(dòng)變形情況，在纖維的邊界上均勻地取17個(gè)坐標(biāo)點(diǎn)，并且在以上6個(gè)時(shí)刻中取了3個(gè)時(shí)刻來分析這17個(gè)坐標(biāo)點(diǎn)在z坐標(biāo)軸上的位移，通過這17個(gè)點(diǎn)相對(duì)噴嘴中心軸線的距離來定量地分析纖維的彎曲變形情況。彎曲變形特性如圖4所示。

圖4 纖維在3個(gè)時(shí)刻的波動(dòng)情況Fig.4 Volatility of fiber at three moments

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

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

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

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

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

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

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

3.2 纖維的受力特性

纖維受到的彎曲影響主要與纖維表面受到的主應(yīng)力有關(guān)，而纖維受到的拉伸變形作用與纖維表面受到的切應(yīng)力有關(guān)；因此，分析纖維在3個(gè)時(shí)刻上受到的主應(yīng)力與切應(yīng)力的分布。結(jié)果如表3、4所示。

表3 3個(gè)時(shí)刻主應(yīng)力在纖維上的分布

表4 3個(gè)時(shí)刻切應(yīng)力在纖維上的分布

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

纖維長度/mm切應(yīng)力/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的時(shí)間上，受到的最大應(yīng)力值都集中在纖維的中間部位，這與圖3中t=0.625 ms時(shí)纖維的彎曲集中在中間位置符合一致，說明纖維在彎曲時(shí)受到的應(yīng)力是最大的。在t=0.79 ms時(shí)，纖維受到2個(gè)應(yīng)力的最大值均位于18.8 mm附近，而在圖3中纖維在該區(qū)域正與壁面接觸，說明纖維在與壁面接觸時(shí)受到的應(yīng)力也比較大，并且由于壁面對(duì)纖維的反作用力會(huì)使得纖維往反方向運(yùn)動(dòng)。在t=1.012 ms時(shí)，與圖3比較可以分析纖維前半部分并沒有發(fā)生明顯的彎曲變形，該部分受到的應(yīng)力相對(duì)穩(wěn)定。

4 結(jié) 論

1) 修正纖維的本構(gòu)方程可以得到纖維應(yīng)力與應(yīng)變對(duì)應(yīng)的函數(shù)關(guān)系；在Arbitrary-Lagrangian-Eulerian坐標(biāo)系下構(gòu)建的纖維與氣流耦合控制方程可以數(shù)值模擬纖維的彎曲變形。

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

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

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

[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] 徐存強(qiáng)，馮志華，董騰中,等.氣流引緯主噴嘴內(nèi)氣流速度與緯紗間摩擦因數(shù)的關(guān)系[J].紡織學(xué)報(bào)， 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ī)織物懸垂屈曲的數(shù)值模擬[J]. 天津大學(xué)學(xué)報(bào)，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的噴氣織機(jī)人機(jī)交互系統(tǒng)設(shè)計(jì)[J]. 紡織學(xué)報(bào)，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

浙江省自然科學(xué)基金資助項(xiàng)目(LZ14E050004, LQ12A02002)；浙江理工大學(xué)流體工程技術(shù)創(chuàng)新團(tuán)隊(duì)資助項(xiàng)目(11132932611309)；浙江理工大學(xué)研究生創(chuàng)新研究項(xiàng)目(YCX13023)

金玉珍(1979—)，女，副教授，博士。研究方向?yàn)榧徔椦b備流體傳動(dòng)。E-mail: gracia1101@foxmail.com。

TS 101.2

A