1、,主要内容,1. 通用有限元软件二次开发简介,2.用户自定义材料本构关系,3. 小结,通用有限元软件二次开发简介：ABAQUS,ABAQUS二次开发工具 用户子程序：Fortran,VC 自定义载荷，边界，本构关系，后处理 脚本语言：PYTHON 对ABAQUS功能进行全面的用户更新,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序 CREEP：定义时相关，黏塑性行为（蠕变和膨胀） DFLOW：定义固结分析中的非均匀空隙流速 DFLUX：定义热传导中的非均匀热流 DISP：定义非均布边界条件 DLOAD：定义非均布载荷 FILM：定义固结沉降分析中的非均匀渗流系数 FRIC：定
2、义接触面中的摩擦行为,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序 GAPCON：定义耦合温度-位移分析或纯热传导分析中接触面的导率 GAPELECTR：定义耦合热电分析中的导率 HARDINI：定义初始等效塑性应变和初始背应力张量 HETVAL：提供热传导分析中的内部热生成 MPC：定义多点约束 ORIENT：提供局部材料方向或运动耦合约束中的局部方向或惯性作用的局部刚体方向,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序 RSURFU：定义刚体面 SDVINI：定义初始求解相关的状态变量 SIGINI：定义初始应力场 UCORR：定义随机响应载荷的交叉
3、相关性质 UEL：定义单元 UEXPAN：定义热应变增量 UEXTERNALDB：管理用户定义的外部数据库并计算模型无关的历史信息,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序 UFIELD：定义预定义场变量 UFLUID：定义静水流体单元的流密度 UFLUIDLEAKOFF：定义空隙压力粘着单元的流体渗漏系数 UGENS：定义壳界面的力学行为 UHARD：定义各向同性或混合硬化的屈服面尺寸和硬化参数 UHYPEL：定义亚弹性应力应变关系,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序 UHYPER：定义超弹性材料 UINTER：定义接触面的面交互行为 U
4、MASFL：定义对流/扩散分析中的质量流动率条件 UMAT：定义材料的力学行为 UMATHT：定义材料热行为 UMESHMOTION：指定自适应网格中的网格运动约束 UMULLINS：定义Mullins效应材料模型的损伤变量,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序 UPOREP：定义初始流体空隙压力 UPRESS：指定等效压应力条件 UPSD：定义随机响应载荷的频率相关性 URDFIL：读结果文件 USDFLD：指定材料点的域变量 UTRACLOAD：指定非均布牵引载荷 UTRS：定义粘弹性材料的减缩时间平移函数 UVRAM：单元输出,通用有限元软件二次开发简介：AB
5、AQUS,ABAQUS用户子程序 UWAVE：定义ABAQUS/AQUA分析波运动 VOIDRI：定义初始空穴比,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例1：不同载荷步间改变弹性模量,Input file *HEADING 用户自定义损伤弹性模型（ USDFLD） *ELEMENT, TYPE=T2D2, ELSET=ONE 1, 1, 2 *NODE 1, 0., 0. 2, 10., 0. *SOLID SECTION, ELSET=ONE, MATERIAL=ELASTIC 1.,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例1,*MAT
6、ERIAL, NAME=ELASTIC *ELASTIC, DEPENDENCIES=1 2000., 0.3, 0., 0.00 1500., 0.3, 0., 0.01 1200., 0.3, 0., 0.02 1000., 0.3, 0., 0.04 *USER DEFINED FIELD *DEPVAR 1,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例1,*BOUNDARY 1, 1, 2 2, 2 *STEP *STATIC 0.1, 1.0, 0.0, 0.1 *CLOAD 2, 1, 20. *END STEP,通用有限元软件二次开发简介：ABAQUS,A
7、BAQUS用户子程序示例1,*STATIC 0.1, 1.0, 0.0, 0.1 *CLOAD 2, 1, 0. *END STEP *STEP, INC=20 *STATIC 0.1, 2.0, 0.0, 0.1 *CLOAD 2, 1, 40. *END STEP,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例1,SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT, TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,NDI,
8、NSHR,COORD,JMAC,JMATYP,MATLAYO, LACCFLA) INCLUDE ABA_PARAM.INC CHARACTER*80 CMNAME,ORNAME CHARACTER*3 FLGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3), 1 T(3,3),TIME(2) DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*), 1 COORD(*),通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例1,C Absolute value of c
9、urrent strain: CALL GETVRM(E,ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP, MATLAYO,LACCFLA) EPS = ABS( ARRAY(1) ) C Maximum value of strain up to this point in time: CALL GETVRM(SDV,ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP, MATLAYO,LACCFLA) EPSMAX = ARRAY(1) C Use the maximum strain as a field variable FIELD(1)
10、 = MAX( EPS , EPSMAX ) C Store the maximum strain as a solution dependent state C variable STATEV(1) = FIELD(1),通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例1,C If error, write comment to .DAT file: IF(JRCD.NE.0)THEN WRITE(6,*) REQUEST ERROR IN USDFLD FOR ELEMENT NUMBER , NOEL,INTEGRATION POINT NUMBER ,NPT END
11、IF C RETURN END,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例2：定义蠕变模型,通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例2：定义蠕变模型,SUBROUTINE CREEP(DECRA,DESWA,STATEV,SERD,EC,ESW,P,QTILD, 1 TEMP,DTEMP,PREDEF,DPRED,TIME,DTIME,CMNAME,LEXIMP,LEND, 2 COORDS,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) C INCLUDE ABA_PARAM.INC C CHARACTER*
12、80 CMNAME C DIMENSION DECRA(5),DESWA(5),STATEV(*),PREDEF(*),DPRED(*), 1 TIME(2),COORDS(*),EC(2),ESW(2),通用有限元软件二次开发简介：ABAQUS,ABAQUS用户子程序示例2：定义蠕变模型,A= SIG0= AN= T1=EXP(QTILD/SIG0) T2=EXP(QTILD/SIG0) DECRA(1) = A*(.5*(T1T2)*AN*DTIME IF(LEXIMP.EQ.1) THEN DECRA(5) = AN*A*(.5*(T1T2)*(AN1.)*DTIME/ 1 SIG0*.
13、5*(T1+T2) END IF RETURN END,通用有限元软件二次开发简介：ANSYS,ANSYS二次开发工具 APDL: 通过参数化模型来自动完成一些通用性强的任务； UIDL：用户界面设计语言，允许用户设计ANSYS图形界面； UPFs：Fortran90函数及例程以扩展或修改程序的功能，包括定义新材料本构，新单元，新的屈服准则，自定义优化算法，将ANSYS作为一个子程序来调用等。,通用有限元软件二次开发简介：ANSYS,ANSYS-APDL *Do循环,*do !起始行 !循环语句块 *enddo !结束行,不允许用label分支语句*if或*go命令跳出do循环语句； 不允许用
14、label将程序跳到另一行，但可以用ifthenelse来实现； do循环结构中，第一次循环后自动禁止命令结果输出；欲得到所有结果输出，在do循环结构中使用/gopr或/go语句； /clear命令不会清除do循环的堆栈，但它会删除所有的参数。可在/clear命令前运行/parsav命令来防止。,通用有限元软件二次开发简介：ANSYS,ANSYS-APDL,*do,i,0,64,1 !采用do循环设置移动载荷 *set,tim,tim+1 time,tim nsel,all fdele,all,all nsel,all nsel,s,loc,x,(64+i)*0.5 !施加5个集中载荷 f,a
15、ll,fy,-load nsel,all nsel,s,loc,x,(64+i-3)*0.5 *enddo,通用有限元软件二次开发简介：ANSYS,ANSYS-APDL,CYC1=2 !循环2圈 *DO,I,1, CYC1,1 !DO循环开始LSEL,S,LOC,Y,60TIME,5NSUBST,5AUTOTS,0SFL,ALL,PRES, PEAK1ALLSELLSWRITELSEL,S,LOC,Y,60AUTOTS,0NSUBST,5TIME,10SFL,ALL,PRES, 0ALLSEL,通用有限元软件二次开发简介：ANSYS,ANSYS-APDL,LSWRITELSEL,S,LOC,Y
16、,60TIME,12.5NSUBST,5AUTOTS,0SFL,ALL,PRES, VALLEY1ALLSELLSWRITELSEL,S,LOC,Y,60NSUBST,5AUTOTS,0TIME,15SFL,ALL,PRES, 0ALLSELLSWRITE *ENDDO !DO循环结束,通用有限元软件二次开发简介：ANSYS,ANSYS-APDL：Repeat,n,1,10,360/(nnode-1),0 !创建nnode个节点，夹角为360/nnode *repeat,nnode,1,0,360/nnode,0 !重复执行上述命令nnode次,通用有限元软件二次开发简介：ANSYS,ANSY
17、S-UIDL：单行参数输入 *ASK， Par，Query，DVAL 其中，Par为参数名称，用于存储用户输入的参数。Query是询问信息，用户可以输入最多包含54个字符串的提示信息以方便正确输入参数。DVAL是用户用空响应时程序自动赋给该参数的缺省值。用户用空格响应时则表示删除该参数。,*ASK,RADIUS,INPUT THE RADIUS OF CIRCLE,4,通用有限元软件二次开发简介：ANSYS,ANSYS-UIDL：多行参数输入 MULTIPRO,START,PROMPT_NUM *CSET,STRT_LOC,END_LOC,PARAM_NAME,PROMPT_STRING,DE
18、F_VALUE MULTIPRO,END 其中，START为第一个参数，用于标识MULTIPRO指令的开始； P ROMPT_NUM为一整型数，等于执行“MULTIPRO”命令行后的*CSET参数输入提示行的数目，至少有一个*CSET命令省略了DEF_VALUE参数或DEF_VALUE为0，才必须用到该参数。,通用有限元软件二次开发简介：ANSYS,ANSYS-UIDL：多行参数输入,MULTIPRO,START,3*CSET,1,3,EX_1,YOUNGS MODULUS(MPa),2.06E5 *CSET,4,6,NUXY_1,POISSIONS RATIO,0.3 *CSET,7,9,D
19、ENS_1,DENSITY(103Kg/mm3),7.85e-9 *CSET,61,62,ENTER THE ATTRIBUTES OF, MATERIAL 1 *CSET,63,64,NOTE:UNIT OF LENGHTH IS MM!, MULTIPRO,END /PREP7 MP,EX,1,EX_1 MP,NUXY,1,NUXY_1 MP,DENS1,1,DENS_1,通用有限元软件二次开发简介：ANSYS,ANSYS-UIDL：多行参数输入,通用有限元软件二次开发简介：ANSYS,ANSYS-UIDL：自定义工具栏,C：Program FilesAnsys Incv90ANSYSap
20、dlstart90.ans !视安装路径而定用文本编辑器打开后，在最后一行追加如下代码： !* /psearch,e:ansys !宏文件存放路径 *ABBR,Config,Kconfig *ABBR,etype,Ketype *ABBR,Mp,Kmp *ABBR,Biso,Kbiso *ABBR,Chaboche,kchaboche *ABBR,Model,Kmodel *ABBR,Load,Kload *ABBR,Solve,KSolve *ABBR,SS_center,Kscenter *ABBR,SS_edge,Kedge *ABBR,Avi_seqv,Kflash !*,通用有限元软
21、件二次开发简介：ANSYS,ANSYS-UIDL：多行参数输入,ANSYS高级工程应用实例与二次开发,通用有限元软件二次开发简介：MARC,MSC.MARC二次开发工具 用户子程序：Fortran 自定义载荷，边界，本构关系，后处理 脚本语言：PYTHON 用户开发新的软件功能，仅调用MARC的求解器。,通用有限元软件二次开发简介：MARC,subroutine ufxord(xord,ncrd,n)implicit real*8 (a-h,o-z) dp dimension xord(ncrd)r=24.d0/3.14159d0t=xord(2)/rxord(2)=r*sin(t)xord(
22、3)=r*cos(t) write (6,1) n,(xord(k),k=1,ncrd) 1 format(i5,3e13.5)returnend,通用有限元软件二次开发简介：MARC,二用户子程序定义本构模型,ABAQUS ：Plastic模型包括了各向同性硬化、随动硬化混合模型，Prager模型，AF模型等 UMAT ANSYS ：刚塑性、双线性、多段线性模型， 塑性Chaboche模型，ANADE模型,非 线性塑性等 USERMAT MARC ：刚塑性、双线性、多段线性模型，粘 塑性Chaboche模型等 HYPERLA2,本构模型简介,在材料力学行为的数学描述中，通常通过本构方程来表征
23、材料的反应，即在本构方程中给出应力作为物体变形历史的函数。这些本构方程即为本构模型。 不同的材料力学行为需要不同的本构模型来描述，如粘性流体、橡胶和混凝土，它们的本构模型截然不同。 在一维固体力学中，本构关系经常归属于材料的应力-应变关系。,本构模型简介,描述固体材料力学行为的本构模型大致归类：,弹性,非弹性,率无关弹性,率相关弹性,率无关非弹性,率相关非弹性,线弹性,非线性弹性,粘弹性,亚弹性,超弹性(伪弹性或拟塑性),双线性弹塑性,多线性弹塑性,非线性弹塑性,各向同性硬化,随动硬化,蠕变,粘塑性,各向同性硬化,随动硬化,固体材料唯象本构模型,仅蠕变,体积膨胀蠕变,非金属塑性,混凝土,岩土,
24、1,2,3,4,非各向同性硬化,非各向同性硬化,超塑性,线弹性本构模型的有限元实现,路径无关,可逆,无能量耗散,单轴状态下：,多轴状态下：,线弹性本构模型的有限元实现,其中，,有81个独立常数，与应力张量的9个分量对应。,应力和应变的对称性要求应力的6个独立分量与应变的6个独立分量有关。独立常数减少为36个。即66矩阵。 又由于矩阵的主对称性，独立弹性常数减少为n(n+1)/2=21个。,线弹性本构模型的有限元实现,由广义胡克定律可知：,线弹性本构模型的有限元实现,本构模型的离散：,线弹性本构模型的有限元实现,SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,S
25、PD,SCD,1 RPL,DDSDDT,DRPLDE,DRPLDT,2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) CINCLUDE ABA_PARAM.INC C,以ABAQUS为例，调用用户材料子程序UMAT,线弹性本构模型的有限元实现,CHARACTER*80 CMNAMEDIMENSION STRESS
26、(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),2 DDSDDT(NTENS),DRPLDE(NTENS),3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)DIMENSION DSTRES(6),D(3,3) C C PROPS(1) E 弹性模量 C PROPS(2) NU 泊松比 C,NTENS：总的应力或应变分量个数，三维问题为6个，平面和轴对称问题为4个，一维问题为
27、1个。 NTENS=NDI+NSHR,线弹性本构模型的有限元实现,C if (NDI.NE.3) THEN C WRITE (7, *) “THIS UMAT BE USED FOR 3D“ C CALL XIT C ENDIF C C ELASTIC PROPERTIESPARAMETER (ONE=1.0D0, TWO=2.0D0, THREE=3.0D0)EMOD=PROPS(1) EENU=PROPS(2) vEBULK3=EMOD/(ONE-TWO*ENU) K EG2=EMOD/(ONE+ENU) 2GEG=EG2/TWO GEG3=THREE*EG 3GELAM=(EBULK3-
28、EG2)/THREE ,C 计算弹性常数,线弹性本构模型的有限元实现,C C ELASTIC STIFFNESS CDO K1=1, NDIDO K2=1, NDIDDSDDE(K2, K1)=ELAMEND DODDSDDE(K1, K1)=EG2+ELAMEND DODO K1=NDI+1, NTENSDDSDDE(K1 ,K1)=EGEND DO C,计算弹性刚度矩阵C=DDSDDE(ntens,ntens),线弹性本构模型的有限元实现,C C CALCULATE STRESS CDO K1=1, NTENSDO K2=1, NTENSSTRESS(K2)=STRESS(K2)+DDSD
29、DE(K2, K1)*DSTRAN(K1)END DOEND DO CRETURNEND,计算当前应力,线弹性本构模型的有限元实现,1/8单元验证（C3D8）,*Heading * Job name: Linearelastic Model name: C3D8 *Preprint, echo=NO, model=NO, history=NO, contact=NO * PARTS *Part, name=PART-1-1 *Node1, 0., 0., 0.2, 1., 0., 0.3, 1., 1., 0.4, 0., 1., 0.5, 0., 0., 1.6, 1., 0., 1.7,
30、1., 1., 1.8, 0., 1., 1. *Element, type=C3D8 1, 1, 2, 3, 4, 5, 6, 7, 8 *Elset, elset=ALLE1, *Solid Section, elset=ALLE, material=ALLE 1., *End Part,线弹性本构模型的有限元实现,1/8单元验证（C3D8）,* ASSEMBLY * *Assembly, name=Assembly * *Instance, name=PART-1-1, part=PART-1-1 *End Instance * *Nset, nset=ALLN, instance=PA
31、RT-1-1, generate1, 8, 1 *Nset, nset=N1, internal, instance=PART-1-14, *Nset, nset=N2, internal, instance=PART-1-15, *Nset, nset=N3, internal, instance=PART-1-11, *Nset, nset=N4, internal, instance=PART-1-14, 5, 8 *Elset, elset=E1, internal, instance=PART-1-11, *Surface, type=ELEMENT, name=S1, intern
32、al E1, S4 *End Assembly,线弹性本构模型的有限元实现,1/8单元验证（C3D8）,*Amplitude, name=AMP-1 0., 0., 1., 1., 2., 0. * * MATERIALS * *Material, name=ALLE *Depvar0, *User Material, constants=2 2.e6, 0.3,线弹性本构模型的有限元实现,1/8单元验证（C3D8）,* BOUNDARY CONDITIONS * * Name: Disp-BC-1 Type: Displacement/Rotation *Boundary N1, 3, 3
33、* Name: Disp-BC-2 Type: Displacement/Rotation *Boundary N2, 2, 2 * Name: Disp-BC-3 Type: Displacement/Rotation *Boundary N3, 1, 1 * Name: Disp-BC-4 Type: Displacement/Rotation *Boundary N3, 2, 2 * Name: Disp-BC-5 Type: Displacement/Rotation *Boundary N3, 3, 3 * Name: Disp-BC-6 Type: Displacement/Rot
34、ation *Boundary N4, 1, 1,线弹性本构模型的有限元实现,1/8单元验证（C3D8）,* - * * STEP: Step-1 * *Step, name=Step-1, inc=200 *Static 0.1, 2., 0.1, 0.1 * * LOADS * * Name: SURFFORCE-1 Type: Pressure *Dsload, amplitude=AMP-1 S1, P, -200. * -,线弹性本构模型的有限元实现,1/8单元验证（C3D8）,* OUTPUT REQUESTS * *Restart, write, frequency=0 * *
35、FIELD OUTPUT: F-Output-1 * *Output, field *Element Output E, S * * HISTORY OUTPUT: H-Output-1 * *Output, history, variable=PRESELECT *End Step,弹塑性本构模型的有限元实现,首先以双线性弹塑性本构模型为例，阐述弹塑性有限元实现的过程。,应变率分解：,弹性本构方程：,屈服条件：,塑性流动法则：,加卸载条件：,一致性条件：,应力率-应变率关系：,连续体弹塑性切线模量：,关联流动法则，一般取,弹塑性本构模型的有限元实现,径向回退算法：,求当前应力：,弹塑性本构模
36、型的有限元实现,径向回退算法：,针对某一具体的屈服面，如：,弹塑性本构模型的有限元实现,连续体弹塑性切线模量：,弹塑性本构模型的有限元实现,径向回退算法1,1.设初始值,2.第K次迭代时检查屈服条件,收敛,4.更新塑性应变和内变量,5.计算连续体弹塑性切线模量,Else 转到第3步。,IF,3.计算塑性参数的增量,k=k+1，转到第2步,弹塑性本构模型的有限元实现,径向回退算法2,1.设初始值,2.检查屈服条件,弹性状态,4.更新塑性应变和内变量,5.计算连续体弹塑性切线模量,Else 转到第3步。,IF,3.计算塑性参数的增量,收敛,IF,Else 转到第3步。,ABAQUS,弹塑性本构模型
37、的有限元实现(各向同性硬化),SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS, 3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)C INCLUDE ABA_PARAM.INCC CHARACTER*80
38、MATERL DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3), 4 DFGRD0(3,3),DFGRD1(3,3),弹塑性本构模型的有限元实现(各向同性硬化),C DIMENSION EELAS(6),EPLAS(6),FLOW(6) PARAMETER (ONE=1.0D0,T
39、WO=2.0D0,THREE=3.0D0,SIX=6.0D0) DATA NEWTON,TOLER/10,1.D-6/ C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC PLASTICITYC J2 FLOW THEORY C CAN NOT BE USED FOR PLANE STRESS C C PROPS(1) E E C PROPS(2) NU V C PROPS(3) SYIELD y C PROPS(4) - HARDENING MODULUS H C-,弹塑性本构模型的有限元实现,C IF (NDI.NE.3) THEN WRITE(6,
40、1) 1 FORMAT(/,30X,*ERROR - THIS UMAT MAY ONLY BE USED FOR , ELEMENTS WITH THREE DIRECT STRESS COMPONENTS) ENDIF C C ELASTIC PROPERTIESC EMOD=PROPS(1) ENU=PROPS(2) SYIEL0=PROPS(3) HARD=PROPS(4) IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499 EBULK3=EMOD/(ONE-TWO*ENU) EG2=EMOD/(ONE+ENU) EG=EG2/TWO EG3=T
41、HREE*EG ELAM=(EBULK3-EG2)/THREE,弹塑性本构模型的有限元实现,C ELASTIC STIFFNESS C DO 20 K1=1,NTENS DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.0 10 CONTINUE CONTINUE C DO 40 K1=1,NDI DO 30 K2=1,NDI DDSDDE(K2,K1)=ELAM CONTINUE DDSDDE(K1,K1)=EG2+ELAM 40 CONTINUE DO 50 K1=NDI+1,NTENS DDSDDE(K1,K1)=EG 50 CONTINUE,弹塑性本构模型的有限元实现,
42、C CALCULATE STRESS FROM ELASTIC STRAINSC DO 70 K1=1,NTENS DO 60 K2=1,NTENS STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) CONTINUE 70 CONTINUE C RECOVER ELASTIC AND PLASTIC STRAINS DO 80 K1=1,NTENS EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS) CONTINUE EQPLAS=STATEV(1+2*NTENS),弹性应变,塑性应变
43、,等效塑性应变,C MISES STRESSC SMISES=(STRESS(1)-STRESS(2)*(STRESS(1)-STRESS(2) + (STRESS(2)-STRESS(3)*(STRESS(2)-STRESS(3) + (STRESS(3)-STRESS(1)*(STRESS(3)-STRESS(1) DO 90 K1=NDI+1,NTENS SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) CONTINUE SMISES=SQRT(SMISES/TWO)IF (SMISES.GT.(1.0+TOLER1)*SYIEL0) THEN C FLOW
44、 DIRECTIONC SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3)/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) CONTINUE DO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUE,弹塑性本构模型的有限元实现,C SOLVE FOR EQUIV STRESS, NEWTON ITERATION C SYIELD=SYIEL0 DEQPL=0.0 DO 130 KEWTON=1,NEWTON
45、RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD) EQPLAS= EQPLAS+DEQPL SYIELD=SYIEL0+HARD*DEQPL IF(ABS(RHS).LT.TOLER2*SYIEL0) GOTO 140 CONTINUE WRITE(6,2) NEWTON 2 FORMAT(/,30X,*WARNING - PLASTICITY ALGORITHM DID NOT , CONVERGE AFTER ,I3, ITERATIONS) 140 CONTINUE,弹塑性本构模型的有限元实现,C CALC STRESS AND
46、 UPDATE STRAINSC DO 150 K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO CONTINUE DO 160 K1=NDI+1,NTENS STRESS(K1)=FLOW(K1)*SYIELD EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL CON
47、TINUE EQPLAS=EQPLAS+DEQPL,弹塑性本构模型的有限元实现,C JACOBIANEFFHRD=EG3*HARD/(EG3+HARD) EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG EFFG3=THREE*EFFG2/TWO EFFLAM=(EBULK3-EFFG2)/THREE DO 220 K1=1,NDI DO 210 K2=1,NDI DDSDDE(K2,K1)=EFFLAM CONTINUE DDSDDE(K1,K1)=EFFG2+EFFLAM CONTINUE DO 230 K1=NDI+1,NTENS DDSDDE(K1,K1)=EFFG CONTINUE DO 250 K1=1,NTENS DO 240 K2=1,NTENS DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1) *(EFFHRD-EFFG3) CONTINUE CONTINUE ENDIF ENDIF,
