邓肯-张非线性模型研究及其在ANSYS中的实现

文章编号:100722284(2010)0320076204邓肯2张非线性模型研究及其

在ANSYS 中的实现

宿　辉1,2,党承华2,崔佳佳2

(1.北京科技大学土木与环境学院,北京100083;2.河北工程大学水电学院,河北邯郸056021)

摘　要:对工程领域使用广泛的邓肯2张非线性本构模型进行了研究,总结了国内外的研究现状及理论成果,针对其无法判定因结压力降低时的加载情况,提出了相应的变形模量的计算方法,同时考虑中主应力、土体抗拉强度的影响对模型进行了修正。利用A PDL 编写程序实现了ANSYS 的材料本构模型的二次开发,运用重启动方法实现单元应力修正后数据重写入数据库,通过试验模拟对比分析验证了模型的适用性。

关键词:邓肯2张模型;抗拉强度;中主应力;应力分析

中图分类号:TU470+.3 文献标识码:A

Duncan 2Chang Nonlinear Elastic Model and R ealization in ANSYS

SU H ui 1,2,DANG Cheng 2hua 2,CUI Jia 2jia 2

(1.School of Civil and Environmen Engineering ,Beijing University of Science and Technology ,Beijing 100083,China ;

2.School of Water Resources and Hydropower ,Hebei University of Engineering ,Handan 056021,Hebei Province ,China )Abstract :Based on the present research situation and theoretical results ,research is done on Duncan 2Chang nonlinear elastic model ,which is applied widely in engineering ,.According to the fact that Duncan 2Chang model can ’t judge the loading situation while con 2solidation pressure decreases .The model is modified cosidering the effect of intermediate principal stress and tensile strength of soil.The secondary development of material constitutive model in ANSYS is accomplished by utilizing the A PDL language.Then the re 2start method is used to modify element stress ,the corrected data is database rewritten.A series of compliance tests verify the accura 2cy and applicability of the modified Duncan 2Chang nonlinear elastic model embedded ANSYS.

K ey w ords :Duncan 2Chang nonlinear elastic model ;tensile strength ;intermediate principal stress ;stress analysis

收稿日期:2009205218

作者简介:宿　辉(19722),男,副教授,博士研究生,主要从事岩土

工程及水工结构教学及研究工作。 岩土力应变关系十分复杂,变形主要是由颗粒间位置的变

化引起的,不同应力水平下相同的应力增量引起的应变增量并

不相同,从而表现出复杂的非线性特征[1]。为了反映土的特

性,人们提出了许多本构模型,以期更好地反映岩土的应力应

变规律,但直到目前为止还没有一种为人们所普遍接受的土体

本构模型。ANSYS 程序是我国科研机构及院校中常用的通用

有限元分析软件,具有完善的前后处理及求解功能,在土木结

构领域应用广泛。作为通用程序,ANSYS 缺少岩土常用的本

构模型,影响了其在土工分析中的应用,不仅使研究人员耗费

大量的精力编写或寻求其他软件,也造成了不必要的浪费,因

此对ANSYS 软件进行二次开发,利用其完善的前后处理功

能[2],建立适合土体的非线性本构模型是非常有必要的。1　邓肯2张模型及修正邓肯2张模型是非线性弹性模型,弹性矩阵中的弹性系数不再是常量,而是随应力状态而改变,因而可以反映土体应力-应变的非线性关系,并在一定程度上反映土变形的弹塑性。邓肯2张模型只有8个参数,参数的确定和计算分析较为简单方便,物理意义明确,只需用常规三轴压缩试验即可确定,因此在工程实践中得到了广泛应用,成为土的最为普及的本构模型之一[3-5]。该模型也存在一些缺陷,如无法反映土的剪胀性、不同应力路径的影响、加卸载判断不明确等,不可避免造成了计算分析误差,长期以来许多学者试图来对其进行修正,以更符合土体的响应[6]。在对邓肯2张模型进行改进和修正时,要认识到有些缺陷是可以通过对模型的修正来改善,而有些缺陷则是由于其理论基础的限制,无法通过对模型本身的改良来消

6

7中国农村水利水电?2010年第3期

除的,不能试图让本构完关系成其理论基础不支持的功能,如

