劉 紅

(安徽理工大學(xué)理學(xué)院，安徽 淮南 232001)

一種基于離散元法的彎管中顆粒流動(dòng)問(wèn)題的數(shù)值模擬

劉 紅

(安徽理工大學(xué)理學(xué)院，安徽 淮南 232001)

本文建立了直角彎管中顆粒流動(dòng)問(wèn)題的數(shù)學(xué)模型，并對(duì)核心模塊的算法進(jìn)行詳細(xì)描述，在此基礎(chǔ)上設(shè)計(jì)了基于離散元法求解彎管中顆粒流動(dòng)問(wèn)題的算法與程序。數(shù)值實(shí)驗(yàn)結(jié)果表明：實(shí)驗(yàn)結(jié)果與實(shí)際物理現(xiàn)象相符，驗(yàn)證了所設(shè)計(jì)的求解算法和所研制的程序的正確性。

離散元法；直角彎管；顆粒

1 引言

管道運(yùn)輸具有一些其他運(yùn)輸機(jī)械無(wú)法比擬的優(yōu)勢(shì),它在工農(nóng)業(yè)和物流運(yùn)輸業(yè)中應(yīng)用非常廣泛[1],目前在高速運(yùn)輸?shù)臄?shù)值模擬方面和管壁磨損方面的研究是較多的, 對(duì)于顆粒在運(yùn)動(dòng)過(guò)程中對(duì)彎管的磨損機(jī)理研究也逐漸被人們重視[2]。彎管中顆粒的流動(dòng)增加了問(wèn)題處理的復(fù)雜度,因此多數(shù)選取球形顆粒(二維選取圓盤(pán)形顆粒)來(lái)模擬。

近年來(lái),計(jì)算機(jī)軟件和硬件的快速發(fā)展[3],為長(zhǎng)時(shí)間大規(guī)模的數(shù)值模擬打下了堅(jiān)實(shí)的基礎(chǔ)。離散單元法在研究氣-固兩相流中固相方面有獨(dú)特的優(yōu)勢(shì), 離散元這一概念最早是1971[4]年Cundall為了分析和處理巖土工程方面的問(wèn)題而提出來(lái)的,對(duì)于不連續(xù)介質(zhì)，離散元法是行之有效的數(shù)值模擬方法[5]。在實(shí)際應(yīng)用方面,離散元技術(shù)在巖土、化工、農(nóng)業(yè)、風(fēng)沙流動(dòng)及散粒材料的運(yùn)輸?shù)阮I(lǐng)域有廣泛應(yīng)用。1986年王泳嘉[6]首次在會(huì)議上向國(guó)內(nèi)巖石力學(xué)與工程界介紹了離散元法。李世海[7]在2012年顆粒材料計(jì)算力學(xué)會(huì)議上發(fā)表了CDEM的最新研究進(jìn)展。Zhang[8-9]等對(duì)流體輸送過(guò)程用離散元方法進(jìn)行了數(shù)值模擬,分析了顆粒之間的相互作用和顆粒對(duì)管道的沖蝕磨損作用,并側(cè)重預(yù)測(cè)輸送管道最容易損壞的位置,對(duì)實(shí)際的流體輸送工程問(wèn)題有很大的參考價(jià)值。Yue等[10]基于CPU-GPU異構(gòu)架構(gòu)將DEM的串行代碼做了并行化處理，提升了其的計(jì)算效率。毛靖儒等[11]采用離散元法對(duì)彎管內(nèi)稀疏氣固兩相流進(jìn)行了模擬，同時(shí)他也研究了固體顆粒對(duì)壁面磨損效應(yīng),指出了壁面處粒子的分布、碰撞角度及濃度分布是固粒對(duì)壁面磨損的主要因素。

本文首先對(duì)Hertz模型和M-D模型進(jìn)行詳細(xì)的總結(jié)并給出了平動(dòng)、轉(zhuǎn)動(dòng)平衡方程的離散格式，然后給出模型的基本假設(shè)，并建立直角彎管中顆粒流動(dòng)問(wèn)題的模型，給出了基于離散元法的核心模塊的算法,最后用直角彎管中顆粒流動(dòng)的數(shù)值實(shí)驗(yàn)對(duì)我們的程序設(shè)計(jì)進(jìn)行了正確性驗(yàn)證，得出結(jié)論。

