一种新型螺旋内槽管的气液固三相流数值模拟研究

， ，

(北京化工大学 化学工程学院，北京 100029)

1 引 言

螺旋内槽管广泛应用于化工换热器行业，特别是在换热器的强化传热过程中。与光滑管相比，螺旋内槽管会产生强烈的扰动，出现二次流的现象，使其内部的流体呈螺旋状流动。二次流现象的出现增强了壁面边界的扰动，从而使管内流体的对流换热系数得到明显的提高，进而提高了管道的传热性能。有学者从实验和数值模拟两个角度，针对不同螺旋内槽管的结构和形状进行了传热性能研究，包括不同的内螺纹形状[1-3]、螺纹的旋转角度[4]以及螺纹的头数[4]等。天然气水合物是水分子与天然气在特定温度和压力条件下形成的一种非化学计量的笼型结构冰状晶体物[5]。在水合物制备过程中，存在的主要问题是水合物不能连续、有效地生成及反应热不能快速移除。在本文提出的一种制备水合物的多管气-液鼓泡式反应结晶器[6]中，反应的管道为圆管，管内为气(天然气)液(水)固(反应生成的水合物)三相流动。为了提高换热效率，及时移走反应热，促进水合物的生成，本文提出利用螺旋内槽管来代替圆管。但是现有的关于螺旋内槽管的研究大多集中在传热性能方面，对于螺旋内槽管内复杂多相流现象的研究有限。

随着CFD(Computational Fluid Dynamics)技术的发展，气液固三相流动的数值模拟得到了迅速发展[7,8]。现有的CFD方法主要包括欧拉法和拉格朗日法，欧拉法主要对连续相流体在欧拉框架下求解N-S方程，对粒子相也在欧拉框架下求解颗粒守恒方程，以空间点为对象；欧拉-拉格朗日法对连续相流体在欧拉框架下求解N-S方程，对粒子相在拉格朗日框架下求解颗粒轨道方程，以单个粒子为对象。欧拉-欧拉-欧拉三流体模型主要用于致密颗粒流，因为其计算时间短，可以看作是连续介质的固相。在处理颗粒流的各种理论方法中，颗粒动力学理论(KTGF)得到了最广泛的应用，该理论主要基于Chapman等[9]所描述的稠密气体的经典动力学理论。

本文采用欧拉-欧拉-欧拉三流体模型结合颗粒动力学的理论对一种新型的水合物制备单元——螺旋内槽管管内的气(天然气)-液(水)-固(水合物)三相流动过程进行数值模拟。考察不同的表观速度(0.3 m/s，0.5 m/s和0.7 m/s)，水合物粒径(500 μm，750 μm和1000 μm)，气泡大小(10 μm，100 μm和1000 μm)，螺距(400 mm和800 mm)，螺纹头数(12和20)及螺纹旋向对于管内三相流动特性的影响。

2 数值模型的建立

2.1 数值模型

为了描述和预测管道的三相流动状态，建立描述预测三相流体力学的数学模型，将气体(天然气)、液体(水)和固体(天然气水合物)假设为三个独立的连续相，计算域空间各点三种流体都有各自不同的速度、温度和密度，这些流体存在于同一空间并相互渗透，但各有不同的体积分数，相互间有滑移。采用欧拉-欧拉-欧拉方法，结合颗粒动力学理论来模拟三相管道系统。

下面对控制方程进行一下简单的概述。当无相间质量传递时，气液固三相的连续性方程为

(1)

式中αq,ρq和uq分别代表三相的体积分数，密度和速度矢量。同时三相体积分数满足:

αl+αg+αs= 1

(2)

三相的动量方程分别为

(3)

(4)

(5)

(6)

式中μq为有效粘度。

固体压力和应力的描述采用颗粒动力学理论，具体计算如方程(7～11)所示。

固体压力ps[10]为

ps=αsρsΘs+2(1+e)α2sg0ρsΘs

(7)

式中e为颗粒-颗粒碰撞恢复系数，是一类经验参数。

体积粘度λs[10]为

(8)

径向分布函数g0[11]为

g0= [1-(αs/αs,max)1/3]-1

(9)

颗粒温度Θs[12]是颗粒脉动速度大小的量度，与颗粒脉动速度的平方成正比，如式(10)所示。

(10)

固体的剪切粘度/有效粘度μs[13]为

