孫 倩，彭天驥，嚴(yán) 安，周志偉,*

(1.清華大學(xué) 核能與新能源技術(shù)研究院，北京 100084；2.中國科學(xué)院 近代物理研究所，甘肅 蘭州 730000)

顆粒物質(zhì)存在于生產(chǎn)和生活的眾多方面。在食品、化工等工業(yè)領(lǐng)域，原材料通常以顆粒形式存儲(chǔ)于筒倉中。我國重點(diǎn)發(fā)展的球床模塊式高溫氣冷堆(簡稱高溫氣冷堆)，幾十萬球形燃料顆粒緊密堆積，在重力驅(qū)動(dòng)下流動(dòng)。中國科學(xué)院近代物理研究所創(chuàng)造性地提出了新型流態(tài)固體(mm級(jí)鎢)顆粒靶概念[1]。然而到目前為止，人們對顆粒流機(jī)理認(rèn)識(shí)不深，沒有完善的本構(gòu)模型能描述顆粒的流動(dòng)[2-3]。

為描述顆粒的流動(dòng)，根據(jù)研究目的和精度需求，提出了很多方法[4-6]，整體上可分為兩類：離散單元法(DEM)和連續(xù)介質(zhì)力學(xué)方法。離散單元法[6]是把整個(gè)介質(zhì)看作由一系列離散的獨(dú)立運(yùn)動(dòng)的單元所組成，單元本身具有一定的幾何和物理特征。其運(yùn)動(dòng)受經(jīng)典運(yùn)動(dòng)方程控制，整個(gè)介質(zhì)演化由各單元的運(yùn)動(dòng)和相互接觸來描述。連續(xù)介質(zhì)模型則忽略顆粒系統(tǒng)的離散性和單個(gè)顆粒的特性，利用連續(xù)性的假設(shè)，用統(tǒng)一的計(jì)算模型代替所有離散的顆粒。從物理角度，離散單元法更符合離散體系本身的特性。但從工程角度，將顆粒介質(zhì)進(jìn)行連續(xù)化處理，發(fā)展適合顆粒介質(zhì)大變形的數(shù)值方法，一方面能減少計(jì)算量，另一方面能復(fù)現(xiàn)顆粒集合的運(yùn)動(dòng)現(xiàn)象，更好地為災(zāi)害防治和工業(yè)應(yīng)用提供參考[7]。

本文基于物質(zhì)點(diǎn)方法，采用μ(I)流變模型，實(shí)現(xiàn)用連續(xù)性方法模擬密集顆粒流動(dòng)的計(jì)算框架，用于顆粒堆積坍塌和重力驅(qū)動(dòng)下顆粒在漏斗流動(dòng)過程中的模擬，并與物理實(shí)驗(yàn)和DEM計(jì)算結(jié)果進(jìn)行對比。

1 連續(xù)體本構(gòu)方程

根據(jù)流動(dòng)特點(diǎn)，顆粒流通常可分為準(zhǔn)靜態(tài)流、快速流和慢速流3種流態(tài)，分別對應(yīng)類固、類氣、類液狀態(tài)。當(dāng)顆粒處于準(zhǔn)靜態(tài)流狀態(tài)時(shí)，可用彈塑性理論對其進(jìn)行描述；對于處于快速流動(dòng)狀態(tài)的顆粒物質(zhì)，人們發(fā)展了顆粒動(dòng)理論來描述該流態(tài)下顆粒物質(zhì)的本構(gòu)關(guān)系[8]。而對于介于快速流與準(zhǔn)靜態(tài)流之間的過渡狀態(tài)，顆粒像流體一樣流動(dòng)，相互間又保持著相對持續(xù)的接觸，塑性理論和顆粒動(dòng)理論都不適用。Midi[9]和Jop等[10]受黏塑性和賓漢塑性流體行為的啟發(fā)，根據(jù)大量的數(shù)值和物理實(shí)驗(yàn)，提出了描述密集顆粒流動(dòng)的本構(gòu)方程，即局部流變模型，較好地復(fù)現(xiàn)了顆粒在不同邊界下的流動(dòng)。該模型在宏觀尺度上對顆粒介質(zhì)的耗散描述是通過有效摩擦系數(shù)μ=τ/p(τ為剪切應(yīng)力，p為正應(yīng)力)實(shí)現(xiàn)的。

