基于有限元法的E-K模型参数敏感性分析
龚羊庆
江西省水利科学研究院
吴海真
江西省水利科学研究院
黄英
昆明理工大学
模型参数的确定是数值分析过程中的重要环节,合理把握参数的选取原则,可大大提高结构分析的可靠度以及参数反分析的计算效率。本文借助有限单元法考察了E-K模型参数变化对一个均质堤坝的位移和应力值的影响程度。通过分析明确了E-K模型中各参数对计算结果的敏感程度,同时提出了模型参数的选取原则。
1 引言
土体力学参数的取值因其复杂的力学特性而存在较大差异。它不仅随荷载大小而异,而且还与加载应力路径有关。此外,土体还会发生屈服,产生不可恢复的塑性变形[1]。诸多因素致使用土体试验参数作为计算输入参数进行数值分析时,所得结果往往与实际情况会有一定的误差,且不同程度地阻碍了数值计算在岩土工程中的应用。针对这一情况,国内外研究工作者开始对计算理论的反演问题开展研究,力图依据在工程现场获取的信息,尤其是根据位移量测信息反演,确定各类未知参数(如弹性模量和体积模量等)。因此,在土工计算中最棘手的问题是如何合理并切合实际地选取介质材料的物理力学参数。
邓肯-张模型在国内外广泛使用近30年,大量的试验成果表明,由于取样制样、试验仪器、试验方法、试验人员操作熟练程度、整理分析资料等诸多因素,使E-K模型中的7个参数(c、φ、k、n、Rf、kb、m)变化较大(其中k值可成倍甚至数量级的差别),用于计算所得结果差别也较大。所以有必要对这些参数进行敏感性分析,以正确掌握这些参数对计算结果的影响,这将有助于把握参数的选取原则,提高结构分析的可靠度以及参数反演分析的计算效率等。
2 敏感性分析方法
邓肯等人曾对E-K模型中的7个参数作了初步讨论,并给出了几种不同类型土参数的取值范围(如表1所示)[2]。由于表1中的数值变化范围较大,不同取值对计算结果的影响没作进一步讨论。本文将借助分析土工结构中应力和变形问题的有限元软件Geo-Studio[3]来考察参数变化对一个均质堤坝(见图1)的位移和应力的影响程度。考察某一参数时,该参数值上下波动,其余6个参数保持试验值。
表1 邓肯-张模型参数取值范围
注 kur为土体回弹模量的试验常数。
图1 堤坝有限单元网格
本文发表于2011年。
3 结果分析及讨论
图1为模拟一座堤坝在填筑过程中的应力和变形分析的有限单元网格图。堤坝填筑土为黏性土,高15m,上下游坡度均为1:2.0。堤坝地基约为20m的砂质黏土层。有关该实例的有限单元分析信息见表2和表3。按表3中堤坝填筑土的试验参数进行有限元分析,得出堤坝的水平位移、垂直位移、竖向总应力、水平方向总应力以及坝轴线处竖向位移、竖向总应力和水平向总应力随高度的变化情况如图2~图5所示。
表2 堤坝有限单元分析的有关信息
续表
表3 堤坝填筑土的物理力学参数
图2 水平位移等值线图(单位:m)
图3 竖向总应力等值线图(单位:kPa)
图4 垂直位移等值线图(单位:m)
由图2~图5可见,坝体最大竖向位移位于坝体中部区域;水平位移发生在上下游坝脚附近区域,这就是工程中常见的上下游坝脚附近隆起或滑动现象;最大竖向应力和水平应力均发生在坝底部区域。
图5 坝轴线处竖向位、竖向总应力和水平总应力随高度的变化
3.1 Et中各参数的影响
邓肯-张模型Et中有5个须经试验确定的参数:摩尔-库仑准则的黏聚力c、内摩擦角φ、试验参数k(相当于为100kPa时对应的Ei值)、初始模量Ei与围压力σ3曲线的斜率n以及破坏比Rf。下面以各参数的变化对坝顶A点的垂直位移SVA、坝脚B点的水平位移SHB和C点的总应力的影响进行敏感性分析。
3.1.1 黏聚力c
通过分析可知,黏聚力c对坝顶A点垂直位移SVA、坝脚B点水平位移SHB和C点竖向总应力的影响如图6所示。
坝顶A点垂直位移SVA和坝脚B点水平位移SHB均随c值的增加而减小,其变化率在c值减小时大于c值增大时,如c值减小50%时,SVA增大5.53%,SHB增大15.83%,c值增大50%时,SVA减少2.93%,SHB减少6.59%。
C点竖向总应力与c值增减成正比,其变化率明显小于位移的变化率(以下各参数均有类似情况),这是因为应力的大小主要由填筑土容重的大小决定,而与各参数变化的关系不大。同时,在计算过程中还发现,c值增减以后,堤坝的位移和应力等值线图与图2~图5所示的堤坝在基本试验值情况下得出的等值线图大致相同。
3.1.2 内摩擦角φ
内摩擦角φ的增减对坝顶A点垂直位移SVA、坝脚B点水平位移SHB和C点总应力的影响见图6~图8。
坝顶A点垂直位移SVA、坝脚B点水平位移SHB随φ值的增加而减小,φ值减小25%时,A点垂直位移SVA和B点水平位移分别增加3.95%和6.27%;φ值增加25%时,A点垂直位移SVA和B点水平位移SHB分别减小2.90%和4.77%。
C点竖向总应力与φ值增减成反比,其变化率明显小于位移的变化率。
3.1.3 参数n
初始模量Ei与围压力σ3成指数关系,n为指数,其影响见图6~图8。n值增减对位移和应力水平的影响较小。当n增加25%时,A点垂直位移SAB和B点水平位移SHB分别减小0.61%和增大3.29%;n值减小25%时,分别增大0.66%和减小3.08%。同样C点竖向总应力与n值增减成反比。因此,n值在整理分析时的局限性要小一些。
3.1.4 破坏比Rf
邓肯-张模型定义破坏比为:
Rf=(σ1-σ3)f/(σ1-σ3)u
式中:(σ1-σ3)f为土的抗剪强度取值;(σ1-σ3)u是土的极限抗剪强度理论值;多数取值Rf=0.70~0.90。Rf的影响见图6~图8。A点垂直位移SVA、B点水平位移SHB和C点总应力均随Rf的增加而增大,但变化率均不大。依据Rf的物理意义,Rf值越大,计算使用的抗剪强度(σ1-σ3)越接近土的理论极限抗剪强度(σ1-σ3)u。土体出现(σ1-σ3)f状态时,已经产生了相应的较大变形,濒于破坏。
3.1.5 参数k
k为σ3=100kPa时的初始模量Ei,但无量纲。文献[1]及大量试验资料表明,由于各种原因,k值可以相差数倍到几十倍。常规三轴试验过程中使用百分表量测轴向应变,试验点绘在普通坐标中求取k值。表面上看是εa→0时的纵坐标倒数,实际上多数情况下仅求取到εa=10-3级的相应值[1],与原本物理意义相去甚远,也是众多资料相差较大的原因之一。
本文算例中堤坝填筑土的k值是由常规三轴试验法求取的,k值对位移和应力的影响见图6~图8。
A点垂直位移SVA和B点水平位移SHB位移均与k值的变化率相反,当k变化较大时,位移与应力也相应地变化较大。
k从256减小至64,SVA增加119.05%,SHB增加204.88%,σV增加5.81%;k从256增至460.8时,SVA减小15.67%,SHB也减小35.83%,σV减小0.67%。
当k值变化率较大时,C点的竖向应力变化较大,其变化率大于其他参数引起的变化,例如当k值减少75%时,C点的竖向应力σV增加5.81%。
3.2 Kt中各参数的影响
3.2.1 参数kb
参数kb为σ3=100kPa时的体积模量Kt,但无量纲。与参数k一样,由于各种原因,kb值可以相差数倍到几十倍。kb值对位移和应力的影响见图6~图8。
A点垂直位移SVA的变化率与kb值的变化率相反,B点水平位移SHB的变化率与kb值的变化率相同,当kb变化较大时,位移也相应地变化较大。
当增加50%时,A点垂直位移SVA和B点水平位移SHB分别减小10.64%和增大13.73%;kb值减小50%时,分别增大32.89%和减小26.30%。
3.2.2 参数m
体积模量Kt与围压力σ3成指数关系,m为指数,其影响见图6~图8。由图可见,m值增减对位移和应力的影响不大;当m增加25%时,A点垂直位移SVA和B点水平位移SHB分别增大3.42%和减小2.11%;m值减小25%时,分别减小2.86%和增大1.85%。
3.3 各参数综合对比
将上述7个参数对A点垂直位移SVA、B点水平位移SHB和C点竖向应力σV的敏感性分析结果绘制成图,详见图6~图8。
图6 参数变化率对SVA的影响
图7 参数变化率对SHB的影响
图8 参数变化率对σV的影响
(1)A点垂直位移SVA随参数m、Rf增减而正向增减,随参数c、φ、k、n、kb变化成反向增减。
(2)B点水平位移SHB随参数kb、Rf和n增减而正向增减,随参数c、φ、k、m变化成反向增减。
(3)C点竖向应力σV随参数φ、n、Rf和kb增减呈反向增减,随参数c、k、m成正向增减,但是σV的变化幅度都非常小,这是因为应力的大小主要由填筑土容重大小决定,而与各参数变化关系不大。
(4)相比而言,参数k和kb对位移的影响最大,其变化率可达30%~40%,甚至80%以上,而且k和kb值本身的变化范围就很大,因此对结果的影响很重要。
参数c、φ和Rf对位移的影响次之,多数变化率在10%以内,虽然在本次的敏感性分析中,它们对位移和应力的影响不大,但是参数c和φ是土体的两个重要的抗剪强度指标,直接关系到土工结构的强度问题,而对于参数Rf,当它增大时,本身就意味着变形增大和相应的应力水平提高。
因此,以上5个参数对计算结果的敏感性都很高,在参数的确定过程中需要给予重视,一方面要求提高试验设备及量测仪器的精度;另一方面要求严格按原本物理力学意义进行试验,同时试验人员还要对计算结果作出合理的分析判断后,根据各种情况提出建议值,使计算结果与工程实际观测尽可能地吻合等。此外,参数n和m对位移的影响最小,因此,在整理分析时的局限性要小一些。
4 结论
E-K模型中的7个参数,除k值以外,由于试验条件和试验方法等原因,相差10%~40%是经常发生的事情。本文已证明,参数的上下浮动可引起位移5%~50%的变化率。至于k值,其物理意义是围压力σ3=100kPa的初始模量Ei[4]。在模型推导过程中明确定义Ei为εa→0时的模量,即要求在试验过程中尽可能地量测到微小应变,土体处于弹性状态时的模量。大多数土体应变εa从10-6增至10-2时,模量E将减小几十倍。多数试验由于量测设备精度不够,数据整理方法问题,实际上只取到εa=10-3数量级的E,与其物理力学意义有一定差距。仅就现行整理分析方法而言,也因经验方法的差异,所得结果经常相差几倍甚至几十倍,引起计算位移相差一倍之多。多数情况下表现为计算结果远大于实测值,主要是k值取值偏小。这就是模型计算结果与实测结果相差较大的主要原因之一。
本构模型中的这些问题,类似于土压力、土坡稳定分析、地基承载力等土力学的普遍问题,参数的影响或带来的误差可能比某些理论和方法要大得多。要克服上述不利因素,须提高仪器的量测精度,改变整理分析方法,对理论、试验和计算的全面深刻理解以及经验的积累。
参考文献
[1]钱家欢,殷宗泽.土工原理与计算[M].北京:中国水利水电出版社,1996:54-56.
[2]郑颖人,等.岩土塑性力学原理[M].北京:中国建筑工业出版社,2002:180-183.
[3]SIGMA/W for finite element stress and deformation analysis(User’s guide)[M].Geostudio,2004.
[4]李广信.关于Duncan双曲线模型参数确定的若干错误做法[J].岩土工程学报.1998,5(1):116.
本文发表于2011年。