重力式挡土墙地震响应分析

重力式挡土墙地震响应分析

1小应变条件下的耦合方法在挡土墙地震模拟中的应用基于颗粒介质力学理论，从薄视力学的角度研究了砂、土和结构之间的相互作用，成为一个新的焦点。由于计算机性能的限制,用离散颗粒模拟砂土时,颗粒数目过多的情况下不可避免地遇到运算无法执行或计算过于缓慢,因此出现了土体离散-连续耦合方法以减少颗粒数量,只在关心的区域进行细观层次模拟。Cai等在临近隧道开挖的局部土体区域用离散颗粒模拟,其他区域用连续模型模拟,但连续模型是按照静力条件设定的。对于挡土墙和砂土的相互作用,已有基于离散颗粒方法的挡土墙主动和被动土压力模拟。挡土墙的地震响应动力分析中,若土体全部采用离散颗粒模拟,运算时间会过长,因此在分析中尝试用耦合方法进行模拟。耦合方法涉及颗粒材料宏细观参数联系的问题,对二维问题已有的模拟方法有:双轴数值试验模拟静模量和循环双轴试验模拟动力性质。但对于动力耦合分析来说,小应变情况下的颗粒集合体动模量的模拟方法仍不完善。本文使用离散-连续耦合动力模型,模拟了地震作用下挡墙的动力响应。对靠近结构物的近场土体采用离散元方法进行精细化分析,砂土与地下结构之间满足非连续动力接触关系,振动过程中砂土对结构物的作用直接以接触力的形式作用于结构物上;远离结构的土体采用连续模型模拟。耦合方法的两个核心问题:(1)保证离散-连续土体在边界上的连续性,这里保证边界处速度和力的连续性,本文提出新的边界处耦合力提取方法;(2)模拟中细观区域颗粒的宏观性质需与连续土体模型一致,为完善离散颗粒宏细观参数联系的模拟方法,首次给出颗粒材料的自振柱模拟,获得了小应变下颗粒集合体动模量随应变衰减的关系。对于结构与砂土的相互作用,模拟中采用与离散-连续土体边界耦合一样的方法,但结构-砂土边界处耦合力直接采用颗粒与结构的接触力。本文使用的耦合方法易于实现,细观土体区域采用颗粒流软件PFC2D模拟,远离结构的土体用有限差分软件FLAC2D模拟,每个循环时步中基于两个软件的数据传输实现宏细观土体的耦合。为避免有限差分算法中结构与土体模量差别过大而导致的时步过小问题,FLAC2D软件只模拟了连续土体区域,对于挡土墙编写了梁单元受碰撞的动力有限元程序嵌入颗粒流软件PFC2D中。通过模拟文献中重力式挡土墙地震响应的离心机试验,说明了耦合的实现过程。耦合方法从细观层次上模拟了结构受砂土作用的动力过程,并有效减少了颗粒数目,节省运算时间。2小应变下颗粒单体的动模量模拟方法离散-连续耦合方法中,必须保证颗粒集合体的宏观性质与连续体模型中的参数一致,这涉及颗粒集合体宏观性质的模拟。砂这类离散颗粒的动模量随应变的增加而衰减,在室内试验中对于应变小于10-4情况下的动模量一般用共振柱试验获得,而对于应变大于10-4动模量一般用动三轴试验获得。对于颗粒集合体的循环双轴模拟方法文献中已有叙述,但对适用于小应变下颗粒集合体的动模量模拟方法尚未见文献报道。连续土体模型的控制方程中,需代入完整的颗粒集合体动模量与应变的关系曲线,因此建立小应变下颗粒集合体的动模量模拟方法是必须的。应当指出,用离散元方法直接模拟室内共振柱试验是非常耗时的,原因在于共振柱方法是采用频率扫描的方法,确定圆柱试样的一阶固有频率,再结合试样尺寸和密度,根据弹性波在一维杆中的传播公式推求试样的动模量。对于单一频率的动载荷,用离散元方法模拟试样的振幅响应较为耗时。实际模拟时需对某一带宽的频率进行扫描,还需模拟不同应变下的动模量,由于计算机性能的限制,计算时间过长,难以实现。1981-1982年间,在共振柱试验的基础上发展了自振柱试验。这种新的试验方法是预先对圆柱试样施加一给定应变,释放应变限制后让试样按其固有特性作自由振动,根据振动的自由衰减曲线获得试样的固有频率,再推求试样的动模量。由于不需要频率扫描,用离散元模拟自振柱试验耗时远小于模拟共振柱试验。这里给出了自振柱试验的离散元模拟方法,结合动三轴模拟结果,得到完整的动模量随应变变化关系曲线。2.1砂级配线的模拟试验中所用的砂为Nevada-120砂,砂土性质见表1,砂土级配见表2。在离散元数值模拟中,由于计算机性能的限制,即使采用宏细观耦合方法,仍须采用将实际砂土粒径放大的办法减少颗粒数量。数值模拟试样和离心机试验所采用的Nevada-120砂的级配曲线如图1所示。模拟时,需使数值试样的宏观性质对应于实际试验中砂的宏观性质。文献虽然是基于干砂的离心机试验,但没有提到干砂试样的三轴试验结果,仅提到饱和试样的三轴试验结果。文献通过对双轴数值试样施加等幅应变,并控制加荷过程中试样面积不变的方法模拟砂土的液化性质。这里采用其加荷过程中试样面积不变的方法,模拟饱和砂土的双轴试验,不断调整颗粒的细观参数,使数值试样在饱和双轴试验中摩擦角与实际饱和砂一致。数值模拟需确定砂土的细观参数有:颗粒孔隙率e、法向接触刚度kn、切向接触刚度ks、颗粒摩擦系数fc和颗粒密度ρs,颗粒的细观参数如表3所示。2.2分散自扰的模拟2.2.1自振柱体系的接触问题自振柱模拟本质上是模拟一端固定、一端自由的一维黏弹性杆,受一初始应变后自由衰减振动的过程。对于如何从自振柱扭转试验中推求动模量,文献予以详细叙述。这里采用二维模拟,基于圆柱扭转推求的公式不适用,因此需得到自振柱压缩试验中固有频率与动模量的关系。一维纵波在杆中的微分方程如式(1)~(3)所示:式中:u(x,t)为杆中坐标为x、时刻t对应的位移;U(x)为位移函数;T(t)为时间函数。式中:A和B为由边界条件确定的参数;ρ为柱体密度;E为柱体模量;ωn为柱体n阶固有频率。式中:n=μωn2/(2E);;μ为阻尼系数;ϕ为由初始条件确定的参数。对于如图2所示的自振柱体系,在x=0端固定,在自由端x=L处加一质量为m的块体,则在x=L处颗粒集合体所受力如式(4)所示,块体的平衡方程如式(5)所示。式中:N为x=L处杆所受外力,杆应变ε(L,t)=∂u(x,t)/∂x|x=L;a为柱体面积;Eε(L,t)为弹性应力;μ∂ε(L,t)/∂tu(x,t)为因黏性产生的应力。式中:m为块体的质量。在x=L处杆受力和质量块受力相等,由式(4)和(5)得将式(1)~(3)代入(6),结合边界条件,可得模拟时,x=L处施加一指定位移,从x=L处的衰减曲线可以得到一阶频率ω1,将其代入式(7),可求解得到颗粒集合体在一指定应变下的动模量。2.2.2颗粒不平衡和围压的共同作用分析采用Newmark积分编写质量块受力运动的程序嵌入离散元软件PFC中,质量块的动力方程形式如式(5)所示。模拟的步骤是:(1)先生成颗粒试样,对颗粒集合体用刚性墙在四边施加给定围压,如图2所示(但不包括质量块,这里质量块是虚拟的)。(2)加围压至颗粒平衡后,固定底边墙(即杆件固定端),对顶墙施加一指定位移(即在杆自由端施加位移,给杆件一初始应变),在这个过程中保持两边侧墙的伺服程序,即保持试样侧边围压不变。(3)指定数值试样顶部墙的运动速度,速度是按照质量块受杆件顶端作用力(即试样顶部墙受颗粒不平衡力)和围压的共同作用时,采用Newmark数值积分求解得到,通过记录顶部墙的位移变化,可得到试样自由衰减振动曲线,从中得到试样固有频率,进而求解颗粒集合体的动模量。数值模拟时,对质量块的运动求解中,需将Newmark积分的步长与颗粒流软件PFC中的步长设为一致。由于是松砂,第3步模拟时,若侧墙仍保持伺服围压算法将对试样造成过大扰动,这时侧墙固定。质量理论中质量块对模拟得到的频率没有影响,但数值试验表明,较大的块体质量可以使得到的自由振动衰减曲线更加平滑,因此更易从中求得频率。数值模拟试样的颗粒细观参数在土样模型介绍中已有叙述,这里模拟时试样的初始高度为0.2m,底部宽度为0.1m。模拟了围压为40、80、160、320kPa时不同初始应变下颗粒集合体的动模量。图3给出初始应变接近10-5,各围压下试样柱体自由衰减振动时,轴向应变随时间的变化关系,应变数值为正时表示试样受压缩。从图中可以看出,当围压大于一定值时,试样柱体表现简谐衰减振动,这时可以从中确定试样的固有频率。但是,当围压过小时,例如图中显示围压为40kPa时,有可能无法从振动曲线中获得频率,这是因为松砂在低围压下过于松散,只能单纯受压缩,以致无法将试样当成一个弹性体。表4~7给出了不同围压、不同初始应变下求得的颗粒集合体动模量。从总的趋势上说,随初始应变加大,其对应的动模量会减小。但是由于是松砂,加大初始应变后会造成颗粒试样变密,而使对应的动模量变大。因此,当初始应变增大,对应的动模量表现出一定的波动性,但总的趋势是减小。2.3小应变循环载荷模拟应变大于10-4时,离散颗粒集合体的动模量用动三轴获取更为合适,基于小应变理论的共振柱或自振柱方法已不合适。采用二维循环双轴数值模拟,可以模拟大应变时颗粒集合体的动模量。双轴循环模拟在文献中已有叙述,这里简要叙述采用的循环双轴试验实现步骤:(1)如同自振柱模拟时先生成颗粒试样,加至预定围压,并循环至平衡,如图2所示(不包括质量块);(2)保持侧墙的伺服围压,对颗粒集合体的底部墙和顶部墙施加方向相反的等幅正弦位移,获得等应变循环载荷的应力-应变关系曲线,由曲线求得试样的割线模量。图4给出围压为80kPa、应变为0.0005时试样在循环载荷下的应力-应变关系,这里均采用第2个循环圈来求得割线模量。表8给出不同围压、不同等幅应变的循环载荷模拟得到的动模量。当围压过小而应变过大时,砂过于疏松,导致施加等幅循环载荷时颗粒与上下墙分离;当围压过小且应变过小时,上下墙与颗粒接触力过小。以上两种情况都难以模拟得到动模量。由上述第2种情况可知,对于松砂来说,小应变情况下的动模量用循环双轴试验模拟是难以实现的。但是这仅对于松砂而言,对于密砂,通过数值试验可以发现用循环双轴试验可得到与自振柱模拟相似甚至更高的动模量。2.4量随应变变化的拟合函数动模量E与应变ε关系的拟合曲线如图5所示。动模量随应变变化的拟合曲线形式有多种,这里采用的拟合曲线形式见文献,拟合函数形式如式(8)所示。式中:C1、C2、C3、C4和5C均为拟合参数。3宏细观察联合算法3.1连续模型耦合在土体离散-连续区域边界上,需保证连续性。这里保证速度和力的连续性。图6显示了一个土体连续模型单元与离散颗粒接触的情况。连续模型单元与颗粒的交界面上承受颗粒碰撞,因此这个交界面可以作为连续模型的力边界;同时运动的交界面约束颗粒,交界面可以作为颗粒的速度边界。因此土体离散-连续区域动力耦合算法可以表述为:(1)在动力分析的每一时步中,用离散颗粒模型计算边界墙应受到的力并将此数据传输至连续模型中,接收连续模型中边界墙两个节点的速度数据并由此给定边界墙的速度;(2)在连续模型中,计算边界墙两个节点的速度并将此数据传输至离散模型中,接收离散模型中边界墙受力数据并将此力施加至连续模型上。不同于已有耦合算法,这里给出确定连续模型在耦合边界上受力的新方法。已有的算法是将颗粒与边界墙的接触力直接施加至连续模型上,但是砂土连续模型实际上是对离散颗粒的平均近似,因此为更好地在边界线附近保持力的连续性,在紧邻边界墙的位置,设一过渡单元,如图6中边界墙与虚线包围区域所示。通过求解过渡单元中颗粒的平均应力来确定连续模型在耦合边界上须施加的力。设过渡单元中颗粒数为N,颗粒的孔隙率为P,颗粒的应力分量为σij(i=1,2;j=1,2),则过渡单元处的平均应力为式中:为过渡单元的平均应力分量。假定过渡单元处的平均应力与边界墙处应力一致,则可根据式(9)确定施加至边界墙上的应力。3.2梁单元的运动方程编写了梁单元受颗粒碰撞的动力有限元程序,将其嵌入离散元软件PFC中。动力有限元算法中不考虑结构的阻尼,其动力方程为式中:{x}为结构各节点的位移向量;[M]为结构的质量矩阵;[K]为结构的刚度矩阵;{f}为外力载荷向量,包括重力和颗粒与结构的接触力。梁单元节点上的自由度包括转动与两个方向上的平动,有限元形函数、Newmark积分算法、质量与刚度矩阵见文献[11-12,14-15]。图7给出了结构与颗粒相互作用图,图中的4个梁单元组成了一个矩形框结构,每一个梁单元对应离散元软件PFC中的墙,这四面墙的运动速度由编写的动力有限元程序求解得到。结构单元与颗粒相互作用的耦合算法同土体离散-连续耦合算法,但这里不设置过渡单元,颗粒碰撞在边界墙上的力直接作用在结构上。在离散元软件中编写结构动力有限元程序的原因在于:土体离散-连续区域耦合中连续区域是在有限差分软件FLAC中模拟的,但在有限差分算法中,结构与土体模量差别巨大会造成运算时步急剧缩小,总运算时间增加。直接将结构动力程序嵌入离散元软件,满足土体耦合算法节省计算时间的目的。4数值模拟和结果分析4.1模型试验概况文献进行了重力式挡土墙地震响应的离心机试验,试验的主要目的是研究挡土墙与砂土在交界面的相互作用,试验假定为平面应变状态。选取模拟的一组试验中,砂土为松砂。挡墙用一长方体来模拟,其内布置加速度和力传感器,挡土墙近似为一刚体,长为100mm,宽为100mm,高为125mm。模型箱的底部和侧边是相对固定的,试验中模型箱的底部施加频率为80Hz的加速度载荷。试验模型图如图8所示。4.2计算边界条件下的平衡数值模拟中挡墙和模型箱中土体尺寸依据图8所示的离心机试验设定。模拟步骤:(1)生成颗粒与挡土墙结构后,让其与连续土体模型在非耦合的情况下,在1倍重力下相互作用至平衡;(2)让挡墙结构、离散颗粒和连续土体模型耦合,先在1倍重力下循环计算至平衡,再在28倍重力下计算至平衡;(3)在连续模型土体区域的底边和两侧边施加频率为80Hz、最大加速度为6.3倍重力加速度的正弦波载荷。挡墙简化为如图7所示的由梁单元组成的结构,离散颗粒在离散元软件PFC2D中模拟,连续土体区域在有限差分软件FLAC2D中模拟。连续土体模型中采用摩尔-库仑模型,模量随应变变化的关系采用图5中围压为80kPa的拟合曲线。离散元软件PFC2D中颗粒用厚度为1m的圆盘来表示,但实际试验中挡墙厚度只有0.1m,因此模拟时将离散元模型与连续模型中的密度都设为原来的0.1倍。图9给出挡墙、离散颗粒和连续土体模型耦合下,在1倍重力下的平衡结果。土体离散-连续模型中,搜索颗粒是否在过渡单元中采用文献中的算法,但每一循环中将所有颗粒搜索一遍将会非常耗时,因此实际模拟时每1000循环步全局搜索一次。4.3土体应力分量的测量圈模拟与实际试验的主要差别在于:(1)模拟时底部施加光滑的正弦波载荷,但实际试验中测得的模型箱底部加速度不如模拟时光滑规则;(2)在细观颗粒尺度上,虽然在颗粒级配、颗粒集合体模量上模拟使用的颗粒尽量与试验接近,模拟使用的颗粒与试验中实际砂颗粒还是有不小区别。实际砂颗粒有棱角,而模拟使用的颗粒截面是光滑规则的圆形;颗粒放大后,离散颗粒在动载荷下孔隙率的变化性质会发生变化,造成颗粒集合体在重力和挡墙作用下的压缩性与实际试验有差别。图10显示了在动力载荷施加至0.2s时,挡墙和离散颗粒的相对运动。试验中挡墙的转动方向与模拟得到的结果不同,这里认为主要是上面提及的模拟与实际试验的差别造成的,特别是加大重力后,模拟试样和实际砂试样的压缩性的差别被放大,例如若实际试样中砂试样的压缩性大于模拟试样,加大重力后动载荷下,挡墙右下侧的土体沉降更易大于挡墙左下侧土体,这种情况下挡墙的转动方向会与模拟得到的结果相反。图11和图12分别显示挡墙右侧测量圈中土体应力分量在动载荷下的时程曲线。测量圈的设置如图9所示,测量圈可以在细观尺度上分析记录圈内颗粒平均特性。模拟时每隔一定循环步在挡墙附近设一测量圈,因每一次设置的测量圈与挡墙有重叠,这里记录的颗粒平均应力分量要乘以一个大于1的系数才是实际离散模型中颗粒的平均分量