τ=μ(I)p

(1)

式中，μ(I)為顆粒物質(zhì)體系的總摩擦系數(shù)，是復(fù)雜力鏈響應(yīng)的宏觀體現(xiàn)，與慣性系數(shù)I有關(guān)。

(2)

μ(I)流變模型認(rèn)為式(1)中的摩擦系數(shù)μ(μ=τ/p，圖1)可表示為：

(3)

式中：μs為顆粒介質(zhì)在零剪切率下的最小摩擦系數(shù)；μ2為高速流態(tài)下摩擦系數(shù)的上限；I0為常數(shù)。

圖1 摩擦系數(shù)μ(I)與慣性系數(shù)I的關(guān)系[10]Fig.1 Friction coefficient μ(I) as a function of inertial number I[10]

結(jié)合針對顆粒固態(tài)的彈塑性理論和適用于顆粒液態(tài)的μ(I)黏塑性模型，即可描述顆粒的固態(tài)行為和液態(tài)行為[11]。定義顆粒介質(zhì)的柯西應(yīng)力張量為σ，應(yīng)力偏張量σ0=σ-(trσ)I/3，顆粒介質(zhì)的運(yùn)動(dòng)速度為v。則速度梯度張量L=gradv，L可分解為旋率張量W=(L-LT)/2和變形率張量D=(L+LT)/2。D可認(rèn)為由其彈性部分De和塑性部分Dp組成，即D=De+Dp。

顆粒介質(zhì)的質(zhì)量守恒方程可表示為：

(4)

動(dòng)量守恒方程可表示為：

(5)

式中：b為顆粒介質(zhì)所受到的體積力；ρ為顆粒介質(zhì)的密度。

焦曼應(yīng)力率為：

(6)

當(dāng)顆粒處于彈性穩(wěn)定態(tài)時(shí)，體系內(nèi)無塑性應(yīng)變。本構(gòu)關(guān)系用彈性變形率表示為：

(7)

其中：E為楊氏模量；υ為泊松比。定義顆粒介質(zhì)的剪切力τ、壓力p和摩擦系數(shù)μ如下：

τ=‖σ0‖，p=-(trσ)/3，μ=τ/p

當(dāng)顆粒介質(zhì)的摩擦系數(shù)μ大于其屈服摩擦系數(shù)μs時(shí)，發(fā)生彈性失穩(wěn)，引發(fā)塑性流動(dòng)，這時(shí)μ(I)流變模型被用作黏塑性流動(dòng)法則。即：

(8)

2 物質(zhì)點(diǎn)計(jì)算方法

Sulsky等[12]將用于流體動(dòng)力學(xué)的質(zhì)點(diǎn)網(wǎng)格法(PIC)擴(kuò)展到固體力學(xué)問題中提出了物質(zhì)點(diǎn)法(MPM)。MPM采用攜帶材料所有信息的物質(zhì)點(diǎn)離散材料區(qū)域，以表征材料區(qū)域的運(yùn)動(dòng)和變形狀態(tài)，并避免了處理對流項(xiàng)；采用規(guī)則的歐拉背景網(wǎng)格計(jì)算空間導(dǎo)數(shù)和動(dòng)量方程，從而實(shí)現(xiàn)了質(zhì)點(diǎn)間的相互作用與聯(lián)系，并避免了網(wǎng)格畸變問題，非常適合處理涉及材料特大變形的問題[13-14]。

以顆粒介質(zhì)為例，物質(zhì)點(diǎn)法將材料區(qū)域離散為1組相對背景網(wǎng)格運(yùn)動(dòng)的質(zhì)點(diǎn)，如圖2所示。每個(gè)質(zhì)點(diǎn)表示1個(gè)材料團(tuán)簇，并攜帶其所有物質(zhì)信息，如質(zhì)量、速度、應(yīng)力和應(yīng)變等。背景網(wǎng)格僅用于求解動(dòng)量方程和物理量的空間導(dǎo)數(shù)，一般取為固定于空間的規(guī)則網(wǎng)格。

