1、主要内容主要内容1.通用有限元软件二次开发简介通用有限元软件二次开发简介2.用户自定义材料本构关系用户自定义材料本构关系3.小结小结通用有限元软件二次开发简介：通用有限元软件二次开发简介：ABAQUSABAQUS二次开发工具二次开发工具用户子程序：用户子程序：Fortran,VC自定义载荷，边界，本构关系，后处理自定义载荷，边界，本构关系，后处理脚本语言：脚本语言：PYTHON对对ABAQUS功能进行全面的用户更新功能进行全面的用户更新通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序CREEP：定义时相关，黏塑性行为（蠕变和膨胀）：定义时相关，
2、黏塑性行为（蠕变和膨胀）DFLOW：定义固结分析中的非均匀空隙流速：定义固结分析中的非均匀空隙流速DFLUX：定义热传导中的非均匀热流：定义热传导中的非均匀热流DISP：定义非均布边界条件：定义非均布边界条件DLOAD：定义非均布载荷：定义非均布载荷FILM：定义固结沉降分析中的非均匀渗流系数：定义固结沉降分析中的非均匀渗流系数FRIC：定义接触面中的摩擦行为：定义接触面中的摩擦行为通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序GAPCON：定义耦合温度：定义耦合温度-位移分析或纯热传导分析位移分析或纯热传导分析 中接触面的导率中接触面的导率
3、GAPELECTR：定义耦合热电分析中的导率：定义耦合热电分析中的导率HARDINI：定义初始等效塑性应变和初始背应力张量：定义初始等效塑性应变和初始背应力张量HETVAL：提供热传导分析中的内部热生成：提供热传导分析中的内部热生成MPC：定义多点约束：定义多点约束ORIENT：提供局部材料方向或运动耦合约束中的局：提供局部材料方向或运动耦合约束中的局部方向或惯性作用的局部刚体方向部方向或惯性作用的局部刚体方向通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序RSURFU：定义刚体面：定义刚体面SDVINI：定义初始求解相关的状态变量：定义初始求
4、解相关的状态变量SIGINI：定义初始应力场：定义初始应力场UCORR：定义随机响应载荷的交叉相关性质：定义随机响应载荷的交叉相关性质UEL：定义单元：定义单元UEXPAN：定义热应变增量：定义热应变增量UEXTERNALDB：管理用户定义的外部数据库并计算：管理用户定义的外部数据库并计算模型无关的历史信息模型无关的历史信息通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序UFIELD：定义预定义场变量：定义预定义场变量UFLUID：定义静水流体单元的流密度：定义静水流体单元的流密度UFLUIDLEAKOFF：定义空隙压力粘着单元的流体渗：定义空
5、隙压力粘着单元的流体渗漏系数漏系数UGENS：定义壳界面的力学行为：定义壳界面的力学行为UHARD：定义各向同性或混合硬化的屈服面尺寸和硬：定义各向同性或混合硬化的屈服面尺寸和硬化参数化参数UHYPEL：定义亚弹性应力应变关系：定义亚弹性应力应变关系通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序UHYPER：定义超弹性材料：定义超弹性材料UINTER：定义接触面的面交互行为：定义接触面的面交互行为UMASFL：定义对流：定义对流/扩散分析中的质量流动率条件扩散分析中的质量流动率条件UMAT：定义材料的力学行为：定义材料的力学行为UMATHT：
6、定义材料热行为：定义材料热行为UMESHMOTION：指定自适应网格中的网格运动约束：指定自适应网格中的网格运动约束UMULLINS：定义：定义Mullins效应材料模型的损伤变量效应材料模型的损伤变量通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序UPOREP：定义初始流体空隙压力：定义初始流体空隙压力UPRESS：指定等效压应力条件：指定等效压应力条件UPSD：定义随机响应载荷的频率相关性：定义随机响应载荷的频率相关性URDFIL：读结果文件：读结果文件USDFLD：指定材料点的域变量：指定材料点的域变量UTRACLOAD：指定非均布牵引载
7、荷：指定非均布牵引载荷UTRS：定义粘弹性材料的减缩时间平移函数：定义粘弹性材料的减缩时间平移函数UVRAM：单元输出：单元输出通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序用户子程序UWAVE：定义：定义ABAQUS/AQUA分析波运动分析波运动VOIDRI：定义初始空穴比：定义初始空穴比通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例1：不同载荷步间改变弹性：不同载荷步间改变弹性模量模量Input file*HEADING用户自定义损伤弹性模型（用户自定义损伤弹性模型（USDFLD）*ELE
8、MENT,TYPE=T2D2,ELSET=ONE1,1,2*NODE1,0.,0.2,10.,0.*SOLID SECTION,ELSET=ONE,MATERIAL=ELASTIC1.通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例1*MATERIAL,NAME=ELASTIC*ELASTIC,DEPENDENCIES=12000.,0.3,0.,0.001500.,0.3,0.,0.011200.,0.3,0.,0.021000.,0.3,0.,0.04*USER DEFINED FIELD*DEPVAR1通用有限元软件二次开发简介通
9、用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例1*BOUNDARY1,1,22,2*STEP*STATIC0.1,1.0,0.0,0.1*CLOAD2,1,20.*END STEP通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例1*STATIC0.1,1.0,0.0,0.1*CLOAD2,1,0.*END STEP*STEP,INC=20*STATIC0.1,2.0,0.0,0.1*CLOAD2,1,40.*END STEP通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户
10、子程序示例用户子程序示例1SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT,TIME,DTIME,CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,NDI,NSHR,COORD,JMAC,JMATYP,MATLAYO,LACCFLA)INCLUDE ABA_PARAM.INCCHARACTER*80 CMNAME,ORNAMECHARACTER*3 FLGRAY(15)DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3),
11、1 T(3,3),TIME(2)DIMENSION ARRAY(15),JARRAY(15),JMAC(*),JMATYP(*),1 COORD(*)通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例1C Absolute value of current 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:CA
12、LL GETVRM(SDV,ARRAY,JARRAY,FLGRAY,JRCD,JMAC,JMATYP,MATLAYO,LACCFLA)EPSMAX=ARRAY(1)C Use the maximum strain as a field variableFIELD(1)=MAX(EPS,EPSMAX)C Store the maximum strain as a solution dependent stateC variableSTATEV(1)=FIELD(1)通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例1C If error,w
13、rite comment to.DAT file:IF(JRCD.NE.0)THENWRITE(6,*)REQUEST ERROR IN USDFLD FOR ELEMENT NUMBER,NOEL,INTEGRATION POINT NUMBER,NPTENDIFCRETURNEND通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例2：定义蠕变模型：定义蠕变模型通用有限元软件二次开发简介通用有限元软件二次开发简介：ABAQUSABAQUS用户子程序示例用户子程序示例2：定义蠕变模型：定义蠕变模型SUBROUTINE CREEP(DECR
14、A,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)CINCLUDE ABA_PARAM.INCCCHARACTER*80 CMNAMECDIMENSION DECRA(5),DESWA(5),STATEV(*),PREDEF(*),DPRED(*),1 TIME(2),COORDS(*),EC(2),ESW(2)通用有限元软件二次开发简介通用有限元软件二次开发简介：ABA
15、QUSABAQUS用户子程序示例用户子程序示例2：定义蠕变模型：定义蠕变模型A=SIG0=AN=T1=EXP(QTILD/SIG0)T2=EXP(QTILD/SIG0)DECRA(1)=A*(.5*(T1T2)*AN*DTIMEIF(LEXIMP.EQ.1)THENDECRA(5)=AN*A*(.5*(T1T2)*(AN1.)*DTIME/1 SIG0*.5*(T1+T2)END IFRETURNEND通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS二次开发工具二次开发工具APDL:通过参数化模型来自动完成一些通用性强通过参数化模型来自动完成一些通用性强的任务；的任务
16、；UIDL：用户界面设计语言，允许用户设计：用户界面设计语言，允许用户设计ANSYS图形界面；图形界面；UPFs：Fortran90函数及例程以扩展或修改程序函数及例程以扩展或修改程序的功能，包括定义新材料本构，新单元，新的屈服的功能，包括定义新材料本构，新单元，新的屈服准则，自定义优化算法，将准则，自定义优化算法，将ANSYS作为一个子程作为一个子程序来调用等。序来调用等。通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-APDL *Do循环循环*do !起始行起始行 !循环语句块循环语句块*enddo !结束行结束行不允许用label分支语句*if或*go命令跳出
17、do循环语句；不允许用label将程序跳到另一行，但可以用ifthenelse来实现；do循环结构中，第一次循环后自动禁止命令结果输出；欲得到所有结果输出，在do循环结构中使用/gopr或/go语句；/clear命令不会清除do循环的堆栈，但它会删除所有的参数。可在/clear命令前运行/parsav命令来防止。通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-APDL*do,i,0,64,1!采用do循环设置移动载荷*set,tim,tim+1 time,timnsel,allfdele,all,all nsel,allnsel,s,loc,x,(64+i)*0.5
18、!施加5个集中载荷f,all,fy,-loadnsel,allnsel,s,loc,x,(64+i-3)*0.5*enddo通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-APDLCYC1=2!循环2圈*DO,I,1,CYC1,1!DO循环开始 LSEL,S,LOC,Y,60 TIME,5 NSUBST,5 AUTOTS,0 SFL,ALL,PRES,PEAK1 ALLSEL LSWRITE LSEL,S,LOC,Y,60 AUTOTS,0 NSUBST,5 TIME,10 SFL,ALL,PRES,0 ALLSEL 通用有限元软件二次开发简介通用有限元软件二次开发
19、简介：ANSYSANSYS-APDLLSWRITE LSEL,S,LOC,Y,60 TIME,12.5 NSUBST,5 AUTOTS,0 SFL,ALL,PRES,VALLEY1 ALLSEL LSWRITE LSEL,S,LOC,Y,60 NSUBST,5 AUTOTS,0 TIME,15 SFL,ALL,PRES,0 ALLSEL LSWRITE*ENDDO!DO循环结束通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-APDL：Repeatn,1,10,360/(nnode-1),0!创建nnode个节点，夹角为360/nnode*repeat,nnode,1
20、,0,360/nnode,0!重复执行上述命令nnode次通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-UIDL：单行参数输入：单行参数输入*ASK，Par，Query，DVAL 其中，Par为参数名称，用于存储用户输入的参数。Query是询问信息，用户可以输入最多包含54个字符串的提示信息以方便正确输入参数。DVAL是用户用空响应时程序自动赋给该参数的缺省值。用户用空格响应时则表示删除该参数。*ASK,RADIUS,INPUT THE RADIUS OF CIRCLE,4通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-UIDL：多行参
21、数输入：多行参数输入MULTIPRO,START,PROMPT_NUM*CSET,STRT_LOC,END_LOC,PARAM_NAME,PROMPT_STRING,DEF_VALUEMULTIPRO,END其中，START为第一个参数，用于标识MULTIPRO指令的开始；P ROMPT_NUM为一整型数，等于执行“MULTIPRO”命令行后的*CSET参数输入提示行的数目，至少有一个*CSET命令省略了DEF_VALUE参数或DEF_VALUE为0，才必须用到该参数。通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-UIDL：多行参数输入：多行参数输入MULTIPR
22、O,START,3*CSET,1,3,EX_1,YOUNGS MODULUS(MPa),2.06E5*CSET,4,6,NUXY_1,POISSIONS RATIO,0.3*CSET,7,9,DENS_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/PREP7MP,EX,1,EX_1MP,NUXY,1,NUXY_1MP,DENS1,1,DENS_1通用有限元软件二次开发简介通用有限元软件二次
23、开发简介：ANSYSANSYS-UIDL：多行参数输入：多行参数输入通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-UIDL：自定义工具栏：自定义工具栏C：Program FilesAnsys Incv90ANSYSapdlstart90.ans!视安装路径而定用文本编辑器打开后，在最后一行追加如下代码：!*/psearch,e:ansys !宏文件存放路径*ABBR,Config,Kconfig*ABBR,etype,Ketype*ABBR,Mp,Kmp*ABBR,Biso,Kbiso*ABBR,Chaboche,kchaboche*ABBR,Model,Kmod
24、el*ABBR,Load,Kload*ABBR,Solve,KSolve*ABBR,SS_center,Kscenter*ABBR,SS_edge,Kedge*ABBR,Avi_seqv,Kflash!*通用有限元软件二次开发简介通用有限元软件二次开发简介：ANSYSANSYS-UIDL：多行参数输入：多行参数输入ANSYS高级工程应用实例与二次开发高级工程应用实例与二次开发通用有限元软件二次开发简介通用有限元软件二次开发简介：MARCMSC.MARC二次开发工具二次开发工具用户子程序：用户子程序：Fortran自定义载荷，边界，本构关系，后处理自定义载荷，边界，本构关系，后处理脚本语言：脚本
25、语言：PYTHON用户开发新的软件功能，仅调用用户开发新的软件功能，仅调用MARC的求解器。的求解器。通用有限元软件二次开发简介通用有限元软件二次开发简介：MARC subroutine ufxord(xord,ncrd,n)implicit real*8(a-h,o-z)dpdimension xord(ncrd)r=24.d0/3.14159d0 t=xord(2)/r xord(2)=r*sin(t)xord(3)=r*cos(t)write(6,1)n,(xord(k),k=1,ncrd)1 format(i5,3e13.5)return end通用有限元软件二次开发简介通用有限元软件
26、二次开发简介：MARC二用户子程序定义本构模型二用户子程序定义本构模型ABAQUS：Plastic模型包括了各向同性硬化、随模型包括了各向同性硬化、随 动硬化混合模型，动硬化混合模型，Prager模型，模型，AF模模 型等型等 UMATANSYS ：刚塑性、双线性、多段线性模型，：刚塑性、双线性、多段线性模型，塑性塑性Chaboche模型，模型，ANADE模型模型,非非 线性塑性等线性塑性等 USERMATMARC ：刚塑性、双线性、多段线性模型，粘：刚塑性、双线性、多段线性模型，粘 塑性塑性Chaboche模型等模型等 HYPERLA2 本构模型简介本构模型简介在材料力学行为的数学描述中，通
27、常通过在材料力学行为的数学描述中，通常通过本构方程来表征材料的反应，即在本构方本构方程来表征材料的反应，即在本构方程中给出应力作为物体变形历史的函数。程中给出应力作为物体变形历史的函数。这些本构方程即为本构模型。这些本构方程即为本构模型。不同的材料力学行为需要不同的本构模型不同的材料力学行为需要不同的本构模型来描述，如粘性流体、橡胶和混凝土，它来描述，如粘性流体、橡胶和混凝土，它们的本构模型截然不同。们的本构模型截然不同。在一维固体力学中，本构关系经常归属于在一维固体力学中，本构关系经常归属于材料的应力材料的应力-应变关系。应变关系。本构模型简介本构模型简介描述固体材料力学行为的本构模型大致归
28、类：描述固体材料力学行为的本构模型大致归类：弹性弹性非弹性非弹性率无关弹性率无关弹性率相关弹性率相关弹性率无关非弹性率无关非弹性率相关非弹性率相关非弹性线弹性线弹性非线性弹性非线性弹性粘弹性粘弹性亚弹性亚弹性超弹性超弹性(伪弹性或拟塑性伪弹性或拟塑性)双线性弹塑性双线性弹塑性多线性弹塑性多线性弹塑性非线性弹塑性非线性弹塑性各向同性硬化各向同性硬化随动硬化随动硬化蠕变蠕变粘塑性粘塑性各向同性硬化各向同性硬化随动硬化随动硬化固体固体材料材料唯象唯象本构本构模型模型仅蠕变仅蠕变体积膨胀蠕变体积膨胀蠕变非金属塑性非金属塑性混凝土混凝土岩土岩土1234非各向同性硬化非各向同性硬化非各向同性硬化非各向同性
29、硬化超塑性超塑性 线弹性本构模型的有限元实现线弹性本构模型的有限元实现路径无关路径无关可逆可逆无能量耗散无能量耗散单轴状态下：单轴状态下：E多轴状态下：多轴状态下：C:线弹性本构模型的有限元实现线弹性本构模型的有限元实现C:C其中，其中，有有81个独立常数，与应力张量的个独立常数，与应力张量的9个分量对应。个分量对应。应力和应变的对称性要求应力的应力和应变的对称性要求应力的6个独立分量与应变的个独立分量与应变的6个独立分量有关。独立常数减少为个独立分量有关。独立常数减少为36个。即个。即66矩阵。矩阵。又由于矩阵的主对称性，独立弹性常数减少为又由于矩阵的主对称性，独立弹性常数减少为n(n+1)
30、/2=21个。个。121323332211665655464544363534332625242322161514131211121323332211222CCCCCCCCCCCCCCCCCCCCC对称 线弹性本构模型的有限元实现线弹性本构模型的有限元实现由广义胡克定律可知：由广义胡克定律可知：IIICG2)21)(1(vvvEe)1(2vEG)21(332vEGKGGGGGGeeeeeeeee222C线弹性本构模型的有限元实现线弹性本构模型的有限元实现本构模型的离散：本构模型的离散：11nnnC:1nnC:1n 线弹性本构模型的有限元实现线弹性本构模型的有限元实现 SUBROUTINE UM
31、AT(STRESS,STATEV,DDSDDE,SSE,SPD,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)C INCLUDE ABA_PARAM.INCC 以以ABAQUS为例，调用用户材料子程序为例，调用用户材料子程序UMAT 线弹性本构模型的
32、有限元实现线弹性本构模型的有限元实现CHARACTER*80 CMNAME DIMENSION STRESS(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(
33、2)NU 泊松比泊松比CNTENS：总的应力或应变分量个数，三维问题为：总的应力或应变分量个数，三维问题为6个，平面和轴对称问题为个，平面和轴对称问题为4个，个，一维问题为一维问题为1个。个。NTENS=NDI+NSHR线弹性本构模型的有限元实现线弹性本构模型的有限元实现Cif(NDI.NE.3)THEN CWRITE(7,*)THIS UMAT BE USED FOR 3DCCALL XITCENDIFCC ELASTIC PROPERTIESPARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0)EMOD=PROPS(1)EENU=PROPS(2)vEBULK3
34、=EMOD/(ONE-TWO*ENU)K EG2=EMOD/(ONE+ENU)2GEG=EG2/TWO GEG3=THREE*EG 3GELAM=(EBULK3-EG2)/THREE C 计算弹性常数计算弹性常数 线弹性本构模型的有限元实现线弹性本构模型的有限元实现CC ELASTIC STIFFNESSCDO K1=1,NDIDO K2=1,NDIDDSDDE(K2,K1)=ELAMEND DODDSDDE(K1,K1)=EG2+ELAMEND DODO K1=NDI+1,NTENSDDSDDE(K1,K1)=EGEND DOC计算弹性刚度矩阵计算弹性刚度矩阵C=DDSDDE(ntens,n
35、tens)GGGGGGeeeeeeeee222C线弹性本构模型的有限元实现线弹性本构模型的有限元实现CC CALCULATE STRESSCDO K1=1,NTENSDO K2=1,NTENSSTRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)END DOEND DOCRETURNEND 计算当前应力计算当前应力C:1nn线弹性本构模型的有限元实现线弹性本构模型的有限元实现1/8单元验证（单元验证（C3D8）*Heading*Job name:Linearelastic Model name:C3D8*Preprint,echo=NO,model=NO,
36、history=NO,contact=NO*PARTS*Part,name=PART-1-1*Node 1,0.,0.,0.2,1.,0.,0.3,1.,1.,0.4,0.,1.,0.5,0.,0.,1.6,1.,0.,1.7,1.,1.,1.8,0.,1.,1.*Element,type=C3D81,1,2,3,4,5,6,7,8*Elset,elset=ALLE 1,*Solid Section,elset=ALLE,material=ALLE1.,*End Part 线弹性本构模型的有限元实现线弹性本构模型的有限元实现1/8单元验证（单元验证（C3D8）*ASSEMBLY*Assembl
37、y,name=Assembly*Instance,name=PART-1-1,part=PART-1-1*End Instance*Nset,nset=ALLN,instance=PART-1-1,generate 1,8,1*Nset,nset=N1,internal,instance=PART-1-1 4,*Nset,nset=N2,internal,instance=PART-1-1 5,*Nset,nset=N3,internal,instance=PART-1-1 1,*Nset,nset=N4,internal,instance=PART-1-1 4,5,8*Elset,elset
38、=E1,internal,instance=PART-1-1 1,*Surface,type=ELEMENT,name=S1,internalE1,S4*End Assembly 线弹性本构模型的有限元实现线弹性本构模型的有限元实现1/8单元验证（单元验证（C3D8）*Amplitude,name=AMP-10.,0.,1.,1.,2.,0.*MATERIALS*Material,name=ALLE*Depvar 0,*User Material,constants=22.e6,0.3 线弹性本构模型的有限元实现线弹性本构模型的有限元实现1/8单元验证（单元验证（C3D8）*BOUNDARY
39、CONDITIONS*Name:Disp-BC-1 Type:Displacement/Rotation*BoundaryN1,3,3*Name:Disp-BC-2 Type:Displacement/Rotation*BoundaryN2,2,2*Name:Disp-BC-3 Type:Displacement/Rotation*BoundaryN3,1,1*Name:Disp-BC-4 Type:Displacement/Rotation*BoundaryN3,2,2*Name:Disp-BC-5 Type:Displacement/Rotation*BoundaryN3,3,3*Name
40、:Disp-BC-6 Type:Displacement/Rotation*BoundaryN4,1,1线弹性本构模型的有限元实现线弹性本构模型的有限元实现1/8单元验证（单元验证（C3D8）*-*STEP:Step-1*Step,name=Step-1,inc=200*Static0.1,2.,0.1,0.1*LOADS*Name:SURFFORCE-1 Type:Pressure*Dsload,amplitude=AMP-1S1,P,-200.*-线弹性本构模型的有限元实现线弹性本构模型的有限元实现1/8单元验证（单元验证（C3D8）*OUTPUT REQUESTS*Restart,wri
41、te,frequency=0*FIELD OUTPUT:F-Output-1*Output,field*Element OutputE,S*HISTORY OUTPUT:H-Output-1*Output,history,variable=PRESELECT*End Step 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现首先以双线性弹塑性本构模型为例，阐述弹塑性有限元实现的过程。首先以双线性弹塑性本构模型为例，阐述弹塑性有限元实现的过程。应变率分解：应变率分解：pe弹性本构方程：弹性本构方程：)(:peCC屈服条件：屈服条件：),(qrp塑性流动法则：塑性流动法则：hq 0),(qf加
42、卸载条件：加卸载条件：0,0,0ff一致性条件：一致性条件：0:fffrCrhCr:.:qf应力率应力率-应变率关系：应变率关系：C:eprCrhCrrCCC:.):():(qepf连续体弹塑性切线模量：连续体弹塑性切线模量：)(rfp 关联流动法则，一般取关联流动法则，一般取1h 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现径向回退算法：径向回退算法：求当前应力：求当前应力：)(:111pnnnC)(:11pnpnnC)(:1pnpnnnCpnnpnnCC1:)(:CpnnnCC1:):(pnnC1*1:rC:1*1nn01nfn*1n1nGHfk3)(:)(10*pynpykkfr
43、C 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现径向回退算法：径向回退算法：rCss:1*11devnnn)(3)(10*pynpykkGf针对某一具体的屈服面，如：针对某一具体的屈服面，如：rsrCss10*10*2:ndevnkGn)23232(10*nG10*3nkGns32ssn nr5.1GHfk3rC:1*11nnn)(pykkf 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现rCrhCrrCCC:.):():(qepf连续体弹塑性切线模量：连续体弹塑性切线模量：rrCG2:G3:rCrGHGGdeep3422rrIIIC 弹塑性本构模型的有限元实现弹塑性本构模型的有
44、限元实现径向回退算法径向回退算法11.设初始值设初始值2.第K次迭代时检查屈服条件检查屈服条件收敛4.更新塑性应变和内变量更新塑性应变和内变量5.计算连续体弹塑性切线模量计算连续体弹塑性切线模量0k00)(:1*pnnC)(3*pkykkGf1TOLfkElse 转到第3步。IFGHfkk33.计算塑性参数的增量计算塑性参数的增量)()()1(kpkpkpn23)(kkp)()1(kpkkk=k+1，转到第2步GHGGdeep3422rrIIIC 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现径向回退算法径向回退算法21.设初始值设初始值2.检查屈服条件检查屈服条件弹性状态4.更新塑性应
45、变和内变量更新塑性应变和内变量5.计算连续体弹塑性切线模量计算连续体弹塑性切线模量0k00C:*)(*pyf1TOLf Else 转到第3步。IF3.计算塑性参数的增量计算塑性参数的增量)()()1(kpkpkpn23)(kkp)()1(kpkkGHGGdeep3/122nnIIICGHGpkykk3)(3*2)(3*TOLGypkykk收敛收敛IFElse 转到第3步。ABAQUS 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现(各向同性硬化各向同性硬化)SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1 RPL,DDSDDT,DRPL
46、DE,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 MATERL DIMENSION STRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRP
47、LDE(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,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)DATA NEWTON,TOLER/10,1.D-6/C UMAT FOR ISOTR
48、OPIC ELASTICITY AND ISOTROPIC PLASTICITYC J2 FLOW THEORYC CAN NOT BE USED FOR PLANE STRESSCC PROPS(1)E EC PROPS(2)NU VC PROPS(3)SYIELD yC PROPS(4)-HARDENING MODULUS HC-弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现C IF(NDI.NE.3)THEN WRITE(6,1)1 FORMAT(/,30X,*ERROR-THIS UMAT MAY ONLY BE USED FOR,ELEMENTS WITH THREE DIRE
49、CT STRESS COMPONENTS)ENDIFCC 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=THREE*EG 1ELAM=(EBULK3-EG2)/THREE 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现C ELASTIC STIFFNESSC DO 20 K1=1,NTE
50、NS DO 10 K2=1,NTENS DDSDDE(K2,K1)=0.0 10 CONTINUE CONTINUEC 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 2050 CONTINUE 弹塑性本构模型的有限元实现弹塑性本构模型的有限元实现C CALCULATE STRESS FROM ELASTIC STRAINSC DO 70 K1=1,NTENS DO 60 K2=1,
