基于介观结构的饱和与非饱和多孔介质有效应力李锡夔;杜友耀;段庆林【摘要】基于描述含液颗粒材料介观结构的Voronoi胞元模型和离散颗粒集合体与多孔连续体间的介-宏观均匀化过程,定义饱和与非饱和多孔介质有效应力.导出了计及孔隙液压引起之颗粒体积变形的饱和多孔介质广义有效应力.用以定义广义有效应力的Biot系数不仅依赖于颗粒材料的多孔连续体固体骨架及单个固体颗粒的体积模量(材料参数),同时与固体骨架当前平均广义有效应力及单个固体颗粒的体积应变(状态量)有关.提出了描述非饱和多孔介质中非混和固体颗粒、孔隙液体和气体等三相相互作用的具介观结构的Voronoi胞元模型.具体考虑在低饱和度下双联(binary bond)模式的摆动(pendular)液桥系统介观结构.导出了基于介观水力-力学模型的非饱和多孔介质的各向异性有效应力张量与有效压力张量.考虑非饱和多孔介质Voronoi胞元模型介观结构的各向同性情况,得到了与非饱和多孔连续体理论中唯象地假定的标量有效压力相同的有效压力形式.但本文定义的与确定非饱和多孔介质有效应力和有效压力相关联的Bishop参数由基于三相介观水力-力学模型,作为饱和度、孔隙度和介观结构参数的函数导出,而非唯象假定.【期刊名称】《力学学报》【年(卷),期】2016(048)001【总页数】12页(P28-39)【关键词】饱和与非饱和多孔连续体;含液离散颗粒集合体;介观结构的Voronoi胞元模型;有效应力;Biot系数;Bishop参数【作者】李锡夔;杜友耀;段庆林【作者单位】大连理工大学工业装备结构分析国家重点实验室,大连116024;大连理工大学工业装备结构分析国家重点实验室,大连116024;大连理工大学工业装备结构分析国家重点实验室,大连116024【正文语种】中文【中图分类】O354.4饱和与非饱和颗粒材料由众多颗粒与其间充满或部分填充液体的孔隙组成.它在宏观尺度通常模型化为饱和或非饱和多孔连续体.Terzaghi[1]提出了饱和多孔介质力学的有效应力定义.基于有效应力原理，Biot[2-5]建立了多孔介质中静力和动力响应流固相耦合作用的控制方程.考虑孔隙液压所引起之颗粒体积变形对饱和多孔介质固体骨架变形的效应，Zienkiewicz和Shiomi[6]引入了Biot系数，提出了饱和多孔介质广义有效应力定义.Kytomaa等[7]将饱和颗粒材料在介观尺度模型化为含液离散颗粒集合体，研究了孔隙液压作用引起之单个颗粒压缩变形对饱和多孔介质本构行为的影响.Bishop[8-9]，Skempton[10]，Li等[11]，Zienkiewicz等[12]将有效应力概念推广到含两相或多相不可混孔隙流体的非饱和多孔介质中.Bishop[8-9]首先把有效应力概念推广到非饱和多孔介质，用平均孔隙压力表示孔隙多相流体效应.把有效应力概念简单地推广到非饱和多孔介质的论据不如它在仅有单一孔隙流体存在的饱和多孔介质中那么清楚，对非饱和多孔介质中各种有效应力定义仍有争论[10,12-13].Gray和Schrefle[14]在热动力学框架内对非饱和多孔介质有效应力概念作出了一些澄清.Gray和Miller[15]把热动力约束平均理论引入到多孔介质系统，所提供的一致性框架有助于对Gray和Schrefle[16]所提出的非饱和多孔介质中的有效应力和总应力张量中的物理量作出解释.在国内，陈正汉等[17]以热力学和连续介质力学的应力理论为基础，导出了一系列非饱和土应力状态变量;赵成刚等[18]基于多相孔隙介质理论提出了非饱和土广义有效应力原理;邵龙潭等[13]阐明了有效应力、粒间应力和土骨架应力的物理意义及其相互关系；刘艳、赵成刚等[19]综述了非饱和土力学理论的研究进展.然而，对于饱和、非饱和多孔介质，还未见基于介观结构特性和介观水力--力学响应状态的有效应力研究的工作发表.利用介--宏观均匀化过程，人们基于介观水力--力学模型致力于发展定义非饱和多孔介质有效应力与有效压力的方法[20-23].有效压力不再由非饱和多孔连续体理论[9,24]中所定义的标量表示，而是一个依赖于由固体颗粒、孔隙液体和气体组成的介观结构、反映毛细力各向异性效应的张量.在介观尺度上，颗粒材料的每个离散颗粒在运动学上除平移自由度还具独立的旋转自由度，在动力学上由一个颗粒到其相邻接触颗粒能传递力偶.为表征离散颗粒介观结构，颗粒材料应模型化为在每个材料点定义具有独立旋转自由度的等效多孔Cosserat连续体[25-29].考虑饱和或非饱和颗粒材料，它在介观尺度模型化为含液离散颗粒集合体；而在宏观尺度上模型化为在一个材料点处同时存在固相和饱和液相、或固相与不可混液、气相的饱和或非饱和等效多孔Cosserat连续体元.尽管如此，孔隙流压对多孔连续体的效应可假定仅对总Cauchy应力的有效Cauchy应力和有效压力的定量剖分有关，而与定义于Cosserat连续体的偶应力无关.定义于宏观连续体一个材料点的有效压力张量依赖于该材料点由固、液、气三相组成的介观结构.本文中指定的饱和或非饱和颗粒材料的介观结构由一个参考颗粒与其直接相邻颗粒以及其间的间隙流体组成的两相或三相Voronoi胞元模型描述.它基于表示干颗粒材料介观结构的单相Voronoi胞元模型[30].湿Voronoi胞元模型不仅包含位于湿Voronoi胞元内的参考颗粒与间隙流体，且包含湿Voronoi胞元外的直接相邻颗粒.它能用以描述包含局部孔隙比、局部液相和气相饱和度、参考颗粒与其直接相邻颗粒接触拓扑的局部介观结构.由湿Voronoi胞元模型导出的有效压力张量依赖于胞元内液桥系统的几何特征.液桥结构的构形随孔隙液体的饱和度增长可从摆动状态演化到索状和毛细状态.本文将具体考虑离散颗粒系统在低饱和度下以具双联（binary bond）模式的摆动（pendular）构形液桥为特征的介观水力--力学模型[31].基于湿颗粒材料的介--宏观均匀化过程与平均场理论定义基于介观水力--力学模型的饱和多孔介质广义有效应力和Biot系数，非饱和多孔介质有效应力、有效压力和Bishop参数.1.1 饱和多孔介质有效应力图1表示饱和多孔介质的二维饱和Voronoi胞元模型.本节用它概念性地描述由一簇浸透孔隙液体的球体组成的三维饱和Voronoi胞元模型.定义Voronoi胞元边界内颗粒为参考颗粒.以m表示半径为r，体积为Vp的参考颗粒的直接接触颗粒数，V表示定义于参考颗粒处的Voronoi胞元体积.以φ表示Voronoi胞元的孔隙度，Voronoi胞元内的孔隙体积则可表示为Vv＝V−Vp＝φV.而Vp可表示为Vp＝（1−φ）V.基于由饱和离散颗粒集合体到饱和多孔连续体的介--宏观均匀化，定义于Voronoi 胞元的饱和多孔连续体的平均Cauchy应力可表示为式中，Vl是Voronoi胞元内孔隙液体体积，对饱和多孔介质有Vl＝Vv;σji是施加于体积为Vp的固体颗粒内任一物质点的Cauchy应力；pl是体积为Vl的孔隙液体中各点孔隙液压，定义压为正.假定Vl中pl为常数，式（1）可简化为式（2）第2个等号右端第1项，根据Gauss散度定理和平衡方程可展开为式中Γp是Voronoi胞元内参考颗粒表面，m是与它直接接触相邻颗粒数，ti是作用于颗粒表面xj处的表面力，nj为该处颗粒表面单位外法线方向.表面力ti包含作用于该处的颗粒间接触力（c＝1,2,···,m）和分布液压pl.式（3）描述了作用于参考颗粒的颗粒间接触力和颗粒表面孔隙液压对于Voronoi胞元总平均Cauchy 应力的贡献.将式（3）代入式（2）得到式中饱和多孔介质有效Cauchy应力'ji定义为式（4）可改写为用总Cauchy应力ji和孔隙液压pl表示的有效Cauchy应力，即它再现了Terzaghi[1]提出的有效应力定义.式（5）表示，在宏观多孔连续体中定义为控制一个材料点处固体骨架力学行为的有效应力表征了该处局部离散颗粒系统的介观结构和力链强度.式（6）表明饱和多孔介质中孔隙液压效应的各向同性. 1.2 饱和多孔介质中广义有效应力利用式（5），略去接触颗粒间的耗散摩擦，由模型化为多孔Cosserat连续体元的Voronoi胞元导出基于介观力学的本构关系为[30]式中是由Voronoi胞元模型定义的等效多孔Cosserat连续体元的应变，线位移空间导数，微转角，微曲率.式（7）中弹性模量张量为[30]式中是参考颗粒与第c个直接相邻颗粒的当前重迭量，是依赖于的Heaviside单位函数，体现了颗粒接触的丧失和产生及其随时间的演变（c＝1,2,···,m）表示第c 个直接相邻颗粒 B的半径.分别是两接触颗粒间相应于切向滑动、滚动摩擦力与法向接触力的刚度系数.代表参考颗粒在接触点c处Voronoi胞元边界的单位切向量和法向量.计及单个颗粒在静水压力pl下的可压缩性，其各向同性体积应变εv可表示为式中Ks,Es,νs是单个颗粒的体积模量、杨氏模量和泊松比；b＝（1−2νs）/Es.静水压力pl作用下初始半径为r0的颗粒的当前半径r可给出为[7]式中定义其中并定义εm，εv压为正.考虑由静水压力引起之颗粒体积应变对基于介观力学的本构关系的影响，由式（7）可给出上式右端第2项考虑了作用于等效多孔Cosserat连续体的有效应力随作用于固体颗粒的静水压力的变化;体现了作用于固体颗粒的静水压力对饱和多孔连续体的模量张量及应变的影响.其中忽略了参考颗粒及其直接相邻颗粒间孔隙液体压力梯度的影响.由此可得其中上式中的计算需要考虑随参考颗粒半径r、直接相邻颗粒半径及基于Voronoi胞元模型的等效多孔Cosserat连续体的体积变化.将式（15）和式（16）代回式（14），且略去下标r0可得定义αji和βji如下并定义与模型化为等效Cosserat连续体的固体骨架应变度量相关联的广义有效应力把式（19）和式（20）代入式（18）得到由式（6）所定义的有效应力可得广义有效应力与总Cauchy应力间存在如下联系由式（9）可得，对式（20）两边同时乘δji可得由此可以假定βji＝βδji，同时由上式可得将上式代回式（22）和式（23）可得由式（19）和式（26）可得基于介观水力--力学的Biot系数，αji为二阶张量，其对广义有效应力的影响为各向异性，且一般地对其剪切分量也有影响.Biot系数αji不仅与介观结构的材料参数，如颗粒半径、固体颗粒的体积模量、颗粒间接触刚度有关，同时与介观结构的模量张量固体颗粒的体积应变、等效连续体的应力状态有关.为与由饱和多孔连续体宏观理论导出的标量 Biot系数比较，需考虑αji的标量形式.对式（19）两边同时乘δji可得由式（9）可得，同时由式（21）可得上式中为对应于多孔连续体体积应变的广义有效应力的平均应力.定义由式（27）～式（29），可以假定将式（28）～式（30）代回式（27）可以得到基于介观结构与响应的标量形式Biot系数一般地，饱和多孔介质固体骨架承受压缩体积应变，即为与宏观饱和多孔连续体理论中由广义Biot公式给出的Biot系数比较，式（31）可改写为一般地由式（8）及式（29）可得由式（32）和式（33）可得多孔连续体广义体积模量为将式（30）代回式（26b）可得如下基于介观水力--力学的总应力与广义有效应力的关系为值得指出，式（32）与在宏观饱和多孔连续体途径中获得的广义有效应力[6,12]具有相同的形式.然而，由文献[6，12]定义的Biot系数α是一个常系数而在式（32）中所定义的Biot系数α是一个不仅依赖于材料参数KT,KS，同时也依赖当前状态变量和εm（即当前时刻的等效连续体的应力状态和颗粒因孔隙液压而产生的体积应变）的变系数.当单个颗粒的平均体积模量KS远大于即因孔隙液体压引起的颗粒体积应变对广义有效应力的影响可略去.由式（6），式（35），广义有效应力退化为经典有效应力即Bishop[8]，Skempton[10]，Bishop 和Blight[9]，Coussy[24]，Zienkiewicz等[12]，Gray和Schrefle[16]在多孔连续体理论框架内把由式（6）所表示的饱和多孔介质有效应力概念推广到非饱和情况.考虑到非饱和多孔介质中孔隙气体（或其他干相孔隙流体（如孔隙油等））的存在，孔隙压力由平均孔隙流体压力定义[8,10-12]即有效应力定义修改为式中依赖于饱和度Sl的加权参数χ称为Bishop参数，Pji和 Qji分别是非饱和多孔连续体理论中定义为对角线二阶张量的有效压力和基质毛细压力张量，即细观力学研究表明[20]，由于液桥分布、液桥力和分支向量分布通常呈各向异性，由非饱和孔隙液体引起的毛细力将导致非饱和多孔介质的诱导各向异性.非饱和多孔介质总Cauchy应力σji应表示为式中控制固体骨架变形的有效应力表征了颗粒间接触力；Pji是表示非饱和多孔介质中各向异性毛细效应的二阶张量，它依赖于局部处液桥介观结构、孔隙气体（干相流体）和孔隙液体（湿相流体）压力差和液桥处表面张力.二维非饱和Voronoi胞元模型如图2所示.它由Voronoi胞元内的参考颗粒及其位于Voronoi胞元几何边界外的直接相邻颗粒、以及Voronoi胞元内的孔隙流体组成.在低饱和度下假定颗粒间孔隙液体以“摆动”（pendular）液桥形态存在，如图2所示.二维非饱和 Voronoi胞元模型描述了模型化为非饱和离散颗粒系统的非饱和颗粒材料介观结构，同时也表示了在参考颗粒处定义的具有介观结构的等效非饱和多孔Cosserat连续体元. 以表示Voronoi胞元内第i个液桥处液体体积，表示Voronoi胞元内围绕参考颗粒两个相邻液桥间的孔隙气体体积.可表示式中以和分别表示和中平均孔隙液压和气压，则Voronoi胞元内Vl和中平均孔隙液压和气压可定义为Voronoi胞元的总Cauchy应力张量在均匀化意义下可定义为利用式（45），式（46）右端第二、第三项可分别表示为式中已定义应用Gauss定理和平衡条件，略去重力效应，式（46）右端第1项积分可表示为式中作用于颗粒表面 xj处的表面力 ti（i＝1,2）包含作用于该处的颗粒间接触力、孔隙液体和气体压力、孔隙液相与气相界面上表面张力，如图2所示.式（50）右端积分可展开和计算如下式中是参考颗粒的第k个直接相邻颗粒在接触点处作用于参考颗粒的颗粒间接触力，是Γp上被参考颗粒与第k个直接相邻颗粒间液桥湿化的表面部分，是把Γp 上两个相邻湿化表面部分和分隔开的干表面部分，如图3所示；是沿连接参考颗粒和第k个直接相邻颗粒形心分支向量的面上外法向单位向量，是中点处外法向单位向量.Tk1与Tk2是第k个液桥半月板液面与颗粒表面两交汇处作用于颗粒的孔隙液相与气相间表面张力向量；是该两交汇处颗粒表面的外法向单位向量，如图4所示.将式（47）和式（48），式（50）和式（51）代入式（46）得到它描述了沿接触网络传递的颗粒间接触力.式（52）表明有效Cauchy应力是总应力ji中代表在宏观尺度控制等效连续体变形和破坏的固体骨架应力的那一部分.假定分别为沿颗粒表面和均匀分布的孔隙液体和气体压力，以γ表示表面张力，将式（53）代入式（52）得到表征了因第k个液桥所引起的毛细压力（基质吸力）各向异性.式（55）中αk1，αk2如图4所示;αk是x坐标轴和连接参考颗粒与其第k个直接相邻颗粒形心的分支向量间的夹角，βk是参考颗粒第k个液桥的半填充角，它们能由描述第k个液桥的介观结构的αk1,αk2确定式（56）中和如图5所示，它们能由描述两个相邻液桥介观结构的数据表示，即如图4所示，式（51）中Tk1和Tk2能表示为式中θ表示接触角.利用式（59），能给出式（54）中张量Fγk它表示由第k个液桥引起的表面张力的各向异性.注意到式（42），式（54）提供了基于介观水力--力学的模型有效应力与总Cauchy应力和有效压力张量之间的关系，即可定义有效压力张量如下对比在非饱和多孔连续体模型中所定义的有效Cauchy应力张量的剖分，基于介观水力--力学模型的有效Cauchy应力张量（54）也能类似地剖分为净应力张量和基质毛细压力张量，即式中上式表明，与多孔连续体模型中定义的各向同性基质毛细压力张量不同，基于介观水力--力学模型的基质毛细压力张量为各向异性.引入表示施加于的虚拟孔隙气压分布的系列积分式（52）改写为为计算式（64）中积分和假定沿颗粒表面Γp的孔隙气体压力按下式近似地分布式中，是定义于参考颗粒形心处的梯度，α表示dΓp的单位外法线向量方向.利用式（65），式（64）中积分可表示为将式（53），（55），（60）和（66）代入式（64）得到基于介观水力--力学模型的有效压力张量和基质毛细压力张量可分别定义为式中定义了组构张量式（68）和式（69）表明孔隙液体对等效连续体水力--力学行为的影响由气--液相界面两侧气液相压力差和流--固--气三相交汇处液体表面张力引起的毛细力分布的各向同性与各向异性部分.织构张量与孔隙液桥分布相关，表征了毛细液压效应的各向异性部分.其各向异性依赖于非饱和等效多孔连续体Voronoi胞元的介观结构，即围绕Voronoi胞元中参考颗粒的周边颗粒和液桥分布.为与如式（40）和式（41）所示在非饱和多孔连续体理论中唯象地假定的有效压力张量与基质毛细压力张量比较，考虑式（68）和式（69）所示和的各向同性情况.为此，假定在非饱和Voronoi胞元模型中参考颗粒的所有直接相邻颗粒具有相同的几何与物理性质，以及：（1）围绕参考颗粒的液桥系统的介观结构对称；（2）非饱和Voronoi胞元内孔隙流体压力均匀地分布，即由式（68）和式（69）可得到各向同性情况下基于介观水力--力学模型的有效压力张量与基质毛细压力张量由双联模式液桥计算模型[32-33]，式（71）和式（72）中表面张力可表示为式中 d是参考颗粒与其一直接相邻颗粒间的间隙（正）或重迭（负），ρ是二维液桥的曲率半径.将式（73）代式（71）得到式中定义了各向同性情况下基于介观--水力--力学模型的有效压力及Bishop参数值得注意的是，如式（75）所示，在退化的各向同性情况下基于介观水力--力学模型的有效压力与非饱和多孔连续体理论中唯象地假定的有效压力定义具有相同的形式.然而，式（75）与式（76）给出的Bishop参数是基于由离散颗粒系统到多孔连续体的介--宏观均匀化过程导出，而非唯象地假定.若假定接触角θ＝0，表面张力γ对有效压力张量的影响将消失.Bishop参数将退化为如式（76）和式（77）所示，基于介观水力--力学模型所导出的Bishop参数不仅与在非饱和多孔连续体中定义的孔隙度φ和饱和度Sl有关，且与描述非饱和Voronoi胞元介观结构的参数m,β,d等有关；虽然式（77）没有显式地表示与d 有关，但表征两个直接相邻颗粒间的液桥构形半填充角β依赖于d.给定颗粒体积，孔隙液体饱和度可由定义于Voronoi胞元的孔隙度、液桥数和单个液桥体积计算式中表示位于Voronoi胞元内的单个双联液桥体积一半的Vb/2可计算如下[33] 式中参考颗粒与其一直接相邻颗粒间的间隙（或重迭）量d与Voronoi胞元内颗粒体积、孔隙度、介观结构参数m之间关系为当颗粒体积和m给定，d可视为孔隙度的函数，即将式（79）代入式（78）和利用式（81）可得到式（81）和式（82）表明，当材料参数r,θ和介观结构参数m给定时，半填充角β将是饱和度Sl和孔隙度φ的函数，因而可令并可得到基于介观水力--力学模型的Bishop参数考虑饱和多孔介质中由孔隙液压引起的固体颗粒压缩体积变形，Zienkiewicz等[6,12]引入了饱和多孔介质的广义有效应力概念和相关联的Biot系数.由饱和多孔介质的多孔连续体理论导出的Biot系数为一仅依赖于多孔连续体体积模量与单个固体颗粒的材料体积模量之比的常数.基于介--宏观均匀化过程和Voronoi胞元模型，本文工作表明，Biot系数应是一个变系数，不仅依赖于饱和多孔连续体和单个固体颗粒的体积模量，同时还与两者的当前状态变量、即与饱和多孔连续体当前平均广义有效应力和静水压力作用下单个颗粒的当前压缩体积应变有关.一般情况下，按所建议的基于介观水力--力学模型计算的Biot系数将比按饱和多孔连续体模型所得Biot系数接近于1.在非饱和多孔连续体理论中定义由不可混孔隙液--气相压力加权平均确定的孔隙液--气相混合体压力为有效压力.孔隙流体有效压力在非饱和多孔连续体有效应力上的效应假定为各向同性.孔隙液--气相压力的权系数由依赖于介质饱和度的称为Bishop参数的标量表示，有效压力张量假定为对角线张量.但将原有饱和多孔连续体中仅有单一饱和孔隙流体的有效应力概念推广到非饱和多孔连续体的论据不十分清楚.在多孔连续体理论框架中所建议的各种不同的非饱和多孔介质有效应力定义存在争议.本文基于介观水力--力学模型所导出的非饱和多孔介质孔隙流体有效压力张量表明，一般地它是表征围绕局部材料点液桥分布、液桥毛细力分布等介观结构和介观结构响应的各向异性张量变量.有效压力的各向异性特征意味着孔隙液--气混合体不仅对非饱和多孔介质有效应力的静水应力分量、同时也对它的剪切应力分量有影响. 考虑退化的各向同性情况，基于介观水力--力学模型所导出的非饱和多孔介质孔隙流体有效压力张量的标量形式与在非饱和多孔连续体模型中唯象假定的有效压力标量形式相同.关键区别在于Bishop参数.本文中给出的Bishop参数是基于由离散颗粒系统到多孔连续体的介--宏观均匀化过程导出，而非唯象地假定.所导出的Bishop参数不仅与在非饱和多孔连续体中定义的孔隙度φ和饱和度Sl有关，且与描述非饱和Voronoi胞元介观结构的参数m,β,d等有关.而在非饱和多孔连续体理论中引入的Bishop参数为唯象地假定，同时它不能反映非饱和多孔介质介观结构和介观水力--力学响应对Bishop参数的影响.颗粒材料实际形态较为复杂.本文旨在从颗粒材料多尺度方法出发，提出基于颗粒材料介观结构和响应的饱和与非饱和多孔连续体的有效应力和广义有效应力表达式，揭示其介观机理.为了具体地应用于工程中一些具体的饱和与非饱和颗粒材料的水力--力学行为模拟，考虑颗粒形状的影响是必要的，将是下一步应要开展的工作.【相关文献】1 Terzaghi K von.The shearing resistance of saturated soils.Proc 1st ICSMFE，1936，1:54-562 Biot MA.General theory of three-dimensional consolidation.J Appl Phys，1941，12:155-1643 Biot MA.Theory of propagation of elastic waves in a flui saturated porous solid.J Acoust Soc Am，1956，28:168-1914 Biot MA.Mechanics of deformation and acoustic propagation in porous media.J Appl Phys，1962，33:1482-14985 Biot MA.Generalized theory of acoustic propagation in porous media.J Acoust Soc of America，1962，34:1254-12646 Zienkiewicz OC，Shiomi T.Dynamic behavior of saturated porous media：the generalized Boit formulation and its numerical solution. Int J Numer Anal Meth Geomech，1984，8:71-967 Kytomaa H，Kataja M，Timonen J.On the effect of pore pressure on the isotropic behavior of saturated porous media.J Appl Phys，1997，24:7148-71528 Bishop AW.The principle of effective stress.Teknisk Ukeblad，1959，106（39）:859-8639 Bishop AW，Blight GE.Some aspects of effective stress in saturated and partly saturated soils.G’eotechnique，1963，13（3）:177-19710 Skempton AW.Effective stress in soils，concrete and rock.//Proc. Conf.on Pore Pressure and Suction in Soils，Butterworth，1960. 4-1611 Li XK，Zienkiewicz OC.Multiphase fl w in deforming porous media and finit element put.Struct.，1992，45（2）:211-227。