聚合物壓力拖拽流的有限元分析及數(shù)值計算

王 超

(中南民族大學 計算機科學學院，武漢430074)

聚合物壓力拖拽流的有限元分析及數(shù)值計算

王 超

(中南民族大學 計算機科學學院，武漢430074)

使用基于Navier-Stock方程的牛頓流體模型對壓力拖拽流下的低剪切速率聚合物流體進行了數(shù)學建模，利用有限元法對計算區(qū)域進行網(wǎng)格離散，并采用插值函數(shù)得到了流體連續(xù)性方程與運動方程的權(quán)函數(shù).通過設定網(wǎng)格結(jié)點數(shù)據(jù)、結(jié)點坐標及相關(guān)邊界條件，建立了流變模型的有限元總體方程.其計算結(jié)果通過Tecplot軟件后處理得到平板間流體流動的速度及壓力分布，數(shù)值計算結(jié)果表明：速度分布在復合條件作用下與在兩種邊界條件單獨作用下完全相等，且壓力分布也符合邊界條件設定，證明了其求解的正確性；通過更改網(wǎng)格劃分及設定方式，此方法亦可用于其他形狀規(guī)則的牛頓流體二維流變模型有限元分析.

Navier-Stock方程；牛頓流體模型；二維流變模型；有限元分析

計算流體力學目前已廣泛應用于航空航天、能源和動力工程、力學、物理和化學、建筑、水利、環(huán)境等領域.高分子流變學作為計算流體力學的一個重要分支[1-3]，其研究成果對高分子凝聚態(tài)物理、乃至高分子科學基礎理論發(fā)展具有重要價值；高分子材料在成型加工中流變狀態(tài)的變化不僅決定制品外觀形狀和質(zhì)量，而且對分子鏈結(jié)構(gòu)的形成和變化有極其重要的影響，是決定高分子制品最終結(jié)構(gòu)、性能的中心環(huán)節(jié)[4,5].

本文以聚合物平行板間牛頓流體壓力拖拽為研究對象，對Navier-Stock方程組進行速度-壓力有限元求解，對平行板間聚合物流體的二維流動模型進行有限元分析，利用數(shù)值計算，得到牛頓流體流變模型下聚合物流體速度-壓力的流動分布[6]，并比較多種邊界條件復合作用下流體的流動特性.

1 實例模型及數(shù)學方程

1.1 幾何模型

以平行板空間內(nèi)LLDPE材料緩慢定常流動為例，將其視為牛頓流體.平行板上表面拖拽速度為u=0.02m/s，入口壓力為P=2000Pa，出口壓力為P=0，平行板下壁面固定.流體黏度取值為LLDPE在230℃時黏度為μ=1300Pa·s.圖1即為牛頓流體在壓力拖拽復合作用下二維流場分布的簡化模型.

圖1 牛頓流體平板間二維流動模型Fig.1 Two dimensional flow model of Newton fluid in flat plate

1.2 數(shù)學方程

考慮低剪切速率條件下，平行板間中聚合物材料流動可近似認為理想牛頓流體的流動，在牛頓流體模型中，要保持穩(wěn)定的流動，所需的剪切應力與剪切速率成正比[7]，即：

(1)

在動量守恒定律的基礎上可建立流體所承受的力和流體運動之間的關(guān)系，上述流體模型的運動方程可以用Navier-Stock方程組來描述：

(2)

(3)

(4)

其中，剪切應力滿足牛頓本構(gòu)方程:

(5)

τxx為作用在x面(與x軸垂直的yz面)上，作用方向為x方向的剪切應力分量；τxy為作用在x面(與x軸垂直的yz面)上，作用方向為y方向的剪切應力分量；τyx為作用在y面(與y軸垂直的xz面)上，作用方向為x方向的剪切應力分量；τyy為作用在y面(與y軸垂直的xz面)上，作用方向為y方向的剪切應力分量.

由于劃分網(wǎng)格進行累加計算時，需進行積分數(shù)值計算，利用二維積分高斯公式可得：

(6)

其中Hi和Hj為積分權(quán)；ξi和ηj為積分點.其數(shù)值查詢高斯求積公式積分點坐標與加權(quán)系數(shù)表可得到.

2 有限元求解及數(shù)值計算

2.1 計算區(qū)域的網(wǎng)格劃分

計算區(qū)域的網(wǎng)格劃分采用四邊形單元[8]，整個計算域的全局網(wǎng)格劃分如圖2所示，將整個單元先劃分為12個四邊形單元，再將單個四邊形由各邊中點二次劃分為四個小區(qū)域，即對應的速度單元比壓力單元高一個階次.整個網(wǎng)格離散單元數(shù)為12，速度單元節(jié)點總數(shù)為63，壓力節(jié)點總數(shù)為20.

