随动双线性硬化本构模型

1、2C =C清华大学研究生精品课程钢筋混凝土有限元教学程序C随动双线性硬化本构模型C Kin ematics hardi ng plastic subrout ine for MSC.MARCC陆新征江见鲸C LU Xin zhe ng , Jia ng Jianj ingC清华大学土木工程系，北京，100084C Dept. Civil En grg. of Tsin ghua Uni versityC last revised: Sep. 2003.C主要参考文献：C 1. MARC V olumn D user subroutine and Special RoutinesC 2. ABA

2、QUS Writ ing UMATs, VUMA Ts, and UELsC 3.江见鲸”钢筋混凝土结构非线性有限元分析 ”C =SUBROUTINE HYPELA(D,G,E,DE,S,TEMP,1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR)c implicit nonec* * * * * *cc user subroutine to define youngs modulus and poissons ratioc as function of stress in non-li near elastic small stra inc material.cc

3、d stress strain law to be formed by userc g cha nge in stress due to temperature effectsc e total strainc de in creme nt of strainc s stress - should be updated by userc temp state variablesc dtemp in creme nt of state variablesc ngens size of stress - stra in lawc n eleme nt nu mberc nn in tegrati

4、on point nu mberc kc layer nu mberc mats material i.d.c ndi nu mber of direct comp onentsc n shear nu mber of shear comp onentscc* * * * * *implicit real*8 (a-h,o-z)INCLUDE ./common/concom !通过concom模块得到当前的计算步数in teger : ngens,nn ,kc,mats ,n di, nshearreal*8 : e(1),de(1),temp(40),dtemp(40),g(1),d(nge

5、ns,ngens),s(1) ! temp(40)需要在前面initial condition里面设置完成in teger : n(2)C* * * * * *C Local ArraysCC EELAS - Elastic Strai ns !弹性应变C EPLAS - Plastic Strai ns !塑性应变C ALPHA - Shift Tensor !硬化参数C Flow - Plasitc Flow Directi ons !塑性流动方向C SIG - Stress at start of in creme nt !增量步开始时的应力C EPSPL - Plastic Strai

6、 ns at start of in creme nt !增量步开始时的塑性应变Creal*8 EELAS( nge ns),EPLAS( nge ns),ALPHA( nge ns)1 ,FLOW( nge ns),SIG( nge ns),EPSPL (n ge ns)c in teger : inc, in csub, n cyclereal*8 E0,ENU,EBULK,EG,ELAM初始弹性模量，泊松比，体积模量，剪切模量in teger : K1, K2 ! 循环变量real*8 : Smises, SYIELD, hard ! VM 应力，屈服应力，硬化模量real*8 : SI

7、GMQEQPL ! 平均正应力，等效应变增量real*8 : effg, efflam, effhard !等效剪切模量，等效lame常数，等效硬化模量Creal*8,parameter : ENUMAX=.49999D0,TOLER=1.D-6CC temp (1) - E ! 弹性模量C temp-NU ! 泊松比C temp-SYIELD !屈服强度C temp (4) - HARD !硬化模量 CC初始状态赋值Cif(in c=0.a nd.i ncsub=0.a nd.n cycle=0) the nTEMP(1)=200e3; !弹性模量 200GPaTEMP(2)=0.27; !

8、泊松比为 0.27TEMP(3)=210.0; ! 屈服强度 210MPaTEMP=20e3; !硬化模量 20GPaend ifE0=TEMP(1)ENU=Mi n( TEMP(2),ENUMAX)EBULK=E0/(1.-2.*ENU)/3.EG=EO心.+ENU)/2.ELAM=(EBULK*3.-EG*2.)/3.C弹性矩阵D=0.do K1=1, NDIdo K2=1, NDID(K2,K1)=ELAMend doD(K1,K1)=EG*2.+ELAMEND dodo K1=NDI+1, NGENSD(K1,K1)=EGend doC数组赋值EELAS=TEMP(4+1:4+NGEN

9、S) EPLAS=TEMP(4+NGENS+1:4+2*NGENS) ALPHA=TEMP(4+2*NGENS+1:4+3*NGENS)CC保存初始应力并计算试算应力Cdo k1=1, NGENSSIG(K1)=S(K1)EPSPL(K1)=EPLAS(K1)EELAS(K1)=EELAS(K1)+DE(K1)do k2=1,NGENS s(k2)=s(k2)+d(k2,k1)*de(k1) end doend doCC 计算 von mises 应力Csmises=(s(1)-alpha(1)-s (2)+alpha (2)*21 +(s (2)-alpha(2)-s (3)+alpha (

10、3)*22 +(s (3)-alpha(3)-s(1)+alpha(1)*2do K1=NDI+1,NGENSsmises=smises+6.*(s(k1)-alpha(k1)*2 end dosmises=sqrt(smises/2.)CC得到屈服应力和硬化模量CSYIELD=temp(3)hard=temp(4)CC判断是否屈服Cif(smises.gt.(1+toler)*syield) thenCC如果屈服，则计算硬化方向dF/dSIGCSIGM=(s(1)+s (2)+s (3) )/3.do K1=1, NDIflow(k1)=(s(k1)-alpha(k1)-SIGM)/smis

11、esend dodo k仁NDI+1, NGENSFlow(k1)=(s(k1)-alpha(k1)/smisesend doCC计算等效应变增量Cdeqpl=(smises-syield)/(EG*3.+HARD)CC Update shift tensor, elastic an plastic strains and stress Cdo k1=1, ndi alpha(k1)=alpha(k1)+hard*flow(k1)*deqpl eplas(k1)=eplas(k1)+3./2.*flow(k1)*deqpl eelas(k1)=eelas(k1)-3./2.*flow(k1)*

12、deqpl s(k1)=alpha(k1)+flow(k1)*syield+SIGM end dodo k1= ndi+1, ngens alpha(k1)=alpha(k1)+hard*flow(k1)*deqpl eplas(k1)=eplas(k1)+3.*flow(k1)*deqpl eelas(k1)=eelas(k1)-3.*flow(k1)*deqpl s(k1)=alpha(k1)+flow(k1)*syield end doCC得到切线刚度矩阵Ceffg=eg*(syield+hard*deqpl)/smisesefflam=(ebulk*3.-effg*2.)/3.effh

13、ard=eg*3.*hard/(eg*3.+hard)-effg*3.do k1=1, ndido k2=1, ndid(k2,k1)=efflamend dod(k1,k1)=effg*2.+efflamend dodo k1= ndi+1, ngensd(k1,k1)=effgend dodo k1=1, ngensdo k2=1, ngensd(k2,k1)=d(k2,k1)+effhard*flow(k2)*flow(k1) end doend doend ifCC保存当前弹性应变，塑性应变和硬化参数Cdo k1=1, ngenstemp(4+k1)=eelas(k1)temp(4+k1+ ngen s)=eplas(k1) temp(4+k1+2* ng