反映土的剪胀性,应力路径影响等,否则会产生理论性的错误。1.1　E -v 模型和E -B 模型比较1970年邓肯等人在康纳(K ondner )提出的用双曲线拟合三轴试验曲线(σ3-σ1)～ε1的假定基础上,提出了双曲线E -v 模型[7],1980年邓肯等又提出了E -B 模型[8],引入体变模量B 代替切线泊松比v t 。2个模型中使用的切线泊松比E t 相同,不同的是E -v 模型使用切线泊松比v T ,E -B 模型使用的是体变模量B 。在E -v 模型中假定轴向应变ε1与侧向应变-ε3存在双曲线关系,从而得到-ε3/ε1与-ε3之间存在近似的直

线关系,从而可得到切线泊松比的表达式:

v t =G -F lg (σ3/p a )

1-D (σ1-σ3)KP 1(σ3p a )n [1-R f (σ1-σ3)(1-sin φ)2c cos φ+2σ3sin

φ]2=v i

(1-D (σ1-σ3)E I (1-R f S ))2

(1) 在E -B 模型中假定体应变εv 和轴向应变ε1为双曲线关系,而(σ1～σ3)～ε1关系也为双曲线,在一定的应力状态下(σ1-σ3)～εv 关系为直线,当σ3为常数时,B 与σ3在双对数坐标

中近似为直线关系,从而可得到体变模量的表达式:

B =K a p a (σ3

p a )m (2) 根据v ′t 和B 的关系可得到切线泊松比表达式(为清楚起

见,E -B 模型推导出的切线泊松比以v ′t 表示):

v ′t =0.5-(1-R f S )2K 6K b (σ3p a )n-m =0.5-E i (1-R f S )26K b p a (σ3

p a )m (3) 从以上分析可以看出,两模型的差异主要取决于切线泊松比的取值。研究表明2个模型的切线泊松比值的差值是非常显著,当σ3及S 水平不高时,E -v 模型的v t 要比E -B 模型v ′t 值为大,即有限元计算得出的变形较小,而σ3则比较大。对式(2)和式(3)进行研究可以得知,应力水平S 对切线泊松比的取值有一定的影响,2个模型中的切线泊松比都随应力水平S 的增大而增大,反映的规律是一致的。关于E -v 模型和E -B 模型哪个更适用的问题,一直存在在不同的意见。由于E -B

模型是在E -v 模型后提出,许多人认为前者更好一些,但E -

B 模型迄今为止只在文献[8]提出,并未在任何正式刊物或会议上发表,朱俊高在和DUNCAN 教授交流后认为,E -B 模型比E -v 模型更优的说法没有理论依据[9]。笔者认为对小粒径

材料如土的计算可采用E -v 模型,而对于在粒径材料如堆石

的计算采用E -B 模型可能更好一些。1.2　考虑中主应力影响的邓肯2张模型修正中主应力能够影响到土的抗剪强度,改变土体的应力-应变特性,对土体的变形有明显的影响[10]。邓肯2张模型采用的是摩尔-库仑破坏准则,这一准则是建立在破坏仅取决于最大和最小主应力这一假定的基础上,没有考虑中主应力σ2对强度和变形的影响,在π平面上的迹线与实验结果误差较大,从

而这将会导致土体变形模量偏低,计算出土体变形偏大[11],只有当中主应力和小主应力相差不大时,邓肯2张模型才能较好

地反映实际的应力、应变状态,两者相差越大,其计算误差也越

大[12]。在理论基础的框架内考虑中主应力的影响来对邓肯2张

模型进行修正,可使是邓肯2张模型更符合实际情况,增强模型

的适用性。对于中主应力的修正可以采用以下方法[13]:即用

来代替小主应力,用(σ2+σ3)/2来代替偏应力σ1-σ3;用球应

力p 和广义剪应力q 分别代替二维计算模型中的σ3和σ1-σ3;

以σ33σ2/σ3来代替σ3;以平面应变内摩擦角φp 代替轴对称应力条件下的内摩擦角φc 。中主应力对土体的强度和变形的影响机理非常复杂,中主应力及模型各参数之间会有相互影响,目前对此问题仍没有明晰的结论。1.3　考虑土体抗拉强度的邓肯2张模型修正土体是一种散粒体材料,相对于承载及抗剪能力而言,其抗拉性能一般不予以考虑。但随着高土石坝的建设,人们开始关注土体的裂缝问题以及土体的抗拉性能。为了在计算中考

虑土的抗拉强度,可以在计算初始切线变形模量时引入土体抗