圖2 計算區(qū)域的網(wǎng)格離散Fig.2 Grid discretization of the computational region

根據(jù)網(wǎng)格離散劃分得到各參數(shù)表(表1～4)，包括速度單元結(jié)點數(shù)據(jù)JMV、壓力單元結(jié)點數(shù)據(jù)JMP、速度單元結(jié)點坐標數(shù)據(jù)JXYV、壓力單元結(jié)點坐標數(shù)據(jù)JXYP、速度邊界條件JBV、壓力邊界條件JBP.JMV和JMP分別存儲速度單元和壓力單元所包含的結(jié)點序號.JXYV與JXYP分別存儲各個速度與壓力各結(jié)點的坐標數(shù)據(jù).JBV存儲邊界結(jié)點號及相應x和y方向的速度.JBP存儲處于壓力邊界條件上的單元號、邊界單元的邊界邊號、邊界邊的外法線方向余弦cosθx、邊界邊的外法線方向余弦cosθy、邊界邊第一點處壓力值、邊界邊第二點處壓力值.其中因網(wǎng)格劃分為正方向，邊界邊的外法線與邊界邊夾角均為π/2的整數(shù)倍，故取值在-1、0和1之中.主程序根據(jù)各參數(shù)離散數(shù)據(jù)進行矩陣運算.

表1 速度單元結(jié)點數(shù)據(jù)JMV

表2 壓力單元結(jié)點數(shù)據(jù)JMV

表3 速度邊界條件JBV

表4 壓力邊界條件JBP

2.2 插值函數(shù)與單元方程的建立

離散單元數(shù)據(jù)的處理應引入插值函數(shù)[9,10]，φ為速度單元插值函數(shù)，表示為方程組：

(7)

ψ為壓力單元插值函數(shù)，可表示為方程組：

(8)

式中，-1≤ξ≤1,-1≤η≤1.

插值函數(shù)φ對x、y的導數(shù)與對ξ、η的導數(shù)之間存在以下關(guān)系(J為雅克比矩陣)：

(9)

將插值權(quán)函數(shù)ψ與連續(xù)性方程相乘組成連續(xù)性方程的加權(quán)余量方程，即：

(10)

將(10)式由整個求解區(qū)域轉(zhuǎn)換為網(wǎng)格劃分的單元區(qū)域.

(11)

又由插值函數(shù)的性質(zhì)可知單元內(nèi)任一點的速度與壓力數(shù)值可表示為：

(12)

將方程組(12)代入(11)并展開，得：

(13)

其簡化形式為:

(14)

(15)

(16)

(17)

(18)

將式(4)代入式(3)并與權(quán)函數(shù)φ相乘，在計算區(qū)域內(nèi)積分，整理得:

(19)

(20)

將各部分分別展開，導入插值函數(shù)并轉(zhuǎn)換到單元區(qū)域可得:

(21)

(22)

以上兩式簡化形式為：

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

綜合上式可得網(wǎng)格離散單元的總體方程為:

(33)

通過迭代累加可組成總體方程，其中B、C、D、F是由對應的單元矩陣合成:

(34)

3 計算結(jié)果與后處理

3.1 計算流程

整個有限元數(shù)值計算的求解流程如圖3所示.

圖3 總體求解流程Fig.3 General solving process

3.2 數(shù)據(jù)后處理

得到總體方程的解后，將各結(jié)點坐標、x與y方向結(jié)點速度，結(jié)點壓力等數(shù)據(jù)導入Tecplot進行后處理，采用無結(jié)構(gòu)數(shù)據(jù)進行處理與顯示，如圖4所示，圖4a為復合邊界作用下速度分布圖，圖4b為上壁面拖拽單獨作用下的速度分布，圖4c為入口壓力單獨作用下的速度分布.圖5為復合邊界作用下壓力等值線.

a)復合邊界作用下速度分布圖；

b)上壁面拖拽作用下速度分布;

c)入口壓力作用下速度分布圖4 流場速度分布與壓力分布圖Fig.4 Velocity distribution and pressure distribution of flow field

圖5 復合邊界作用下壓力等值線Fig.5 Pressure contour line under complex boundary conditions

