基于D_P准则的三维弹塑性有限元增量计算的有效算法

基于D-P准则的三维弹塑性有限元

增量计算的有效算法

A practical3D ela sto2pla stic incremental method in FEM

ba sed on D-P yield criteria

杨　强,陈　新,周维垣

(清华大学水利系,北京　100084)

摘　要:针对岩土材料常用的D-P准则,提出了一种新的增量分析方法,不用形成弹塑性增量矩阵,直接导出了符合正交流动法则的转移应力的解析解。该方法无论是对小步长还是大步长加载均有良好的收敛性。当采用精细的步长划分时,它就是严格意义上的理想弹塑性增量计算。在大步长情况下,在收敛域内最大载荷低于结构真实的极限承载力;对应的应力场是一个静力容许应力场;同时由于正交流动法则在平均意义下得到满足,收敛域内最大载荷接近结构真实的极限承载力。按此法所得结果接近真解且偏于安全。将整个计算模型装入三维非线性有限元程序TFI NE中,对某拱坝进行了超载分析。

关键词:转移应力;极限载荷;点安全度

中图分类号:T U452 文献标识码:A 文章编号:1000-4548(2002)01-0016-05

作者简介:杨　强(1964-),男,云南人。1988年在清华大学获硕士学位,1996年在奥地利Innsbruck大学获博士学位,现为清华大学水利系教授。主要从事水工结构及岩石力学方面的研究工作。

Y ANG Qiang,CHE N X in,ZH OU Wei2yuan

(Department of Hydraulic Engineering,Tsinghua University,Beijing100084,China)

Abstract:In this paper,focused on popularly used D-P yield criteria in geomaterials,a new incremental method in which the stresses to be trans2 ferred according to normal flow rule are directly derived without forming elasto2plastic increment matrix,was proposed.This method converges for either small load steps or large load steps.When very small load steps are used,the method is equivalent to standard elasto2plastic incremental method.When large load steps are used,the maximum load applied is lower than limit load in structure,the calculated stress field is an static ad2 missible one.As normal flow rule is satisfied in average,the maximum load is close to limit load.The soltion calculated by the method is on the safe side and close to real solution.The method was embedded into a3D nonlinear FEM software named TFINE,and overloading analysis was performed on an arch dam.

K ey words:the stresses to be transferred;normal flow rule;limit load

1　引 言Ξ

岩土材料具有很复杂的本构特性,如各向异性、硬化、软化等,目前描述岩土材料的本构模型非常多。但在实际工程三维有限元计算分析中,尤其是在岩体工程里,大量使用的仍是最简单D-P准则及理想弹塑性分析。其主要原因是参数选取不易。如在二滩高拱坝建设中,做了大量坝肩岩体现场大型抗剪试验,但具体到某一岩级,试验点数仍然很有限,且离散性很大。很难完全依赖试验确定参数,一般都要进行工程类比,对中、小工程工程类比更是参数确定的主要手段。最终一般只能给出岩体的抗剪参数f,c值。在这种情况下,从工程实用角度来说,追求本构关系的精致、完备并无太多实用意义。

相对而言,在岩土工程三维非线性有限元分析里,计算收敛性是一个较大的问题。弹塑性增量计算要采用精细的步长划分,才能确保计算收敛到正确解。在岩土工程,尤其是岩体工程里,荷载量级都很大,如高拱坝对水荷载的极限承载力可达上亿吨,而这对两岸高陡边坡所承受的的自重荷载来说,还只是一个小数,又如高地应力区大型地下洞室、高边坡(如三峡船闸高边坡)开挖过程中的释放荷载量级也十分巨大。若采用精细的步长划分,计算量将很大。岩体地质构造复杂,三维网格划分时经常会有畸形单元。由于地址缺陷或加固措施导致相邻单元材料性质差异过大,再加上高水平的荷载,各种因素交互影响,使得在计算过程中,经常出现局部发散现象,使得增量计算难以进行下去,最终结果可信度低,也难以从计算结果判断何时结构丧失稳定性。而对岩土工程来说,往往更关注结构的稳定性和极限承载力,而非应力和位移分布。

Ξ基金项目:国家自然科学基金资助项目(59879005);清华大学基础研究基金资助项目

收稿日期:2001-04-12

　第24卷　第1期岩　土　工　程　学　报V ol.24　N o.1 2002年 1月Chinese Journal of G eotechnical Engineering Jan.,　2002　

本文针对目前三维非线性有限元分析岩土材料常用的D -P 准则,提出了一种新的增量分析方法,它直接导出了符合正交流动法则的转移应力的解析解,不用形成弹塑性增量矩阵。该方法无论是对小步长还是大步长加载均有良好的收敛性。当采用精细的步长划分时,它就是严格意义上的理想弹塑性增量计算。当采用大步长增量计算,所得到是在平均意义上符合正交流动法则的一个静力容许应力场;在收敛域内,所能施加的最大荷载接近极限载荷的下限。

2　转移应力的解析解

在三维有限元弹塑性增量分析里,D -P 准则由于