2 離散單元法

考慮到運(yùn)動(dòng)顆粒之間的碰撞并使用相互作用方程，顆粒受到系統(tǒng)中其它顆粒對(duì)其的接觸作用,遵循質(zhì)心運(yùn)動(dòng)定理且滿足牛頓第二運(yùn)動(dòng)定律，則有

(1)

(2)

(3)

其中E1(E2)和,μ1(μ2)分別表示顆粒的彈性模量和泊松比。

粒間的法向接觸力增量ΔF與兩顆粒之間的重疊量Δα關(guān)系如下:

(4)

在實(shí)際的接觸問(wèn)題中,顆粒間存在著相對(duì)滑動(dòng)和摩擦作用,且這些作用是不可以忽略的,因此Hertz模型不能夠準(zhǔn)確地描述接觸力,由于相互接觸產(chǎn)生的切向接觸力的作用,使得兩顆粒在切向力方向上發(fā)生相對(duì)滑移,滑移從接觸表面開(kāi)始之后蔓延到顆粒的內(nèi)部,Mindlin-Deresiewic模型考慮了切向力與切向位移的聯(lián)系,可以更加真實(shí)地描述這一過(guò)程。切向力增量為

(5)

其中，δ，Δδ和T分別表示兩個(gè)顆粒接觸表面切向位移，切向位移增量和切向力。通常采用向前差分(中心差分)格式對(duì)方程左端的導(dǎo)數(shù)項(xiàng)進(jìn)行離散，用顯示格式對(duì)平動(dòng)和轉(zhuǎn)動(dòng)平衡方程的進(jìn)行求解 ,因此需要選取合適的時(shí)間步長(zhǎng)使得兩顆粒間的作用力在Δt的時(shí)間內(nèi)變化非常小，這樣計(jì)算的穩(wěn)定性和準(zhǔn)確性才能得到保證。通常我們由Rayleigh波的波速來(lái)確定時(shí)間步長(zhǎng)。假設(shè)顆粒的材質(zhì)為彈性固體, ρ,G，ν，R分別表示顆粒的密度，剪切模量、泊松比、不同尺寸的顆粒組成的系統(tǒng),時(shí)間步長(zhǎng)由下式計(jì)算得出

(6)

式(6)是在顆粒間相對(duì)速度較小時(shí)計(jì)算得出的,因此這種確定時(shí)間步長(zhǎng)的方法比較適用于顆粒運(yùn)動(dòng)不太劇烈的系統(tǒng)。依據(jù)實(shí)際情況選取合適的時(shí)間步長(zhǎng)來(lái)確保穩(wěn)定性,可以在時(shí)間步長(zhǎng)前乘一個(gè)常數(shù)C (0

3 直角彎管中顆粒流動(dòng)模型的建立與算法描述

3.1 物理模型的建立和基本參數(shù)

二維平面內(nèi)的直角彎管,實(shí)際上也就是同心圓環(huán)的四分之一圓周,為此我們將圓環(huán)進(jìn)行剖分,采用一小段線段代替圓周上所對(duì)應(yīng)的一段圓弧。根據(jù)二維平面圓的參數(shù)方程(7)即可得到每一面墻的兩個(gè)端點(diǎn)坐標(biāo)。

(7)

用球形顆粒代表真實(shí)顆粒,用墻表示管道壁,建立離散元模型模擬了直角彎管中顆粒在重力作用下的流動(dòng),表1為本次數(shù)值模擬的參數(shù)，圖1表示初始狀態(tài).

表1 顆粒和墻的物理參數(shù)

圖1 初始狀態(tài)圖

對(duì)彎管顆粒流動(dòng)的離散元模擬,首先要制備基本穩(wěn)定的顆粒體系。在顆粒生成的團(tuán)塊區(qū)域隨機(jī)的生成1000個(gè)顆粒,并讓顆粒在重力作用下自然沉降,達(dá)到穩(wěn)定狀態(tài),如圖1所示，并將該時(shí)視為零時(shí)刻。下面主要介紹求解過(guò)程中最重要的顆粒運(yùn)動(dòng)模塊，接觸判斷模塊和接觸力求解模塊。

運(yùn)動(dòng)模塊的主要功能為計(jì)算顆粒在每一時(shí)刻的速度與位移。從tn時(shí)刻到tn+1時(shí)刻對(duì)顆粒i的運(yùn)動(dòng)情況可由算法1描述。

算法1：

顆粒間接觸力求解模塊利用前面介紹的接觸模型計(jì)算顆粒的法向接觸力和切向接觸力。算法2為求解顆粒i與顆粒j之間接觸力的算法。

高同型半胱氨酸（Hcy）是半胱氨酸和蛋氨酸代謝的一種重要中間產(chǎn)物，是一種含硫非必需氨基酸，能夠損傷血管內(nèi)皮細(xì)胞，是動(dòng)脈粥樣硬化、血栓、糖尿病腎病的獨(dú)立危險(xiǎn)因素[7-8]。Hcy由腎臟合成與代謝，當(dāng)腎小球?yàn)V過(guò)功能、腎小管代謝功能損害時(shí)，血液Hcy水平升高，可以與高血糖、高糖化血紅蛋白協(xié)同損傷腎臟血管內(nèi)皮細(xì)胞、腎小球基底膜細(xì)胞，使腎小球?yàn)V過(guò)膜孔徑增大，造成腎濾過(guò)功能損害，加重蛋白尿程度。