拉强度:

E i =Kp a (σ3-σt p a )n (4)σt =-2c ×tan (45°

-φ2)(5)式中:σt 为土体的抗拉强度。

因此,切线变形模量的表达式如下:E t =E i (1-R f S )2=Kp a (

σ3-σt p a )n [1-R f (1-sin φ)(σ1-σ3)2c cos φ+2σ3sin φ]2(6)

2　邓肯2张模型卸载-再加载影响分析

2.1　卸载-再加载模量

在单调加载时,土的应力-应变曲线呈双曲线关系,而卸

载-再加载的性质则接近线弹性,为了反映土变形的可恢复和

不可恢复部分,邓肯2张模型在弹性理论框架内采用了不同于

初始加载模量E t 的卸载-再加载模量E ur 。此时考虑土的抗

拉强度,则E ur 可表示为:

E ur =K ur p a (σ3-σt p a )n (7)式中:K ur 为lg (E ur /p a )-lg (σ3/p a )直线的截距,n 为其斜率。

在卸载-再加载范围不大的情况下,可假定致密土至松散土的

K ur 取值范围为1.2K 至3.0K 。2.2　卸载-再加载判断准则邓肯2张模型在加卸载时使用了不同的变形模量,能够部

分反映土变形的塑性[14]。一般而言,土的卸载-再加载可分

为两种情况:偏差应力σ1-σ3的降低以及约束应力σ3的变化,

为综合考虑以上两种情况(尤其是σ3的变化),邓肯等于1984

年提出了加载函数的概念,但邓肯卸载-再加载准则中,遗漏

了土体固结压力降低的情况,造成了判断准则的矛盾。根据邓

肯的判断准则,当约束应力下降时,由S 1>S m 并不能保证(σ1

-σ3)1一定大于(σ1-σ3)m,此时变形模量也无法确定。为此,可以将邓肯2张模型的加-卸载准则补充如下:当S>S m情况下,若(σ1-σ3)1<(σ1-σ3)m,此时的变形模量可接下式计算:

E′t=(1-R f S)2E′i=

(1-R f (1-sinφ)(σ1-σ3)

2c cosφ+2σ3sinφ

)2Kp a(

σ3

m

p a

)n(8)

3　ANSYS中的二次开发

3.1　弹性常数动态变化在ANSYS中的实现

邓肯2张模型是弹性非线性模型,根据广义虎克定律来建立刚度矩阵,包含在矩阵中的弹性常数不再是常量,而是随应力状态相关的变量,即对于每一个计算单元,其弹性常数E、v 随单元的应力状态的变化而改变。当土体处于某一应力状态σ时,若施加微小的应力增量Δσ,则用此应力状态下的弹性常数形成刚度矩阵D,来计算相应的应变增量Δε。

在ANSYS中定义非线性材料的弹性常数一般使用数据表,即需事先输入相关的数据序列,计算过程中使用内插法求得材料的所需参数[15],2张模型中的弹性常数在计算过程中是随应力状态改变而改变的,与每一荷载步施加完成后的单元真实应力状态(σ1和σ3)有关,无法事先确定输入。为解决此问题,利用ANSYS内嵌的参量化设计语言A PDL进行的二次开发,编写Macro文件,通过逐级加载,在每一个载荷步完成后根据计算结果,修正弹性常数E t及v t,然后进行下一步载荷步的运算,从而实现弹性常数根据应力状态动态变化的要求。

3.2　基于APDL的邓肯2张模型

本文使用E-v模型,中主应力采用第3种修正方法,使用APDL来建立邓肯2张模型。由于邓肯2张模型中刚度矩阵不是常数,求解采用中点增量法,将荷载分为一系列的荷载增量{Δσ},用初始应力分量(即上一荷载步末的应力分量)计算弹性常数E′T和v′t,形成单元刚度矩阵组合成总刚度矩阵K′,求解出荷载步末的应力和应变,得到该级荷载的平均应力即中点应力(本级应力与初始应力取平均值),根据邓肯张模型求出其对应的E t和v t,然后求出弹性矩阵K,如此反复迭代即可求得最终结果。在计算过程中连续使用增量近似,不可避免会产生增量累积误差,可能导致最终结果不收敛,为避免此问题可使用N EWTON2RA P HSON平衡迭代来达到平衡收敛。