简单实用,且和当前地质方面所能提供参数相适应,一直是目前应用最广的岩土材料屈服准则,其形式为

f =αI 1+

J 2-k ≤0

(1)

式中　J 2=

16

[(σ1-σ2)2+(σ2-σ3)2+(σ3-σ1)2],I 1=σ1+σ2+σ3,其中σ1,σ2,σ3为主应力。α和k 可通过拟合莫尔-库仑准则而得。例如,在π平面上,若

D -P 准则为库仑六边形的外接园,则

α=2sin φ3(3-sin φ),k =6c cos φ3(3-sin

φ)(2)式中　

φ和c 为材料的摩擦角和粘聚力。设某一高斯点在某一加载步或迭代步前初始应力为σ0且满足f (σ0)≤0。对某一个加载步或迭代步,由位移法求得该点应变增量为Δ

ε,它对应于弹性试应力σ1=σ1

ij =σ0+D :Δε,这里D 为弹性张量;若f (σ1)>0,则需进行应力调整。若此加载步或迭代步中塑性

应变增量为Δεp

,则调整后的应力为σ=σij =σ0+D :(Δε-Δεp

)=σ1-D :Δεp

。将正交流动法则d εp

=d λ

5f 5σ近似写成增量形式Δεp =Δ

λ5f 5σ,并以σ1确定5f Π5σ的代表值。则由条件

f (σ)=0,σ=σ1-ΔλD :5f 5σσ=σ1

(3)

即可确定调整后应力

σ=σij =(1-n )σ1

ij +p δij (4)

这里

n =

w

μJ 2

,p =-mw +

1

3nI 1,m =α(3λ+2μ

),w =f

3αm +μ

(5)

其中J 2,I 1,f 均由σ1确定。λ,μ为拉梅常数,即λ=E ν(1+ν)(1-2ν),μ=

E

2(1+ν

)(6)式中　E ,ν为杨氏模量和泊松比。容易证明由式(4)

确定的调整后应力σ必在屈服面上,即f (σ)=0。故对每个增量步或迭代步,应力转移值为

Δσ=σ1-σ=n σ1

ij -p δij (7) 求得由式(7)确定的转移应力Δσ后,可采用最普通的常刚度迭代进行应力转移计算。该方法避免了复杂的塑性矩阵的运算,直接求得满足正交流动法则的应力转移值。显然当采用精细的步长划分时,正交流动法则d εp

=d λ(5f Π5σ)得到充分满足,它就是严格意义上的理想弹塑性增量计算。以σ1确定5f Π5σ对计算收敛起到重要作用,此时转移应力的方向指向屈服锥体的中心线,所以调整后应力总能退回到屈服面上来。若以σ0确定5f Π5σ,如果σ1和σ0相差较大,则有可能调整后应力无法退回到屈服面上来。这实际上也是一般弹塑性计算当增量步过大就不收敛的重要原因之一。一般弹塑性计算也不进行校核。本文建议的方法在任何步长下都能确保调整后应力在屈服面上。大量的数值计算表明该方法具有良好的收敛性,所以可以用计算是否收敛来作为整体结构是否失稳的判据。

在大步长情况下求得的应力场与精细的步长划分下的所得是有差异的。我们可以从极限承载力的角度来探讨一下这个问题。在比例加载下,结构所能承受的最大荷载为极限承载力。在精细的步长划分下,理想弹塑性增量计算的最大荷载就是严格意义上的极限承载力。极限分析的下限定理认为:“一个满足平衡条件并且到处都不破坏材料屈服条件的应力场,是一个静力容许应力场;跟静力容许应力场对应的外荷载是极限荷载的下限;最高的下限便是极限荷载”。显然,在大步长情况下,在收敛域内最大载荷低于结构真实的极限承载力;对应的应力场是一个静力容许应力场;同时由于正交流动法则在平均意义下得到满足,收敛域内最大载荷接近结构真实的极限承载力。所以按此法所得的极限承载力接近真解且偏于安全。

3　其他考虑

岩土材料为低抗拉材料,故还应考虑抗拉条件:

σ1≤σt ,σ2≤σt ,σ3≤σt

(8)式中　σt 为材料单轴抗拉强度。在程序流程上,首先判断抗拉条件是否满足,若不满足,则调整应力满足之;其次再按D -P 准则判断调整应力。整个计算模型被装入三维非线性有限元程序TFI NE 里。在TFI NE 里,材料有三种破坏模式Ξ

:

7

1Ξ对理想弹塑性模型而言,屈服即为破坏。

　第1期杨　强,等1基于D -P 准则的三维弹塑性有限元增量计算的有效算法

(1)剪切破坏(相对于D -P 准则,式(1));(2)拉坏(式(8));

(3)拉-剪破坏,即材料既拉坏又剪坏。TFI NE 的塑性屈服区3种图案如图1,分别代表这

三种破坏模式。

图1　三种破坏模式图例

Fig.1　Patterns of three failure m odes

对弹性区域可由点安全系数P

来描述该点偏离

屈服面或安全程度。点安全系数的定义方式可以有很多。TFI NE 采用与D -P 准则相一致的点安全系数,即