因速度等值線平行于平板，在平板任一點上作垂直于平板的垂線，則可得到平板中所有結(jié)點的速度值，表5為以平板中點為垂線選取各結(jié)點的速度數(shù)值，由以上數(shù)據(jù)可知，上壁面拖拽單獨作用與入口壓力單獨作用數(shù)據(jù)的累加結(jié)果等于兩種邊界條件復合作用下的速度分布數(shù)值.流體的性質(zhì)復合牛頓流體的累加特性，且由圖5可知，平板間壓力分布從入口開始逐漸減小，垂直于平板方向的壓力等值線數(shù)值無變化，復合邊界條件的設定.

表5 單一邊界作用疊加與復合邊界作用下平板中線結(jié)點的速度對比

4 結(jié)果與討論

本文使用速度-壓力有限元法求解聚合物流體二維流動模型，主要針對形狀規(guī)則的牛頓流體流動情況.使用四邊形單元進行網(wǎng)格離散，根據(jù)Navier-Stock連續(xù)性方程，引入插值函數(shù)與雅克比矩陣完成加權(quán)余量方程與矩陣單元方程的建立，并組合成總體方程求解得到最終的速度、壓力數(shù)值.該方法同樣適用于大多數(shù)規(guī)則形狀中流體的流動，需在網(wǎng)格離散程序結(jié)點拓撲分布中生成相應的求解區(qū)域及結(jié)點信息.同樣可增加求解的網(wǎng)格數(shù)量，但對于牛頓流體的求解不存在迭代計算，因此對數(shù)據(jù)最終精度的影響不大.

[1] Reddy J N, Gartling D K. The finite element method in heat transfer and fluid dynamics[M]. 2nd. Boca Raton: CRC Press Inc, 2001.

[2] 宋維寧,林 祥, 任冬云. 一種廣義牛頓流體的廣義黏度模型[J]. 中國塑料，2009, 23(12): 26-30.

[3] 楊曉東, 劉保臣, 劉春太，等. 高分子聚合物熔體Cross黏度模型的改進[J]. 高分子材料科學與工程，2010, 26(11): 172-174.

[4] 張任良, 狄勤豐, 王新亮，等. 基于Shan-Chen模型的格子Boltzmann方法在微流動模擬研究中的應用[J]. 力學與實踐, 2012, 34(2): 10-18.

[5] Ohta M, Nakamura T, Yoshida Y, et al. Lattice Boltzmann simulations of viscoplastic fluid flows through complex flow channels[J]. Journal of Non-Newtonian Fluid Mechanics, 2011, 166(7): 404-412.

[6] 張 敏, 孫 勝, 賈玉璽，等. 聚合物共擠出過程的擠出脹大有限元分析[J]. 高分子材料科學與工程，2006, 22(5): 36-40.

[7] 徐佩弦. 高聚物流變學及其應用[M]. 北京：化學工業(yè)出版社，2003.

[8] 畢 超. 計算流體力學有限元方法及其編程詳解[M]. 北京：機械工業(yè)出版社, 2013.

[9] Fallah K,Khayat M,Borghei MH,et al.Multiplerelaxation-time lattice Boltzmann simulation of non-Newtonian flows past a rotating circular cylinder[J]. Journal of Non-Newtonian Fluid Mechanics, 2012, 177-178(1): 1-14.

[10] Chai ZH, Shia BC, Guo ZL,et al. Multiple-relaxation-time lattice Boltzmann model for generalized Newtonian fluid flows[J]. Journal of Non-Newtonian Fluid Mechanics, 2011, 166(5): 332-342.

The Research on Finite Element Analysis and Numerical Calculation of Pressure Drag Stream of Polymer

Wang Chao

(College of Computer Science, South-Central University for Nationalities, Wuhan 430074, China)

Using Newton fluid model based on Navier-Stock equation, a mathematical model of the low shear rate polymer fluid under pressure drag and drop flow was carried out. The finite element method is used to calculate the grid discretization, and the interpolation function is used to obtain the weight function of fluid continuity equation and fluid equation of motion. The finite element equation of the rheological model is established by setting the grid node data, the node coordinates and the related boundary conditions. The velocity and pressure distribution of fluid flow in the flat plate were obtained by Tecplot software. The numerical results show that the velocity distribution under the composite condition is equal to that under separate action of two boundary conditions, the pressure distribution is also consistent with the boundary conditions, which proves the correctness of the solution. By changing the meshing and setting method, this method can also be used in the finite element analysis of two dimensional rheological model of other shape rules of Newton fluid.

Navier-Stock equation；Newton fluid model；two dimensional rheological model；finite element analysis

2016-03-15

王 超(1980-)，男，講師，博士，研究方向：聚合物成型與控制，E-mail: wangchaofly@126.com

國家自然科學青年基金資助項目(51403160)

TQ315

A

1672-4321(2016)04-0070-06