饱和多孔介质土动力学研究进展

饱和多孔介质土动力学研究进展

饱和土是工程中常见的土壤，是流固定相介质。饱和土壤肥力特性远比单相固体介质复杂，因此对土壤肥力特性的研究主要是一种简单的单株。结果表明，采用两相多孔介质模型比使用单一固体介质模型更合理，土壤肥力的动态反应更均匀、更真实。关于饱和土壤肥力方程的研究始于20世纪50年代。bio首先建立了液体饱和多孔介质波的传播理论，这是后来各种研究中关于饱和多孔介质动态特性和动态分析的基础。bio理论的本质是将连续介质力学应用于液体饱和多孔介质的两层材料体系中的轨迹效应，并考虑了方程和固体框架之间的应变性关系及其运动。多孔介质力学的理论和应用已经成为岩石工程、地震工程、工程振动、地球物理等科学领域的一项重要课题，也是土壤振动和土球物理工程的一个重要分支。以下是饱和多孔介质力学的基本方程、饱和土壤中弹性波的传播特性，以及饱和土壤肥力分析的时间值解和地震反应分析等各个方面的发展历史和进步。1动力方程的推导自Biot提出描述饱和多孔介质动力特性的基本方程以来,此基本方程已经广泛地应用于土动力学分析中.继Biot之后,国内外很多学者从不同角度对该问题作过研究,但所提出的表达式有些混乱,容易让人产生误解.事实上,多孔介质的动力方程完全可以通过混合物的均匀化理论建立起来,整个推导过程是从有效应力原理、应力应变关系、质量守恒方程以及动量守恒方程出发,这样的推导方法不仅理论上比较严密,且容易推广到多相介质.1.1ter-qp的应力应变关系根据混合物理论的空间平均方法,饱和多孔介质土的总应力(不考虑孔隙液体的粘滞性)表示为σij=(1-n)σsij-nδijpw(1)式中:σsij为作用于固相的应力;pw为孔隙水(流体)压力;δij为Kronecker符号;n为孔隙率.根据Terzaghi有效应力原理可以得到固相应力σsij与有效应力σ′ij之间的关系式为(1-n)σsij=σ′ij-(1-n)δijpw(2)如果考虑土颗粒本身由于孔隙水压力作用所引起的压缩性,基于有效应力的应力应变关系可以表示为(假设土骨架是各向同性理想弹性体)σ′ij=2μεij+λδijεkk+(1-α)δijpw(3)式中:α=1-KT/Ks为Biot参数,KT和Ks分别为土骨架和土颗粒的体积弹性模型,ΚΤ=λ+23μ,μ和λ为弹性理论中的Lame常数;εij表示为土骨架的应变.同样,固相应力σsij和总应力σij可分别表示为(1-n)σsij=2μεij+λδijεkk-(α-n)δijpw(4)σij=2μεij+λδijεkk-αδijρw(5)1.2孔隙水体变模量的计算固相和液相的质量守恒方程可以分别表示为对于固相∂[(1-n)ρs]∂t+(1-n)ρs˙usi,i=0(6)对于液相∂(nρw)∂t+nρw˙uwi,i=0(7)式中:ρs和ρw分别为固相(土颗粒)和液相(孔隙水)的密度;usi和uwi分别为固相和液相的绝对位移.考虑到土颗粒和孔隙水本身的压缩性,由式(6)和式(7)相加可得-(1-n)˙σsii3Κs+n˙ρwΚw+(1-n)˙usi,i+n˙uwi,i=0(8)将式(4)代入式(8)可以得到˙pwΜ+(α-n)˙usi,i+n˙uwi,i=0(9)式中:Kw为孔隙水的体变模量;1Μ=α-nΚs+nΚw.引入e=usi,i和ε=uwi,i,由式(9)和式(4)可以得到-npw=Qe+Rε(10)(1-n)σsij=2μεij+δij(Ae+Qε)(11)其中:R=n2M;A=λ+(α-n)2M;Q=n(α-n)M.式(10)和式(11)与Biot所提出的应力应变关系式是完全一致的.1.3土体密度的测试固相和液相的动量守恒方程可以分别表示为对于固相-(1-n)ρs⋅⋅usi+(1-n)σsij,j+nRwi+(1-n)ρsbi=0(12)对于液相-nρw⋅⋅uwi-npw,i-nRwi+nρwbi=0(13)式中:bi为单位质量的体积力.有时更方便的形式是采用整个土体的动量守恒方程,即以总应力表示的动量守恒方程代替固相的动量守恒方程,由式(12)与式(13)相加可以得到σij,j-(1-n)ρs⋅⋅usi-nρw⋅⋅uwi+ρbi=0(14)以上各式中,Rwi表示固相和液相之间的耦合作用,ρ=(1-n)ρs+nρw为土体的密度.利用Darcy定律可以得到kρwgRwi=˙wwi(15)式中:k为渗透系数;˙wwi为液相的相对速度(相对于固相而言),定义如下:˙wwi=n(˙uwi-˙usi)(16)1.4美—饱和多孔介质土动力学基本方程的建立饱和多孔介质土动力学基本方程可以分别用us-uw-pw形式、us-uw形式、us-ww形式以及us-pw形式来描述.为了书写方便,下面分别用u代替us,U代替uw,ρf代替ρw,pf代替pw以及w代替ww,并且不考虑重力的影响.1.4.1biol相对粘滞性和渗透系数由式(14),(5),(13),(15),(16)以及式(9)可以得到μΔ2u+(λ+μ)grad(divu)-αgradpf=(1-n)ρsu⋅⋅+nρfU⋅⋅(17)-ngradpf+b(u˙-U˙)=nρfU⋅⋅(18)ndivU˙+(α-n)divu˙+p˙fΜ=0(19)式中:b=n2ρfgk=n2ηkp,其中η和kp分别为Biot定义的流体粘滞性和渗透系数.1.4.2u+[++-n2grad质量分数的计算由式(12),(4),(13),(15),(16)以及式(9)可以得到μΔ2u+[λ+μ+(α-n)2Μ]grad(divu)+n(α-n)Μgrad(divU)=(1-n)ρsu⋅⋅+b(u˙-U˙)(20)n(α-n)Μgrad(divu)+n2Μgrad(divU)=nρfU⋅⋅-b(u˙-U˙)(21)1.4.3u+++22grad质量分数的计算将式(16)代入式(17)、(18)和(19)得出μΔ2u+(λ+μ+α2Μ)grad(divu)+αΜgrad(divw)=ρu⋅⋅+ρfw⋅⋅(22)αΜgrad(divu)+Μgrad(divw)=ρfu⋅⋅+ρfnw⋅⋅+ρfgkw˙(23)1.4.4pf的推导计算结果表明,对于像地震荷载这样非高频的情况,可以采用以固相的位移u和孔隙水压力pf为基本变量的所谓u-pf形式,这种u-pf形式是由Zienkiewicz等首先提出的.具体的推导过程如下:将ρfw⋅⋅和ρfnw⋅⋅项略去不计,并消去w,则有μΔ2u+(λ+μ)grad(divu)-αgradpf=ρu⋅⋅(24)αgrad(divu˙)-k—Δ2pf+p˙fΜ-k—ρfdivu⋅⋅=0(25)式中:k—=k/ρfg.2表中弹性波的传播特征2.1饱和土介质中压缩波种类的分析在文献中,采用u-U形式推导了饱和土中的体波波速表达式,并假设土颗粒是不可压缩的.这里采用更常用的u-w形式,并进一步考虑土颗粒的压缩性.仿照文献的推导方法,通过引入势函数,可以得到下列表达式1vΡ121vΡ22=ρ(ρfn-iρfgkω)-ρf2(λ+2μ)Μ(26)1vΡ12+1vΡ22=(λ+3μ+α2Μ)(ρfn-iρfgkω)-2αΜρf+ρΜ(λ+2μ)Μ(27)1vs2=1vs02+ibnρfμ(ib-nρfω)(28)式中:vS0=μ(1-n)ρs.由以上各式可知,饱和土介质中存在2种压缩波P1和P2(波速用vP1和vP2表示),以及一种剪切波S(波速用vS表示).3种体波均与振动频率ω有关,因此都具有弥散性.下面给出2种极限情况下的体波:2.1.1s不同体积的压缩波此时k→∞,有b→0,由式(26)、(27)得1vΡ121vΡ22=ρρfn-ρf2(λ+2μ)Μ‚1vΡ12+1vΡ22=(λ+3μ+α2Μ)ρfn-2αΜρf+ρΜ(λ+2μ)Μ(29)vS=μ(1-n)ρs(30)若不计土颗粒压缩,取α=1,M=Kf/n.这里可以举一个例子,取μ=2.61×107Pa,λ=2.61×107Pa,n=0.27,Kf=2.0×109Pa,ρs=2.65×103kg·m-3,ρf=1.0×103kg·m-3,求解方程(29)可以得到:vP1=2100.15m·s-1,vP2=135.84m·s-1,而只在水中传播的速度vf=Κfρf=1414.21m·s-1,干燥条件下的压缩波vΡ0=λ+2μ(1-n)ρs=201.19m·s-1,由此可见,第一压缩波的传播速度大于只在水中传播的速度,而第二压缩波的速度低于干燥条件下的压缩波.总之,在孔隙流体可以自由流动的情况下,3种体波均与振动频率ω无关,不具有弥散性.2.1.2饱和土中弹性体波的散射性此时k→0,有b→∞,在这种条件下只存在一个P波和一个S波,它们的波速分别为vΡ1=λ+2μ+a2Μρ‚vS1=μρ(31)文献对饱和土中弹性体波的弥散性作了深入的探讨,讨论了波的频率、饱和土的孔隙率以及渗透性对各种波的速度和波速衰减的影响.2.2饱和土中瑞利波的计算Jones和Chiang研究了饱和土中的瑞利波问题,但在建立特征方程时,在压缩波势函数中只考虑一种压缩波,因而势函数不是问题的通解,导致特征方程有误.Tajuddin建立了考虑2种压缩波的瑞利波特征方程,但算例中的结果是错误的.夏唐代等以及陈龙珠等基于文献所提出的饱和土波动理论,推导出饱和土中瑞利波特征方程,并进一步讨论了半空间饱和土中瑞利波的弥散特性,值得注意的是文献所提出的饱和土波动理论与Biot理论是有所不同的.3zoikiefficir形式的数值分析方法地震反应分析无疑是饱和多孔介质土动力学理论研究的一个重要应用.对饱和土体地震液化和变形分析常常采用时域有限单元法分析技术,Biot动力方程的第一个时域有限元数值解法是Ghaboussi和Wilson提出的,当时以固相的位移和相对于固相的液相w位移为基本变量.计算结果表明像地震荷载这样非高频的情况下可以采用以固相位移u和孔隙水压力p为基本变量的所谓u-p形式.这种u-p形式由Zienkiewicz等首先提出,此后进一步推广到非线性大变形问题的分析中.另一种比较常用的形式是Prevost以及Zienkiewicz和Shiomi所用的u-U形式(U为液相的绝对位移),另外Zienkiewicz和Shiomi还给出完整的u-U-p方程的有限元形式.由于u-p方程相对比较简单,因此可以说是目前最常用的.从数值方法角度来看,主要有以下三方面的问题:其一,所谓的不可压缩问题,由于孔隙流体的压缩性通常是忽略不计的,当渗透系数比较小时,u-p混合有限元形式将要求对于位移u的插值形函数至少比对于孔隙压力p的插值形函数高一阶,否则会造成严重的数值不稳定,当然这种不同阶的插值势必给有限单元法编程带来麻烦,也使一些简单而又实用的单元形式不能被采用,因此有必要寻找一种新的办法使得任何单元形式都适应.其二,求解u-p形式的时域有限单元法往往采用完全隐式的时间积分格式,虽然这种积分格式是无条件稳定的,但需要形成包含刚度矩阵、质量矩阵、耦合矩阵、渗透矩阵以及孔隙水压缩矩阵所组成的总体矩阵,这样当分析的问题复杂或自由度较多时,其计算工作量将很大,因而极大地限制了它的应用与发展.为克服上述缺点,有必要寻找一种交替解法,这种新的解法将不需要形成总体矩阵,而且还必须是无条件稳定的,否则,当考虑材料非线性时,解法的稳定性将很难得到保证.其三,饱和土工构筑物在地震力作用下可能会发生液化而造成较大的位移,甚至导致构筑物倒塌,1971年2月9日美国加州大地震中严重破坏的LowerSanFernando土坝就是一个最典型的例子.针对这种大位移的有限元分析技术不仅要采用大变形的分析方法,而且还要结合自适应网格技术来跟踪变形的发展.从数值分析方法的验证来看,验证无非有2种途径.第一条途径是采用离心机试验结果;另外一条途径则是采用历史上受地震破坏的实例.3.1显式—饱和土体动力学u-p方程的不可压缩问题目前常用的求解饱和土动力学u-p方程的时域有限元解法主要采用无条件稳定的完全隐式整体解法,但这种解法存在着一个明显的缺陷,当渗透系数较小时(一般来说孔隙流体的压缩性可以忽略不计),联合方程组的整体矩阵会出现一个数值接近于零的子矩阵,这样的一个子矩阵的存在将要求对位移u和孔压p采用非等价的形函数,即要求位移u的形函数的价数比孔压p的形函数至少高一价,这样会带来很多数值分析上的困难.克服这个困难的方法有2种:第一种方法是采用一种基于矢量分解的交替解法,交替解法可以通过显示—隐式时间积分格式(半显式格式)或者隐式—隐式时间积分格式来实现,数值分析结果表明有效的方法必须是先对基本微分方程进行时间离散,然后再进行空间离散.值得一提的是这种新的显式—隐式交替解法(半显式解法)的稳定性比Zienkiewicz等提出的半显式解法有着很大的改善,其时间步长的控制与土的渗透系数以及孔隙水的压缩系数无关,这样大大改善了数值解法的有效性;所提出的隐式—隐式交替解法已被证明是无条件稳定的,使得数值计算的时间步长的选择更加自由,解决了所谓的不可压缩问题,同时大大提高了计算效率,计算结果尤其是孔压值显示出良好的稳定性.上述交替解法对于线弹性问题是非常有效的,但当用于弹塑性分析时,由于其时间积分的稳定性无法得到严格的证明,可能会造成其稳定性无法得到保证.为了解决这个问题,文献又进一步提出了一种整体解法,这种新的整体解法是通过对有限单元法空间离散前的微分方程直接修正的方法来实现,通过修改基本方程来有效解决u-p形式存在的不可压缩问题;同时,其时间积分格式的稳定性与目前常用的完全隐式整体解法完全一致.3.2交替方向乘子法正如前面已提到的,饱和土动力分析u-p形式通常采用时域有限元的直接解法求解,但当涉及到自适应网格或三维问题时直接解法太费时,因此交替解法是一种应该考虑的选择.交替解法存在的最大问题是时间积分的稳定性.Zienkiewicz等曾提出一种无条件稳定的交替解法,但只适用于对孔压采用非连续形函数的情况,这样就大大地限制了它的应用.文献所提出的一种能保证时间积分无条件稳定的交替解法,不仅其稳定性可以得到严格证明,而且适用于任何单元形式.通过对著名的美国LowerSanFernando土坝地震反应的自适应网格有限单元法分析,表明所提出的交替解法将大大减少计算的工作量.可以预见,对于自由度很多的三维分析所提出的交替解法将显得更加有效.3.3基于自适应网格的地震反应1971年2月9日的美国加州大地震造成了LowerSanFernando土坝的严重破坏,上游边坡面出现了大规模的滑坡.这一近于灾难性的破坏引发了对此坝广泛的现场调查和理论分析.从现有的分析结果可以看出,常规的有限元技术已经很难合理地描述这种由于动力应变局部化引起的渐进破坏过程,因此有必要引入一种新的有限元技术.文献提出了一种基于类似位移梯度的误差指示的自适应网格有限元技术,这种自适应网格技术不仅可以对网格变形比较大的区域进行网格重新形成,而且还可以将单元沿剪切带发展方向拉伸单元.通过对LowerSanFernando土坝的地震反应分析,合理地解释了地震以后由于大位移而造成的坝体滑坡现象,通过自适应分析可以清楚地观测到地震后位移的增量,证明了这种新的有限元技术是有效的,为解决饱和土体