随动双线性硬化本构模型.doc

C =C 清华大学研究生精品课程钢筋混凝土有限元教学程序C 随动双线性硬化本构模型 C Kinematics harding plastic subroutine for MSC.MARCC 陆新征 江见鲸C LU Xinzheng, Jiang JianjingC 清华大学土木工程系，北京，100084C Dept. Civil Engrg. of Tsinghua UniversityC last revised: Sep. 2003.C 主要参考文献：C 1. MARC Volumn D user subroutine and Special RoutinesC 2. ABAQUS Writing UMATs, VUMATs, 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-linear elastic small strainc material.cc d stress strain law to be formed by userc g change in stress due to temperature effectsc e total strainc de increment of strainc s stress - should be updated by userc temp state variablesc dtemp increment of state variablesc ngens size of stress - strain lawc n element numberc nn integration point numberc kc layer numberc mats material i.d.c ndi number of direct componentsc nshear number of shear componentscc* * * * * *implicit real*8 (a-h,o-z)INCLUDE ./common/concom ! 通过concom模块得到当前的计算步数integer : ngens,nn,kc,mats,ndi,nshearreal*8 : e(1),de(1),temp(40),dtemp(40),g(1),d(ngens,ngens),s(1) ! temp(40)需要在前面initial condition 里面设置完成integer : n(2)C* * * * * *C Local ArraysC-C EELAS - Elastic Strains ! 弹性应变C EPLAS - Plastic Strains ! 塑性应变C ALPHA - Shift Tensor ! 硬化参数C Flow - Plasitc Flow Directions ! 塑性流动方向C SIG - Stress at start of increment ! 增量步开始时的应力C EPSPL - Plastic Strains at start of increment ! 增量步开始时的塑性应变Creal*8 EELAS(ngens),EPLAS(ngens),ALPHA(ngens)1 ,FLOW(ngens),SIG(ngens),EPSPL(ngens)c integer : inc, incsub, ncyclereal*8 E0,ENU,EBULK,EG,ELAM ! 初始弹性模量，泊松比，体积模量，剪切模量integer : K1, K2 ! 循环变量real*8 : Smises, SYIELD, hard ! VM应力, 屈服应力，硬化模量real*8 : SIGM,DEQPL ! 平均正应力, 等效应变增量real*8 : effg, efflam, effhard ! 等效剪切模量，等效lame常数，等效硬化模量Creal*8,parameter : ENUMAX=.49999D0,TOLER=1.D-6CC temp (1) - E ! 弹性模量C temp (2) - NU ! 泊松比C temp (3) - SYIELD ! 屈服强度C temp (4) - HARD ! 硬化模量 C -C 初始状态赋值Cif(inc=0.and.incsub=0.and.ncycle=0) thenTEMP(1)=200e3; ! 弹性模量 200GPaTEMP(2)=0.27; ! 泊松比为 0.27TEMP(3)=210.0; ! 屈服强度 210MPaTEMP(4)=20e3; ! 硬化模量 20GPaend ifE0=TEMP(1)ENU=Min(TEMP(2),ENUMAX)EBULK=E0/(1.-2.*ENU)/3.EG=E0/(1.+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+NGENS)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,NGENSs(k2)=s(k2)+d(k2,k1)*de(k1)end doend doC C 计算von mises 应力Csmises=(s(1)-alpha(1)-s(2)+alpha(2)*21 +(s(2)-alpha(2)-s(3)+alpha(3)*22 +(s(3)-alpha(3)-s(1)+alpha(1)*2do K1=NDI+1,NGENSsmises=smises+6.*(s(k1)-alpha(k1)*2end dosmises=sqrt(smises/2.)CC 得到屈服应力和硬化模量CSYIELD=temp(3)hard=temp(4)CC 判断是否屈服Cif(smises.gt.(1+toler)*syield) then CC 如果屈服，则计算硬化方向 dF/dSIGCSIGM=(s(1)+s(2)+s(3)/3.do K1=1, NDIflow(k1)=(s(k1)-alpha(k1)-SIGM)/smisesend dodo k1=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 stressCdo k1=1, ndialpha(k1)=alpha(k1)+hard*flow(k1)*deqpleplas(k1)=eplas(k1)+3./2.*flow(k1)*deqpleelas(k1)=eelas(k1)-3./2.*flow(k1)*deqpls(k1)=alpha(k1)+flow(k1)*syield+SIGMend dodo k1=ndi+1, ngensalpha(k1)=alpha(k1)+hard*flow(k1)*deqpleplas(k1)=eplas(k1)+3.*flow(k1)*deqpleelas(k1)=eelas(k1)-3.*flow(k1)*deqpls(k1)=alpha(k1)+flow(k1)*syieldend doCC 得到切线刚度矩阵Ceffg=eg*(syield+hard*deqpl)/smisesefflam=(ebulk*3.-effg*2.)/3.effhard=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=