(11)

方程(3～5)的Mi,g,Mi,l和Mi,s是相间动量交换项，包括不同的相间作用力，如曳力、浮力、虚假质量力、BASSET力、MAGNUS力、SAFFMAN力和热泳力等[14]，在这些力中，曳力通常要比其他的力大得多[15]。因此在模拟中，相间作用力只考虑曳力的作用。本文液固和气固作用力模型采用Gidaspow模型[16]，气液之间采用schiller-naumann模型[17]。同时为了使上述流动方程组封闭，需要引入湍流模型。由于湍流的复杂性，很难用一种模型对所有的湍流问题都适用，目前标准k-ε湍流双方程模型在解决工程实际问题方面得到广泛应用，包括液-液旋流分离过程[18]、液固两相的旋风分离过程[19]、液滴在气体中的剪切破碎过程[20]、循环流化床(CFB)提升管内的气固两相流动过程[21]和气液固三相搅拌过程[22]等，因此本文采用标准k-ε模型对其进行计算。

2.2 模拟对象及数值方法

模拟对象为一种新型的螺旋内槽管，以20头螺旋内槽管为例，建立其管内流域三维几何模型如图1(a)所示，物理参数列入表1。网格划分全部采用结构化网格，如图1(b)所示。三相流模拟采用三维非稳态隐式求解，液相为不可压缩连续相，气相为单一气泡，固相采用颗粒流理论，单一粒径；压力-速度耦合采用 phase -coupled SIMPLE方法，动量和颗粒温度方程的对流项的离散格式采用二阶迎风格式，体积分数方程采用QUICK格式。

2.3 边界条件及操作

模拟压力为水合物生成实验压力8 MPa，液相为水相，气相为天然气，固相为水合物。进口条件为速度进口，气液固三相以特定速度流动，重力方向与流动方向相反；出口条件为自由出口；气相和液相采用无滑移壁面条件；固相采用Johnson等[23]提出的壁面条件处理。与光滑管相比，由于多头螺旋内槽的作用，使得水合物颗粒与管壁的碰撞频率增大，同时在螺旋内槽相互连接处也会增大能量损耗。目前尚未发现内螺旋管的相关研究，根据分析，可适当减小颗粒-壁面碰撞恢复系数对其进行修正，在本次工作中选用的颗粒-壁面碰撞恢复系数为0.9，镜面反射系数为0.05。假设模拟过程无化学反应，无热量传递。气液固三相的具体物性列入表2。

图1 20头螺旋内槽管三维几何模型及网格划分

Fig.1 Computational domain and mesh of the 20-head spiral-grooved tube

表2 三相物性

Tab.2 Physical properties of the three phases

物性名称物性值液相密度ρl/(kg·m-3)1000液相粘度μl/(kg·(m·s)-1)1.151×10-3颗粒密度ρs/(kg·m-3)900颗粒粒径ds/μm500,750,1000气体密度ρs/(kg·m-3)66气泡大小dg/μm10,100,1000三相表观流速up/(m·s-1)0.3,0.5,0.7进口固相体积分数αs0.35进口气相体积分数αg0.15颗粒-颗粒碰撞恢复系数0.9 固相最大体积分数αs,max0.63

3 结果与讨论

3.1 网格无关性分析

为了得到网格无关的解，以20头螺旋内槽管(螺距为400 mm)为例，在满足拟流体假设的前提下，划分了不同的网格数(48240，84420和253260)，同时对其流动情况进行数值模拟验证。为了便于比较，选取出口横截面沿X轴方向上的固含率作为评判指标，图2给出了不同操作条件下(包括不同的管内平均流速、颗粒粒径、气泡大小和三相体积分数)，不同网格数量的模拟结果。可以看出，三条线的值基本重合，说明当网格数量达到一定程度后，网格数量的增加对于结果基本没有影响。可适当选择较大的网格以节省计算时间，因此选用84420个网格进行计算。

3.2 气液固流速对于管内流动状态和气液固分布的影响

在螺距为400 mm，水合物颗粒粒径和气泡大小分别为750 μm和1000 μm条件下，探究气液固三相流速分别为0.3 m/s,0.5 m/s和0.7 m/s时，气液固三相在新型螺旋内槽管的流动状况和分布情况。