为完成上述过程,在每一个荷载步结束后,要多次退出求解器进入后处理器,提取计算值再次进行求解,并且需要返回前处理器修改参数,为了避免将前次计算结果丢失,必须要在保留计算结果的基础上进行下一荷载步的计算。

3.3　单元破坏后的修正

在增量计算过程中,随着荷载增量的不断增加,可能会在某些单元出现计算应力超过极限应力状态而产生破坏的情况,即拉裂破坏和剪切破坏。实际情况下应力是不可能超过极限破坏强度的,一旦出现上述情况必须对应力进行修正,修正过程中假定主应力方向不变。对于拉裂修正,如果计算σ3小于土的抗拉强度σt,即认为土体开裂不能再承受拉应力,修正σ3为0,同时假定σz不变,对其他应力分量进行修正;对于剪坏修正,当计算应力的摩尔圆超过库仑破坏线时,将摩尔圆修正到与破坏线相切,也修定σz不变,对其他应力分量进行修正。

在ANSYS运行中不能直接修改现有的数据库,必须采用重启动的方法[16],先进行初始分析,然后进入后处理器,对计算结果进行修正,然后采用修改了的结果等进行重启动分析。在重启动分析中要把多点重启动转换为单点重启动,在每次重启动前应将当前工作目录中的“.RST”、“.LD HI”和“. RNNN”文件删除,否则后续荷载步可能不能正常运行。

4　数值验证

为了检验所开发的邓肯-张本构模型的正确性,采用了室内三轴试验和数值模拟三轴试验对比的方法进行验证。试验采用邯郸粉质黏土,试验尺寸<39.1×80mm,进行了饱和固结排水三轴剪切试验,得到邓肯-张模型参数为:c=54kPa,φ= 27°,R f=0.84,K=262,n=0.35,G=0.366,F=0.184,D= 4.18,K ur=1014,n ur=0.41。根据轴对称理论建立试样1/4模型,映射划分单元,共产生2821个结点,2250个6面体单元,底面竖向(Z轴方向)约束,如图1所示。根据前节理论,笔者基于ANSYS的编制了修正邓肯2张非线性A PDL程序,进行三轴压缩仿真模拟试验,试验采用3种不同的围压条件(100, 200,400kPa)进行。表1记录了程序试样在400kPa围压下顶面圆心所在位置(45号节点)在各级荷载下的位移以及室内试验的实测值,通过比较可看出误差均在10%以内,说明二次开发的邓肯2张模型具有较高精度。为进一步对模型进行验证,重点对偏应力-轴向应变以及体积应变的关系进行了分析,并与室内试验结果进行了对比分析,见图2和图3。从图中可以看出,在低围压条件下,模拟结果和实验结果非常吻合,高围压下误差略有增大,但总体上结果也是相符的,表明在ANSYS 基础上进行的邓肯-张模型开发是完全可行的

。

图1　试样有限元模型

表1　400kP a围压下45号节点位移变化表

荷载/kPa50100150200300400500

位移/

mm

计算值0.5430.7910.962 1.174 1.720 2.498 3.643

实测值0.5720.843 1.001 1.287 1.88 2.591 3.771误差/% 5.34 6.57 4.059.639.30 3.72 3.51

图2　

不同围压下偏应力与轴向应变曲线图3　不同围压下体积应变与轴向应变曲线5　结　论本文对邓肯-张模型进行了综合分析,指出了E -v 模型

和E -B 模型的差异,在考虑土体抗拉强及中主应力影响的基础上对邓肯2张模型进行了修正。传统的邓肯卸载-再加载准则中,没有考虑到土体固结压力降低的情况,这是一个不足,在本次的修正中也对其进行了补充。在ANSYS 二次开发过程中,通过重启动解决了单元应力修正后写入数据库的难题,利用A PDL 开发的邓肯-张模型具有较好的准确性和稳定性,弥补了ANSYS 缺少岩土本构模型的缺陷,通过数值计算和试验结果的对比分析,可以确定修正邓肯-张模型较好地模拟了三轴试验,从而使ANSYS 这款软件强大的非线性求解平台能更好地应用于岩土的非线性特性分析,扩展了其用途。□参考文献:

[1]　李广信.高等土力学[M ].北京:清华大学出版社,2004.[2]　梁天锡,廖文冬,陈　强,等.基于二次开发进行猫弹性结构动力

学分析[J ].西安交通大学,力学与实践.2006,28(4):41-45.