圖2 物質(zhì)點(diǎn)方法[13]Fig.2 Material point method[13]

將連續(xù)體離散為質(zhì)點(diǎn)后，其密度近似為：

(9)

物質(zhì)點(diǎn)法將連續(xù)體離散為一系列質(zhì)點(diǎn)，質(zhì)點(diǎn)攜帶的質(zhì)量在運(yùn)動(dòng)過程中保持不變，故質(zhì)量守恒方程(4)自動(dòng)滿足，因此主要求解動(dòng)量方程。動(dòng)量方程(5)的求解是基于其更新拉格朗日式的控制方程的弱形式[13]。

(10)

本文基于物質(zhì)點(diǎn)方法，根據(jù)μ(I)理論模型添加本構(gòu)關(guān)系式(7)、(8)，對密集顆粒的流動(dòng)進(jìn)行模擬。本文的整體計(jì)算流程如圖3所示。

3 顆粒介質(zhì)坍塌模擬

顆粒坍塌是顆粒流動(dòng)的一個(gè)典型現(xiàn)象[15]。清華大學(xué)孫其誠曾針對顆粒坍塌過程做過物理實(shí)驗(yàn)[16]。實(shí)驗(yàn)裝置如圖4所示，實(shí)驗(yàn)區(qū)域是1個(gè)60 cm×20 cm×20 cm的箱子，顆粒的初始堆積體積為10 cm×18 cm×20 cm。實(shí)驗(yàn)中使用了約24 000顆粒徑5 mm的陶土顆粒，顆粒的楊氏模量E為5.5 GPa，泊松比為0.3，顆粒介質(zhì)的表觀密度為2 200 kg/m3，陶粒堆積形成錐角的平均值為22°。

圖3 MPM計(jì)算流程Fig.3 Calculation procedure of MPM

采用MPM模擬圖4的實(shí)驗(yàn)，材料參數(shù)取μs=tan(22°)=0.404，μ2=0.64，I0=0.28[9-11]。模擬中共采用11 520個(gè)物質(zhì)點(diǎn)。在x=0和x=60 mm處設(shè)置對稱邊界，約束x方向動(dòng)量；底部為固定邊界，側(cè)壁采用對稱邊界。顆粒初始堆積的離散間距為2.5 mm，網(wǎng)格邊長為5 mm。DEM模擬采用開源軟件LIGGGHTS，選用Hertz-Mindlin接觸模型。參照文獻(xiàn)[15]的做法，在底部鋪設(shè)1層固定的微小顆粒來模擬實(shí)驗(yàn)中的底部粗糙砂紙。

圖4 顆粒坍塌實(shí)驗(yàn)裝置示意圖[3,16]Fig.4 Experimental setup of granular collapse[3,16]

圖5示出了實(shí)驗(yàn)、MPM模擬和DEM模擬坍塌過程中的顆粒介質(zhì)速率分布情況?？煽吹?，從0 s撤掉擋板后，顆粒體在重力作用下加速崩塌；到0.16 s時(shí)，顆粒速度達(dá)到最大值；到0.32 s時(shí)，重力的加速作用逐漸減弱，由于壁面和內(nèi)部的摩擦力，顆粒開始逐漸減速；到t=0.72 s時(shí)，大部分顆粒均已重新恢復(fù)靜止?fàn)顟B(tài)。

圖5 顆粒坍塌實(shí)驗(yàn)和MPM及DEM數(shù)值模擬結(jié)果對比Fig.5 Results of experiment, MPM and DEM for granular collapse

從圖5可明顯看出，MPM模擬速率分布和實(shí)測速率分布以及DEM模擬結(jié)果基本相同。

4 漏斗流動(dòng)

4.1 卸料速度

漏斗的一個(gè)重要參數(shù)是卸料速度。該參數(shù)影響顆粒材料在工業(yè)中的運(yùn)行狀態(tài)。因此，能否預(yù)測漏斗的卸料速度是判斷一個(gè)模型能否具有預(yù)測漏斗流動(dòng)的重要標(biāo)準(zhǔn)。Beverloo根據(jù)大量的實(shí)驗(yàn)數(shù)據(jù)，得到了矩形漏斗的卸料速度經(jīng)驗(yàn)關(guān)系式[8]：