模拟得到流速up=0.7 m/s时，轴面处的水合物和天然气的流动状态及分布，结果如图3和图4所示，近进口和出口处的状态放大如图右侧所示。可以看出，在流体流入管内后，流体状态迅速达到稳定。壁面处的固体颗粒和天然气的体积分数都比管中心的少，同时水相集中在管壁。这主要是由于内槽管对流体的扰动使管内产生了复杂的二次流现象，在产生的离心力的作用下，由于水合物和天然气的密度比水的密度小，导致管壁处的水合物和天然气含量比管内的含量低得多。在up=0.7 m/s时，壁面处的水合物和天然气的最小体积分数分别为0.144和0.048。管径向存在一个明显的速度梯度，管中心的速度最大，同时壁面处由于壁面的作用速度很快减小为0。另外，由于螺旋内槽的作用，流速加快，管内的流速高于进口处的流速。与进口的0.7 m/s相比,管中心的水合物和天然气的最大流速分别为1.04 m/s和1.053 m/s。

在不同流速下，出口横截面处的流动状态及分布情况如图5～图7所示。可以看出，在不同的表观流速下，气液固三相流动状态存在一定的差别。当流速从0.3 m/s 增大到0.7 m/s 时，管壁处水合物和天然气的最小体积分数分别从0.261和0.103 减小到0.144和0.048，同时管壁处的水相的体积分数从0.628增大到0.807。从图5(d)，图6(d)和图7(d)可以看出，管中心的水合物的最大流速分别为0.461 m/s，0.754 m/s 和1.040 m/s。为了更好地进行比较，三个进口流速下，出口截面处沿着X轴方向的水合物和天然气的体积分数分布如图8 所示。可以看出，进口流速越大，壁面处的水合物和天然气的体积分数越小；由于天然气的密度小于水合物和水的密度，天然气更多集中在管中心，越靠近管壁含量越少，同时水合物的含量在壁面处最少，管中心其次，最多集中在管壁和管中心之间。

图2 出口截面处沿X轴方向的水合物的体积分数分布

Fig.2 Hydrate volume fraction along theXaxis direction in the outlet cross-section

图3 轴面处水合物的体积分数分布及流动状态(up=0.7 m/s)

Fig.3 Contours of axial surface hydrate flow state distribution (up=0.7 m/s)

图4 轴面处天然气的体积分数分布及流动状态(up=0.7 m/s)

Fig.4 Contours of axial surface natural gas flow state distribution (up=0.7 m/s)

图5 出口截面处的流动状态分布(up=0.3 m/s)

Fig.5 Contours of outlet cross-section flow state distribution (up=0.3 m/s)

图6 出口截面处的流动状态分布(up=0.5 m/s)

Fig.6 Contours of outlet cross-section flow state distribution(up=0.5 m/s)

图7 出口截面处的流动状态分布(up=0.7 m/s)

Fig.7 Contours of outlet cross-section flow state distribution(up=0.7 m/s)

3.3 水合物颗粒粒径对于管内流动状态和气液固分布的影响

在流速为0.3 m/s，螺距为400 mm，气泡大小为1000 μm，水合物颗粒粒径分别为500 μm和1000 μm的条件下，对螺旋内槽管的出口截面处流动状态进行模拟，结果如图9和图10所示。比较 图5(b)、图9(b)和图10(b)可以看出，颗粒粒径对壁面处的水合物含量有较大的影响，颗粒粒径越大，壁面处的水合物含量越小，管中心的含量越高，最小体积分数从0.327减s小到0.211。这主要是因为水合物和水之间的密度差比较小，粒径小的水合物不容易与水进行分离。但是，比较 图5(c)、图9(c)和图10(c)，颗粒粒径对于壁面处的天然气含量影响较小。比较图5(d)、图9(d)和图10(d)可得，水合物颗粒粒径的变化对于水合物的流速的影响不大，流速在0.45 m/s到 0.47 m/s 之间。

为了更好地比较，三种水合物颗粒粒径下，出口截面处沿着X轴方向的天然气和水合物的体积分数分布如图11所示。可以较为直观地对比得到，颗粒的粒径越大，壁面处的水合物含量越小，而对于天然气的分布基本上没有什么影响。

图8 出口截面处沿X轴方向的天然气和水合物的体积分数分布(ds=750 μm，dg=1000 μm)