算法2：

Step 1求顆粒之間的重疊量α：

2) 計(jì)算兩球形顆粒間的重疊量α=ri+rj-d,式中rl為顆粒l(l=i，j)的半徑。

Step 2 計(jì)算法向接觸力fn+1、法向力Fn+1、接觸圓半徑an+1以及法向接觸力增量ΔFn+1。

Step 4求合力以及合力矩：

顆粒之間的法向力是有Hertz模型求得的,其算法描述如下。

算法3：

Step 3 生成法向力大?。篎n+1=fn+1+pn+1。

顆粒間切向力大小是由M-D模型計(jì)算出來(lái)的,其算法描述如下。

算法 4：

Step 3 切向接觸力大?。簍n+1=tn+Δtn+1。

4 直角彎管中顆粒流動(dòng)數(shù)值模擬

打開(kāi)擋板,顆粒在重力作用下,沿著彎管向下流動(dòng),并在最終出口向外流出, 圖2給出流動(dòng)過(guò)程中五個(gè)不同時(shí)刻狀態(tài)圖,其中t1=0.23842s，t2=0.47684s，t3=0.71526s，t4=0.95368s，t5=1.1921s。

圖2 五個(gè)不同時(shí)刻的顆粒流動(dòng)狀態(tài)

如圖2所示,在初始狀態(tài),顆粒由于重力的作用堆積在擋板上面。打開(kāi)擋板之后,顆粒在重力作用下首先向下落,由于在下落過(guò)程中顆粒之間會(huì)相互碰撞，進(jìn)而產(chǎn)生接觸力使得顆粒的運(yùn)動(dòng)速度發(fā)生變化,在這一過(guò)程中有靠近管壁面的部分顆粒也會(huì)與管壁發(fā)生接觸,壁面附近有少量顆粒集中,這就導(dǎo)致了不同顆粒有不同的運(yùn)動(dòng)狀態(tài),所以在豎直管下落的階段,從整體上看顆粒由集密變得稀疏。

當(dāng)顆粒流前端接觸到彎管,顆粒與顆粒、顆粒與彎管壁之間發(fā)生碰撞并出現(xiàn)聚團(tuán)現(xiàn)象。從圖2的第3個(gè)和第4個(gè)圖可以看出,顆粒流進(jìn)入彎管后,由于向心力作用大部分顆粒沿著弧形彎管外壁運(yùn)動(dòng)；顆粒從豎直管進(jìn)入彎管的過(guò)程中,顆粒在與壁面發(fā)生碰撞之前其運(yùn)動(dòng)的方向改變較小,顆粒與壁面接觸后并在向心力的作用下在彎管外壁附近移動(dòng),當(dāng)然也有少部分顆粒受到彎管外壁的彈力作用往彎管內(nèi)壁方向運(yùn)動(dòng),這一過(guò)程中有的顆粒也會(huì)與彎管內(nèi)壁發(fā)生接觸碰撞；圖2上很少顆粒在彎管內(nèi)壁附近運(yùn)動(dòng)。上述結(jié)果說(shuō)明彎管外管壁承是顆粒流的沖擊和磨損作用的主要作用區(qū)域。