[3]　周葆春,冯　波,汪　墨,等.武汉黏土Duncan 本构模型参数研究

[J ].岩土力学,2008,29(7):1788-1792.

[4]　何春保,沈建华,蔡　健.Duncan 2Chang 非线性地基模型的反演

分析[J ].岩土工程学报,2009,31(4):634-638.

[5]　牟声远;王正中.堆石料邓肯张模型的参数敏感性与统计分析

[J ].中国农村水利水电,2009,(3):97-100.

[6]　Kjell Karksrud and Lars Andresen.Loads on Braced Excavations

in Soft Clay [J ].International Journal of Geomechanics ,2005,5

(2):107-113.

[7]　Duncan J M ,Chang C Y.Nonlinear analysis of Stress and St rain

in Soils[J ].Journal of t he Soil Mechanics and Foundations Divi 2sion.1970,96(SM5):1629-1653.

[8]　Duncan J M ,Wong K S ,Mabry P.Strengt h ,stress 2strain and

bulk modulus parameters for finite element analysis of stress and movement s in soil masses [R ].University of California ,Berke 2ley ,1980.

[9]　朱俊高,周建方,邓　肯.E -v 模型和E -B 模型的比较[J ].水利

水电科技进展,2008,28(1):4-7.

[10]　Beiman A ,Nagy E J.Improvement of a Large Analytical Model

Using Test Data[J ].AIAA Journal ,1983,21(8):1168-1173.

[11]　Alvin K F.Met hod for Determining Minimum -Order Mass and

Stiffness Mat rices from Model Test Data [J ].AIAA Journal ,1995,33(1):128-135.

[12]　刘增荣,崔伟华,王　鑫.土体非线性本构模型的修正与土层特性分析[J ].西北工业大学学报,2008,26(4):454-457.[13]　毛德志,朱俊高.邓肯-张模型几种改进方法的比较[J ].岩土力学,2004,26(6):971-974.[14]　张　云,薛禹群,吴吉春,等.上海第四纪土层邓肯-张模型的参数研究[J ].水文地质工程地质,2008,(1):19-22.[15]　辛丽萍,王均星,陈　炜,等.三维饱和-非饱和瞬态渗流的ANSYS 模拟[J ].2009,(7):89-92.[16]　何益斌,夏栋舟,黄欣荣,等.土-桩-框架结构动力相互作用非线性有限元模拟[J ].沈阳建筑大学学报(自然科学版),2008,24(3):428-432.(上接第75页)　纵向最小;垂向雷诺切应力,随水深的增大呈递减趋势。

□参考文献:[1]　Von Karman T H ,Rubachh.On t he mechanism of resistance in

fluid[J ].Physikalische Zeit schrift ,1912,13(5):351-358.

[2]　N E 伊杰里奇克.水力摩阻手册[M ].北京:航空发动机编辑部,

1985:256-300.

[3]　Williamson C H K.Vortex dynamics in t he cylinder wake [J ].

Annual Review of Fluid Mechanics ,1996,28:477-539.

[4]　Williamson C H K ,G ovardhan R.Vortex 2in 2duced vibrations

[J ].Annual Review of Fluid Mechanics ,2004,36:413-455.

[5]　唐士芳,李　蓓.桩群阻力影响下的潮流数值模拟研究[J ].中国港湾建设,2001,(5):25-29.

[6]　张　玮,王　元.圆柱绕流涡系演变的DPIV 测试[J ].空气动力学报,2002,15(4):381-387.[7]　李　玲,李玉梁,祁雪春,等.浅水层圆柱绕流流动模式的实验研

究[J ].自然科学进展,2000,10(6):554-559.[8]　ШамовГИ.Речнынаносы,Гидрометеоиздат,Ленинград.1959.[9]　卢金友,徐海涛,姚仕明.天然河道水流紊动特性分析[J ].水利学报,2005,36(9):1029-1034.[10]　杨克君,刘兴年,曹叔尤,等.植被作用下的复式河槽漫滩水流紊动特性[J ].水利学报,2005,36(10):1263-1268.[11]　朱红钧,赵振兴.明渠河道紊流特性实验研究[J ].中国农村水利水电,2007,(9):12-15.[12]　卢金友,徐海涛,姚仕明.天然感潮河道水流紊动特性分析[J ].海洋工程,2005,23(3):70-77.