Fig.8 Natural gas and hydrate volume fraction along theXaxis direction in the outlet cross-section (ds=750 μm，dg=1000 μm)

图9 出口截面处的流动状态分布 (ds=500 μm)

Fig.9 Contours of outlet cross-section flow state distribution (ds=500 μm)

图10 出口截面处的流动状态分布(ds=1000 μm)

Fig.10 Contours of outlet cross-section flow state distribution (ds=1000 μm)

3.4 气泡大小对于管内流动状态和气液固分布的影响

在流速为0.7 m/s，螺距为400mm，颗粒粒径为750 μm，气泡大小分别为10 μm和100 μm的条件下，对螺旋内槽管的出口截面处流动状态进行模拟，结果如图12和图13所示。比较图7(b)、图12(b)和图13(b)，可以看出，气泡大小对于壁面处的水合物含量有较大的影响，气泡越大，壁面处的水合物含量越小，最小体积分数从0.346减小到0.144，这可以更直观地从图14(b)看到。比较图7(c)、图12(c)和图13(c)可知，气泡大小对于壁面处的天然气含量有明显的影响，气泡越大，壁面处的天然气含量越小，管中心的含量越高，最小体积分数从0.148减小到0.048，同样可以直观地从图14(a)看到。 比较图7(d)、图12(d)和图13(d)可得，气泡大小的变化对于水合物的流速的影响不大，最大流速在1.102 m/s到1.040 m/s之间。

3.5 螺距大小对于管内流动状态和气液固分布的影响

在流速为0.7 m/s，气泡大小和水合物颗粒粒径分别为1000 μm和750 μm，螺距大小为800 mm的条件下，对气液固三相在新型螺旋内槽管的流动状况和分布情况进行了研究，结果如图15所示。比较图7(b)和图15(b),可以看出，螺距大小的变化对于水合物的分布有较大的影响，螺距增大一倍，壁面处水合物的含量明显增多，水合物的分布不均匀，最小体积分数增大到0.180，但在流动方向的背风面体积分数达到0.22以上。同样，比较图7(c)和图15(c)可知，螺距大小对于壁面处的天然气含量有明显的影响，螺距增大一倍，导致壁面处的天然气的含量分布不均匀，壁面处的天然气最小体积分数增大到0.076，但在流动方向的背风面体积分数却在0.09以上，这可以很直观地从图16(a)看到。同时，螺距增大一倍，水合物在管中心的最大流速从1.04 m/s减小到0.711 m/s。出现上述现象主要是因为随着螺距的增大，螺旋流强度减小，导致二次流现象的减弱，进而加大了气液固三相分离的难度。

图11 出口截面处沿X轴方向的天然气和水合物的体积分数分布(up=0.3 m/s)

Fig.11 Natural gas and hydrate volume fraction along theXaxis direction in the outlet cross-section(up=0.3 m/s)

图12 出口截面处的流动状态分布(dg=10 μm)

Fig.12 Contours of outlet cross-section flow state distribution (dg=10 μm)

图13 出口截面处的流动状态分布(dg=100 μm)

Fig.13 Contours of outlet cross-section flow state distribution (dg=100 μm)

3.6 螺纹头数对于管内流动状态和气液固分布的影响

在流速为0.7 m/s，螺距大小为400 mm，气泡大小和水合物颗粒粒径分别为1000 μm和750 μm的情况下，考察不同的螺纹头数对于流动状态和气液固分布的影响。头数为12时，出口截面处的流动状态分布如图17所示。对比图7(b)和图17(b)，可以看出，螺纹头数为12时，壁面处的水合物体积分数为0.122，远远高于头数为20时的0.048；同时，壁面处的水合物的分布极其不均匀，特别是在流动方向的背风面，水合物的体积分数在0.2以上，壁面处的三相分离效果差。比较图7(c)和图17(c)可知,壁面处的天然气分布具有同样的规律，在流动方向的背风面，天然气的体积分数在 0.09以上。螺纹头数为12时，管中心的最大流速为 0.744 m/s，低于头数为20时的1.044 m/s，这主要是因为头数为12时，管内产生的螺旋流不够强烈，三相分离效果不好。

图14 出口截面处沿着X轴方向的天然气和水合物的体积分数分布(up=0.7 m/s，ds=750 μm)

Fig.14 Natural gas and hydrate volume fraction along theXaxis direction in the outlet cross-section(up=0.7 m/s，ds=750 μm)

