ABAQUS-材料本构模型及编程

PNEWDT

作为新的时间增量

计算；这对时间 相关的材料如聚合物等有用；如果设为大余

1

的数，则下一个增 量步加大

DTIME

为

DTIME

X

PNEWDT

。可以更新。

其他变量含义可参看手册，暂时用不到。

C

INCLUDE

'ABA_'

定义了一些参数，变量什么的，不用管

C

CHARACTER*8

CMNAME

C

DIMENSION

STRESS( NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),

1

DDS DDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS ),

2

PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3, 3),

3

DFGRD0(3,3),DFGRD1(3,3)

矩阵的尺寸声明

C

C

LOCAL

ARRAYS

C

--------------- -------------------------------------------------

C

EELAS

-

ELASTIC

STRAINS

C

EPLAS

-

PLASTIC

STRAINS

C

FLOW

-

DIRECTION

OF

PLASTIC

FLOW

C

--------------- -------------------------------------------------

C

局部变量，用来暂时保存弹性应变、塑性应变分量以及流动方向

DIMENSION

EELAS(6),EPLAS(6),FLOW(6)

C

PARAMETER(ZERO=0.D0,ONE=1.D 0,TWO=2.D0,THREE=3.D0,SIX=6.D0,

1

ENUMAX=.4999D0,NEWTON=10,TOLER=1.0D-6)

C

C

-- -------------------------------------------------- ------------

C

UMAT

FOR

ISOTROPIC

ELASTICITY

AND

ISOTROPIC

MISES

PLASTICITY

C

CANNOT

BE

USED

FOR

PLANE

STRESS

C

--------------- -------------------------------------------------

C

PROPS(1)

-

E

C

PROPS(2)

-

NU

C

PROPS(3..)

-

SYIELD

AN

HARDENING

DATA

C

CALLS

HARDSUB

FOR

CURVE

OF

YIELD

STRESS

VS.

PLASTIC

STRAIN

C

--------------- -------------------------------------------------

精品资料

____________ __________________________________________________ ________________________________________________

C

C

ELASTIC

PROPERTIES

C

获取杨氏模量，泊松比，作为已知量由

PROPS

向量传入

EMOD=PROPS(1)

E

ENU=PROPS(2)

ν

EBULK3=EMOD/(ONE-TWO*ENU)

3K

EG2=EMOD/(ONE ENU)

2G

EG=EG2/TWO

G

EG3=THREE*EG

3G

ELAM=(EBULK3-EG2)/THREE

λ

DO

K1=1,NTENS

DO

K2=1,NTENS

DDSDDE(K1,K2)=ZERO

END

DO

END

DO

弹性部分，

Jacobian< /p>

矩阵很容易计算

注意，在

ABAQUS

中，剪切应变采用工程剪切应变的定义，所以剪切部分模量是< /p>

G

而不是

2G!

C

C

ELASTIC

STIFFNESS

C

DO

K1=1,NDI

DO

K2=1,NDI

DDSDDE(K2,K1)=ELAM

END

DO

DDSDDE(K1,K1)=EG2 ELAM

END

DO

DO

K1=NDI 1,NTENS

DDSDDE(K1,K1)=EG

END

DO

C

C

RECOVER

ELASTIC

AND

PLASTIC

STRAINS

AND

ROTATE

FORWARD

C

ALSO

RECOVER

EQUIVALENT

PLASTIC

STRAIN

C

读取弹性应变分量，塑性应变分量，并旋转（调用了

ROTSIG

），分别保存在

EELAS

和

EPLAS

中；

CALL

ROTSIG(STATEV(

1),DROT,EELAS,2,NDI,NSHR)

CALL

ROTSIG(STATEV(NTENS 1),DROT,EPLAS,2,NDI,NSHR)

读取等效塑性应变

EQPLAS=STATEV(1 2*NTENS)

精品资料

____________ __________________________________________________ ________________________________________________

先假设没有发生塑性流动，按完全弹性变形计算试算应力

C

C

CALCULATE

PREDICTOR

STRESS

AND

ELASTIC

STRAIN

C

DO

K1=1,NTENS

DO

K2=1,NTENS

STRESS(K2)=STRESS(K2) DDSDDE(K2,K1)*DSTRAN(K1)

END

DO

EELAS(K1)=EELAS(K1) DSTRAN(K1)

END

DO

C

计算

Mises

应力

C

CALCULATE

EQUIVALENT

VON

MISES

STRESS

C

SMISES=(STRESS(1)-STRESS(2))**2 (STRESS(2)-STRESS(3))**2

1

(STRESS(3)-STRESS(1))**2

DO

K1=NDI 1,NTENS

SMISES=SMISES SIX*STRESS(K1)**2

END

DO

SMISES=SQRT(SMISES/TWO)

C

根据当前等效塑性应变，调用

HARDSUB

得到对应的屈服应力

C

GET

YIELD

STRESS

FROM

THE

SPECIFIED

HARDENING

CURVE

C

NVALUE=NPROPS/2-1

CALL

HARDSUB(SYIE L0,HARD,EQPLAS,PROPS(3),NVALUE)

C

C

DETERMINE

IF

ACTIVELY

YIELDING

C

如果

Mi ses

应力大余屈服应力，屈服发生，计算流动方向

IF

(.(ONE TOLER)*SYIEL0)

THEN

C

C

ACTIVELY

YIELDING

C

SEPARATE

THE

HYDROSTATIC

FROM

THE

DEVIATORIC

STRESS

C

CALCULATE

THE

FLOW

DIRECTION

C

SHYDRO=(STRESS(1) STRESS(2) STRESS(3))/THREE

DO

K1=1,NDI

FLOW(K1)=(STRESS(K1)-SHYDRO)/SMISES

END

DO

DO

K1=NDI 1,NTENS

FLOW(K1)=STRESS(K1)/SMISES

精品资料