P =(k -αI 1)Π

J 2

(9)

显然对任何应力状态Ξ

,必有P ≥1。进一步的细分包括:①弹性区,P >1;②

剪切破坏区,P =1;③拉坏

区,P ≥1;④拉-剪破坏区,P =1。

4　计算实例

用本文提出的方法对平面应变状态下的受内压厚壁圆筒进行弹塑性分析,圆筒内径a =100mm ,外径b =200mm ,厚度取为10mm 。利用对称性,取1Π4圆筒进行计算,有限元网格为10×20×1,见图2。取α=0,

k =σs /3,则D -P 准则退化为Mises 准则。材料参数

取为:弹模E 0=200G Pa ,泊松比μ=0.4,屈服强度σs =970MPa 。

图2　厚壁圆筒计算网格

Fig.2　Mesh of thick cylinder

平面应变状态下受内压厚壁圆筒服从Mises 准则

时,其弹塑性分析的理论解(参见文献[1]):弹性极限

为p e =

σs 31-a 2

b 2=420MPa ,塑性极限为p t =2σs

3ln

b

a

=776MPa 。塑性区由内向外扩展,弹塑性区分界面半径r p 与内压p 的关系为

p =

2σs 3

ln

r p a

+

σs 3

1-

r 2p b

2

(10)

根据以上计算条件,用本文提出的基于D -P 准则的新的增量分析方法,对平面应变状态下的受内压

厚壁圆筒进行了弹塑性分析。为了讨论大步长下的计算收敛性,在圆筒进入弹塑性状态后分别采用两种荷载步长ΔP 1=70MPa 和ΔP 2=35MPa 进行计算,发现两种荷载步长的计算结果基本相同,内压p =420MPa 时进入塑性状态,内压p =770MPa 时计算发散,为极限荷载。采用较大步长时,计算迭代5步即收敛,反映了该方法在采用大步长时有良好的收敛性。 圆筒在不同荷载步下的点安全度及塑性屈服区见图3。不同荷载步下弹塑性区分界面半径 r p 与理论解

r p 列于表1,可以看出误差很小,反映了该方法具有良

好的计算精度。

表1　弹塑性区分界面半径 r p 与理论解r p

T able 1　Calculated values r p and analytical s olution r p

内压p ΠMPa 有限元计算解

r p Πmm 理论解r p /mm 误差/%

490108.0109.0 1.0560116.3119.5 2.7630134.3132.8 1.1700149.0

150.00.7

图3　厚壁圆筒点安全度及屈服区

Fig.3　P oint safety factors and yield regions of thick cylinder

81Ξ一般多采用莫尔-库仑准则定义点安全系数,但这和D -P 准则

不配套,无法保证P ≥1。

岩　土　工　程　学　报　2002年　

图4　高程 360m 点安全度及屈服区

Fig.4　P oint safety factors and yield regions at elevation 360

m

图5　上游坝面点安全度及屈服区

Fig.5P oint safety factors and yield regions on upstream dam sur face

9

1　第1期杨　强,等1基于D -P 准则的三维弹塑性有限元增量计算的有效算法

图6　下游坝面点安全度及屈服区

Fig.6P oint safety factors and yield regions on downstream dam sur face

用本文提出的方法对某拱坝进行了三维有限元极

限承载力分析,基本荷载组合为自重+温度+渗流+水载。坝高为81m ,坝顶高程为 416m ,坝底高程为 335m 。节点总数6712,单元总数2641,实体单元总数2354,夹层单元总数287,坝体单元总数432。计算区域包括坝体和基础。共有12种材料分区,材料参数由设计部门给定,此处略。

拱坝在基本荷载组合及超1～4倍水载时高程 360m 平切面点安全度和屈服区图见图4,上、下游面点安全度及屈服区图分别见图5,6。在超载过程中,每一级加载步按1倍水载施加(大步长),在3倍超载下,收敛速度减慢,在4倍超载下,计算发散。故结构超载系数为4～5。和相应的地质力学模型试验对比,超载系数和破坏模式均很接近。

5　结 论

本文提出了一种新的三维有限元非线性增量分析

方法,直接求得满足正交流动法则的应力转移值而不用形成弹塑性增量矩阵,调整后应力总能退回到屈服面上来。当采用精细的步长划分时,它就是严格意义上的理想弹塑性增量计算。在大步长情况下,在收敛域内最大载荷低于结构真实的极限承载力;对应的应力场是一个静力容许应力场;同时由于正交流动法则在平均意义下得到满足,收敛域内最大载荷接近结构真实的极限承载力。所以按此法所得结果接近真解且偏于安全。由于该方法对小步长和大步长加载均有良好的收敛性,克服了一般弹塑性计算当增量步过大时不收敛的缺点,所以可以用计算是否收敛来作为整体结构是否失稳的判据,拱坝超载计算实例表明该法非常实用有效。参考文献:

[1]徐秉业,刘信声.应用弹塑性力学[M ].北京:清华大学出版

社,1995.

02岩　土　工　程　学　报　2002年