图15 出口截面处的流动状态分布(pitch=400 mm)

Fig.15 Contours of outlet cross-section flow state distribution(pitch=400 mm)

图16 出口截面处沿着X轴方向的天然气和水合物的体积分数分布(up=0.7 m/s)

Fig.16 Natural gas and hydrate volume fraction along theXaxis direction in the outlet cross-section(up=0.7 m/s)

为了更直观地比较，不同头数下，螺旋内槽管出口截面处沿着X轴方向的水合物和天然气的体积分数分布如图18所示。减小螺纹头数，壁面处的水合物和天然气的含量增大，但是对于管中心水合物的体积分数影响不大。

3.7 螺纹旋向对于管内流动状态和气液固分布的影响

在流速为0.7 m/s，螺距大小为400 mm，头数为12，气泡大小和水合物颗粒粒径分别为1000 μm和750 μm的情况下，改变螺纹旋向为进口逆时针旋转，考察管内流动状态和气液固分布变化，如 图19 所示。对比图17和图19，可以看出，改变螺纹的旋向，管内三相的分布状况恰好相反，而对于水合物和天然气的分离效果影响较小，同时流速分布基本没有发生变化。

图17 出口截面处的流动状态分布(12-head)

Fig.17 Contours of outlet cross-section flow state distribution(12-head)

图18 出口截面处沿着X轴方向的水合物和天然气的体积分数分布(up=0.7 m/s)

Fig.18 Hydrate and natural gas volume fraction along theXaxis direction in the outlet cross-section(up=0.7 m/s)

图19 出口截面处的流动状态分布(进口逆时针)

Fig.19 Contours of outlet cross-section flow state distribution(anticlockwise rotation)

4 结 论

本文建立了欧拉-欧拉-欧拉三流体模型结合颗粒动力学理论对一种新型螺旋内槽管进行了模拟研究，考察了流体表观流速、颗粒粒径、气泡大小、螺距、螺纹头数及螺纹旋向对于管内流体流动状态和分布的影响。得到的主要结论如下。

(1) 内槽管对流体的扰动使管内产生了复杂的二次流现象，由于气液固三相间的密度差，在产生的螺旋流和离心力的作用下，水合物和天然气在管中心位置聚集，同时管壁处的含量减小，同时由于螺旋内槽的作用，流速加快，管内的流速高于进口处的流速。

(2) 流体表观流速越大，壁面处的水合物和天然气的体积分数越小。由于天然气的密度小于水合物和水的密度，更多的天然气集中在管中心，越靠近管壁含量越少，同时水合物的含量壁面最少，管中心其次，最多集中在管壁和管中心之间。

(3) 颗粒的粒径越大，壁面处的水合物含量越小，而对于天然气的分布基本上没有什么影响。气泡大小对于管内天然气和水合物的含量分布有明显的影响，气泡越大，壁面处的水合物和天然气的含量越小，管中心的含量越高。

(4) 螺距大小对于管内天然气和水合物的含量分布有较大的影响，螺距越大，螺旋流强度越小，壁面处的水合物和天然气的含量越高，分布越不均匀。对比了不同螺纹头数对于管内流动状态和三相分布影响，螺纹头数越多，螺旋流强度越大，气液固三相分离效果越好，壁面处的水合物和天然气的含量越小；同时螺纹旋向的改变对于三相的分离效果影响较小。

参考文献(References):

[1] Liu X Y,Jensen M K.Geometry effects on turbulent flow and heat transfer in internally finned tubes[J].JournalofHeatTransfer,2001,123(6):1035-1044.

[2] Kim J H,Jansen K E,Jencen M K.Analysis of heat transfer characteristics in internally finned tubes[J].NumericalHeatTransfer,PartA:Applications,2004,46(1):1-21.

[3] 张 政,张建文.螺旋内槽管内的层流流动与传热的数值模拟[J].高校化学工程学报,1995,9(2):125-132.(ZHANG Zheng,ZHANG Jian-wen.Numerical simulation of laminar fluid flow and heat transfer in spiral internal finned tube[J].JournalofChemicalEngineeringofChineseUniversities,1995,9(2):125-132.(in Chinese))