(11)

本文模擬的矩形漏斗尺寸如圖6所示，H=1 m，L=1 m?？紤]到縱向深度W?D，故縱向上采用周期性邊界條件。漏斗沿中軸面對稱(圖中虛線所示)，模擬時(shí)只模擬一半幾何，中軸面處采用對稱邊界，側(cè)面壁面采用無摩擦邊界，底部壁面采用完全摩擦邊界。模擬時(shí)網(wǎng)格尺寸為0.01 m，材料參數(shù)與坍塌算例的一致。

圖6 矩形漏斗示意圖Fig.6 Schematic diagram of silo

圖7為D=0.14 m、H=1 m時(shí)漏斗卸料的瞬態(tài)速度分布。在漏斗出口正上方，顆?？焖倭鲃?dòng)；靠近壁面區(qū)域，顆粒稠密堆積，幾乎不流動(dòng)；漏斗中心軸附近的顆粒緊密堆積，緩慢流動(dòng)。

圖8為不同填充高度和出口尺寸下漏斗的剩余質(zhì)量變化。從圖8a可看出，當(dāng)漏斗出口寬度D=0.14 m時(shí)，不同高度下，漏斗的剩余質(zhì)量均隨時(shí)間線性減少，說明卸料過程是以穩(wěn)定的質(zhì)量流量進(jìn)行的。不同初始填充高度下，直線的斜率相同，說明初始填充高度對卸料速率沒有影響，這與經(jīng)驗(yàn)結(jié)果一致。從圖8b可看出，不同出口尺寸下，漏斗均以穩(wěn)定的質(zhì)量流量卸料，出口尺寸越大，斜率越大。

圖7 t=0、1、3、5、15 s時(shí)漏斗卸料的速度分布Fig.7 Velocity throughout discharge process at t=0, 1, 3, 5 and 15 s

圖8 不同填充高度和出口尺寸下漏斗剩余質(zhì)量變化Fig.8 Mass remaining in silo for different initial mass and different hole sizes

圖9 卸料率與出口尺寸D的關(guān)系Fig.9 Variation of discharge rate versus outlet size

4.2 ADS顆粒流靶

在加速器驅(qū)動(dòng)次臨界系統(tǒng)(ADS)中，散裂靶是耦合加速器與次臨界堆的核心部件。中國加速器驅(qū)動(dòng)嬗變研究裝置(CiADS)中的散裂靶采用中國科學(xué)院ADS研究團(tuán)隊(duì)研發(fā)的固體顆粒流方案[17]。該方案以流動(dòng)的固體顆粒為靶材，同時(shí)充當(dāng)自身的冷卻介質(zhì)，與質(zhì)子束流散裂反應(yīng)產(chǎn)生中子的同時(shí)，利用其流動(dòng)特性將高功率密度的束流沉積熱帶出束靶耦合的反應(yīng)區(qū)，兼具了固態(tài)靶和液態(tài)靶的優(yōu)勢。顆粒流靶作為一種新概念散裂靶，尚處于設(shè)計(jì)階段[18]，將顆粒流靶的模型簡化成1個(gè)錐形漏斗，幾何參數(shù)如圖10所示。

分別采用DEM和MPM對靶區(qū)的顆粒流動(dòng)進(jìn)行模擬。顆粒材料模擬參數(shù)列于表1。t=0 s時(shí)漏斗內(nèi)填滿顆粒材料，流動(dòng)開始后，顆粒材料會(huì)從漏斗入口不斷注入。顆粒流經(jīng)靶區(qū)后流出，在靶區(qū)會(huì)形成一個(gè)穩(wěn)定流動(dòng)的流場分布。兩種方法模擬得到的穩(wěn)定流動(dòng)z方向速度分布如圖11所示。圖12為MPM計(jì)算得到的顆粒質(zhì)量流量隨時(shí)間變化的結(jié)果?？煽闯觯?dāng)t=1 s時(shí)，流動(dòng)趨于穩(wěn)定，質(zhì)量流量為285.43 kg/s，DEM的穩(wěn)定質(zhì)量流量為288.5 kg/s。