顆粒離開(kāi)彎管后進(jìn)入水平管道,由于重力的原因顆粒都集中在下管壁流動(dòng)。可以看出顆粒間的碰撞基本恢復(fù)均勻。最后顆粒離開(kāi)管道做拋物線運(yùn)動(dòng),落在水平地面上，水平速度的不同使得它們的拋物線形狀也有所不同。

5 結(jié) 論

本文首先對(duì)彎管結(jié)構(gòu)進(jìn)行了數(shù)學(xué)建模,并用離散元法模擬了直角彎管中顆粒在重力作用下的流動(dòng)現(xiàn)象，最后的實(shí)驗(yàn)結(jié)果表明設(shè)計(jì)的程序能夠準(zhǔn)確地模擬彎管中顆粒流動(dòng)這一模型。

[1] Stoss HA．Pneumatic Conveying[M]．New York：Wiley Interscience，1983．

[2] 盧洲，劉雪東，潘兵．基于CFD-DEM方法的柱狀顆粒在彎道中輸送過(guò)程的數(shù)值模擬[J]．中國(guó)粉體技術(shù)，2011，17(5)：1-2．

[3] 石鐘慈．第三種科學(xué)方法—計(jì)算機(jī)時(shí)代的科學(xué)計(jì)算[M]．清華大學(xué)出版社，2000．

[4] Cundall P A．A computer model for simulating Progressive large scale movements in blocky R- rock systems[C]// Muller Led．Proc Symp Int Soc Rock Mechanics．Rotterdam：Balkama A A，1971 (1)：1-8．

[5] 王泳嘉，邢紀(jì)波．離散單元法同拉格朗日元法及其在巖土力學(xué)中的應(yīng)用[J]．巖土力學(xué)，1979，16(2)：1-6．

[6] 王泳嘉．離散單元法-一種適用于節(jié)理巖石力學(xué)分析的數(shù)值方法[A]//第一屆全國(guó)巖石力學(xué)數(shù)值計(jì)算機(jī)模型實(shí)驗(yàn)討論會(huì)會(huì)議論文集[C]，1986：32-37．

[7] 李世海，馮春，劉曉宇，等．用于描述地質(zhì)體塊體顆粒的連續(xù)及非連續(xù)計(jì)算模型-CDEM最新進(jìn)展[A]．顆粒材料計(jì)算力學(xué)研究進(jìn)展會(huì)議錄[C]，2012．

[8] Hao Zhang，Yuanqiang Tan，Dongmin Yang，et al． Numerical investigation of the location of maximum erosive wear damage in elbow：Effect of slurry velocity，bend orientation and angle of elbow[J]．Technology．2012，217(3)：467- 476．

[9] Yuanqiang Tan，Hao Zhang．Numerical simulation of concrete pumping process and investigation of wear mechanism of the piping wall[J]．Tribology international．2012，46(1)：137-144．

[10] Xiaoqiang Yue，Hao Zhang，Congshu Luo，et al．Parallelization of a DEM code Based on CPU-GPD Heterogeneous Architecture[J]．Paraller Computational Fluid Dynamics，2014，405：149-159．

[11] 柳成文，毛靖儒，俞茂錚．90°彎管內(nèi)稀疏氣固兩相流及固粒對(duì)壁面磨損量的數(shù)值研究[J]．西安交通大學(xué)學(xué)報(bào)，1999，33(9)：53-57．

[12] 孫其誠(chéng)，王光謙．顆粒物質(zhì)力學(xué)導(dǎo)論[M]．北京：科學(xué)出版社，2009．

[責(zé)任編輯：崔海瑛]

(英文摘要略)

DOI 10.13356/j.cnki.jdnu.2095-0063.2016.06.020

Numerical simulation of particle flow in bend pipe based on discrete element method

LIU Hong

(School of Science in Anhui University of Science and Technology, Huainan 232001, China)

劉紅(1992-)，女，河南信陽(yáng)人，碩士研究生，主要從事Petri網(wǎng)理論與應(yīng)用及并行計(jì)算方面的研究。

TQ022.3

A

2095-0063(2016)06-0089-05

2016-06-28

DOI 10.13356/j.cnki.jdnu.2095-0063.2016.06.019