[4] 陶文铨.传热与流动问题的多尺度数值模拟:方法与应用[M].北京：科学出版社,2009.(TAO Wen-quan.MultiscaleNumericalSimulationofHeatTransferandFlow:MethodsandApplications[M].Beijing:Science Press,2009.(in Chinese))

[5] Sloan E D.ClathrateHydratesofNaturalGas[M].Marcel Dekker,1998.

[6] 张建文,杨 林.一种制造天然气水合物的方法:CN 103571557A[P].2014-02-12.(ZHANG Jian-wen,YANG Lin.Method for producing natural gas hydrate:CN 103571557A[P].2014-02-12.(in Chinese))

[7] Li W L,Zhong W Q.CFD simulation of hydrodynamics of gas-liquid-solid three-phase bubble column [J].PowderTechnology,2015,286:766-788.

[8] Jia X Q,Wen J P,Zhou H L,et al.Local hydrodynamics modeling of a gas-liquid-solid three -phase bubble column [J].AIChEJournal,2007,53(9):2221-2231.

[9] Chapman S,Cowling T G.TheMathematicalTheoryofNon-uniformGases:AnAccountoftheKineticTheoryofViscosity,ThermalConduction,andDiffusioninGases[M].Cambridge University Press,2010.

[10] Lun C K K,Savage S B,Jeffrey D J,et al.Kinetic theories for granular flow:inelastic particles in Couette flow and slightly inelastic particles in a general flow-field[J].JournalofFluidMechanics,1984,140:223-256.

[11] Ding J,Gidaspow D.A bubbling fluidization model using kinetic theory of granular flow [J].AIChEJournal,1990,36(4):523-538.

[12] Hartge E U,Lars R,Wischnewski R,et al.CFD -simulation of a circulating fluidized bed riser [J].Particuology,2009,7(4):283-296.

[13] Zhou X Y,Gao J S,Xu C M,et al.Effect of wall boundary condition on CFD simulation of CFB risers[J].Particuology,2013,11(5):556-565.

[14] Olmos E,Gentric C,Vial C,et al.Numerical simulation of multiphase flow in bubble column reactors.Influence of bubble coalescence and break-up[J].ChemicalEngineeringScience,2001,56(21-22):6359-6365.

[15] Sokolichin A,Eigenberger G,Lapin A.Simulation of buoyancy driven bubbly flow:Established simplifications and open questions[J].AIChEJournal,2004,50(1):24-45.

[16] Gidaspow D.MultiphaseFlowandFluidization[M].Academic Press,1994.

[17] Tabib M V,Roy S A,Joshi J B.CFD simulation of bubble column—An analysis of interphase forces and turbulence models[J].ChemicalEngineeringJournal,2008,139(3):589-614.

[18] 陆耀军,周力行,沈 熊.液-液旋流分离管中强旋湍流的k-ε数值模拟[J].计算力学学报,2000,17(3):267-272.(LU Yao -jun,ZHOU Li-xing,Shen Xiong.Numerical simulation of strongly swirling turbulent flows in liquid-liquid hydrocyclones using k-ε model [J].ChineseJournalofComputationalMechanics,2000,17(3):267-272.(in Chinese))

[19] Ghasemi N,Sohrabi M,Khosravi M,et al.CFD simulation of solid-liquid flow in a two impinging streams cyclone reactor:Prediction of mean residence time and holdup of solid particles[J].ChemicalEngineeringandProcessingIntensification,2010,49(12):1277-1283.

[20] 楼建锋,洪 滔,朱建士.液滴在气体介质中剪切破碎的数值模拟研究[J].计算力学学报,2011,28(2):210-213.(LOU Jian-feng,HONG Tao,ZHU Jian-shi.Numerical study on shearing breakup of liquid droplet in gas medium [J].ChineseJournalofComputationalMechanics,2011,28(2):210-213.(in Chinese))

[21] Kong L,Zhang C,Zhu J.Evaluation of the effect of wall boundary conditions on numerical simulations of circulating fluidized beds[J].Particuology,2014,13:114-123.

[22] Murthy B N,Kasundra R B,Joshi J B.Hollow self-inducing impellers for gas-liquid-solid dispersion:Experimental and computational study[J].ChemicalEngineeringJournal,2008,141(1-3):332-345.

[23] Johnson P C,Jackson R.Frictional-collisional constitutive relations for granular materials,with application to plane shearing[J].JournalofFluidMecha-nics,1987,176:67-93.