圖10 顆粒流靶簡化模型Fig.10 Geometrical model of granular flow target

表1 顆粒流靶材料參數(shù)Table 1 Parameter for granular target

從圖11a可看出，由于壁面的摩擦作用，靠近壁面的一層顆粒會(huì)有“貼壁”現(xiàn)象，在壁面處形成較大的速度梯度。連續(xù)性方法模擬出的結(jié)果在主流區(qū)能與DEM結(jié)果吻合較好，但在靠近壁面的邊界處沒有模擬出顆粒在壁面處的“貼壁”現(xiàn)象。分析原因有如下兩點(diǎn)。

1) 連續(xù)性方法中，顆粒與壁面相互作用的理論模型不完善。顆粒與壁面之間的相互作用較復(fù)雜，目前大部分學(xué)者關(guān)注的仍為顆粒與顆粒之間的相互作用。對于顆粒與壁面的相互作用，目前的處理方式是給定屈服摩擦系數(shù)μw(μw為常數(shù))，當(dāng)壁面剪切力與壓力的比值超過此值時(shí)，顆粒開始流動(dòng)。這種處理方式過于簡化，且DEM中顆粒與壁面的摩擦系數(shù)，與將顆粒視為連續(xù)體宏觀上體現(xiàn)出來的與壁面之間的摩擦系數(shù)含義不同。關(guān)于DEM中顆粒-顆粒摩擦系數(shù)、顆粒-壁面摩擦系數(shù)與連續(xù)性方法中μ(I)流變模型參數(shù)、顆粒介質(zhì)-壁面摩擦系數(shù)μw兩套參數(shù)之間如何對應(yīng)還需進(jìn)一步研究。

圖11 顆粒流靶流動(dòng)結(jié)果對比Fig.11 Velocity of granular flow target

圖12 MPM計(jì)算的質(zhì)量流量隨時(shí)間的變化Fig.12 Variation of mass flow rate calculated by MPM with time

2) MPM中接觸算法不成熟。目前程序中處理顆粒與壁面的相互作用是將壁面視為剛體材料處理，考慮顆粒材料與剛體的相對滑動(dòng)，需在物質(zhì)點(diǎn)方法中引入接觸算法。接觸算法中，邊界處約束過強(qiáng)會(huì)導(dǎo)致數(shù)值不穩(wěn)定。

5 結(jié)論

1) 本文基于物質(zhì)點(diǎn)方法，結(jié)合針對顆粒固態(tài)的彈塑性理論和適用于顆粒液態(tài)的μ(I)流變模型，實(shí)現(xiàn)了用連續(xù)性方法模擬密集顆粒流動(dòng)的框架。

2) 在此框架下，模擬了顆粒介質(zhì)坍塌流動(dòng)，計(jì)算結(jié)果與物理實(shí)驗(yàn)和DEM模擬結(jié)果均吻合較好。

3) 對顆粒在重力驅(qū)動(dòng)下漏斗中的流動(dòng)進(jìn)行了模擬，卸料率結(jié)果與Beverloo經(jīng)驗(yàn)公式吻合較好。對于ADS顆粒流靶流動(dòng)，主流區(qū)速度分布與DEM結(jié)果吻合較好，但無法模擬出壁面處顆粒的“貼壁”現(xiàn)象，經(jīng)過分析認(rèn)為這是顆粒與壁面相互作用的理論模型不完善以及MPM中接觸算法不成熟導(dǎo)致的，這些問題需要在后續(xù)研究中加以解決。

4) 可在本文探究的連續(xù)性方法模擬密集顆粒流動(dòng)框架的基礎(chǔ)上，完善邊界條件和本構(gòu)模型，為用連續(xù)性方法處理大量顆粒的工程問題提供借鑒。

本文所用程序基于清華大學(xué)張雄教授所提供的開源三維物質(zhì)點(diǎn)法程序(MPM3D)改編而成，文中的顆粒坍塌實(shí)驗(yàn)數(shù)據(jù)由清華大學(xué)水利系孫其誠老師提供，在此對兩位老師一并表示感謝。