25.3.4 VUMAT
User subroutine to define material behavior.定义材料本构⽤户⼦程序Product: ABAQUS/Explicit
Warning: The use of this user subroutine generally requires considerable expertise. You are cautioned that theimplementation of any realistic constitutive model requires extensive
development and testing. Initial testing on a single-element model with prescribed traction loading is strongly recommended.注意：⽤户⼦程序的使⽤通常需要⼀定的专长。⽤户需要知道执⾏任何实际的本构模型需要⼤量的试验数据。强烈建议⽤户对⽤户⼦程序进⾏在指定拉⼒作⽤下单个单元的验证测试。
The component ordering of the symmetric and nonsymmetric tensors for the three-dimensional case using C3D8R elementsis different from the ordering specified in “Three-dimensional solid element library,” Section 14.1.4, and the ordering used inABAQUS/Standard.
C3D8R单元三维轴对称及⾮轴对称张量成分顺序与“Three-dimensional solid elementlibrary,” Section 14.1.4及ABAQUS/Standard中指定的顺序不同。References
“User-defined mechanical material behavior,” Section 12.8.1*USER MATERIALOverview
User subroutine VUMAT:⽤户⼦程序VUMAT
is used to define the mechanical constitutive behavior of a material;⽤来定义材料的⼒学本构关系
will be called for blocks of material calculation points for which the material is defined in a user subroutine (“Material datadefinition,” Section 9.1.2);
可以被⽤户⼦程序定义材料计算点调⽤
can use and update solution-dependent state variables;可以使⽤和更新结果依赖状态变量
can use any field variables that are passed in;可以使⽤传⼊的任何场变量
is described further in “User-def ined mechanical material behavior,” Section 12.8.1; and 在“User-defined mechanicalmaterial behavior,” Section 12.8.1中详细论述；cannot be used in an adiabatic analysis.可以被⽤于绝热分析
Component ordering in tensors 张量组成顺序
The component ordering depends upon whether the tensor is symmetric or nonsymmetric.张量组成顺序取决于其是否为对称或⾮对称张量。
Symmetric tensors 对称张量
For symmetric tensors such as the stress and strain tensors, there are ndir+nshr components, and the component order isgiven as a natural permutation of the indices of the tensor. The direct components are first and then the indirect components,beginning with the 12-component. For example, a stress tensor contains ndir direct stress components and nshr shear stresscomponents, which are passed in as
对于如同应⼒及应变张量等的对称张量，含有ndir+nshr分量，分量的序号按张量索引号的⾃然排序给出。⾸先是直接分量，然后是从12分量开始的间接分量。例如，包含ndir正应⼒分
engineering components; this is different from user subroutine UMAT in ABAQUS/Standard, which uses engineeringcomponents.
Nonsymmetric tensors ⾮对称张量
For nonsymmetric tensors there are ndir+2*nshr components, and the component order is given as a natural permutation ofthe indices of the tensor. The direct components are first and then the indirect components, beginning with the 12-component. For example, the deformation gradient is passed as
对于⾮对称张量含有ndir+2*nshr分量，分量的顺序按照张量索引号的⾃然排序给出。⾸先是
Initial calculations and checks 最初计算和检查
In the data check phase of the analysis ABAQUS/Explicit calls user subroutine VUMAT with a set of fictitious strains and atotalTime and stepTime both equal to 0.0. This is done as a check on your constitutive relation and to calculate theequivalent initial material properties, based upon which the initial elastic wave speeds are computed.
在ABAQUS/Explicit调⽤⽤户⼦程序VUMAT分析的数据检查阶段，⼩应变、总时间及时间步都为0。这作为对⽤户本构关系的⼀个检查，基于计算得到的初始材料波速来计算等效初始材料属性。Defining local orientations 定义局部⽅向
All stresses, strains, stretches, and state variables are in the orientation of the local material axes. These local material axesform a basis system in which stress and strain components are stored. This represents a corotational coordinate system inwhich the basis system rotates with the material. If a user-specified coordinate system (“Orientations,” Section 2.2.5) is used,it defines the local material axes in the undeformed configuration.
所有的应⼒、应变、延伸及状态变量均按局部材料轴的⽅向。这些局部材料轴形成⼀个应⼒与应变分量存储的基本系统。即这
个基本系是随着材料联合转动的坐标系。如果使⽤⽤户指定坐标系，则它在⽆变形结构中定义局部材料轴。Special considerations for various element types不同单元类型的特殊考虑
The use of user subroutine VUMAT requires special consideration for various element types.⽤户⼦程序VUMAT的使⽤需要对不同的单元类型进⾏特殊的考虑。Shell and plane stress elements壳及平⾯应⼒单元
You must define the stresses and internal state variables. In the case of shell or plane stress elements, you must definestrainInc(*,3), the thickness strain increment. The internal energies can be defined if desired. If they are not defined, theenergy balance provided by ABAQUS/Explicit will not be meaningful.
⽤户必须定义应⼒和初始状态变量。在壳或平⾯应⼒单元的情况下，⽤户必须定义应变包括（*,3），厚度应变增量。如果需要的话还需要定义初始能量。如果没有定义，那么ABAQUS/Explicit提供的能量平衡将没有意义。Shell elements壳单元
When VUMAT is used to define the material response of shell elements, ABAQUS/Explicit cannot calculate a default valuefor the transverse shear stiffness of the element. Hence, you must define the element's transverse shear stiffness. See “Shellsection behavior,” Section 15.6.4, for guidelines on choosing this stiffness.
当使⽤VUMAT定义壳单元的材料响应时，ABAQUS/Explicit不能计算单元的缺省横向剪切刚度。因此，⽤户需要定义单元的横向剪切刚度。关于选择横⾏剪切刚度的详细资料请参考“Shell section behavior,” Section 15.6.4Beam elements梁单元
For beam elements the stretch tensor and the deformation gradient tensor are not available. For beams in space you mustdefine the thickness strains, strainInc(*,2) and strainInc(*,3). strainInc(*,4)
is the shear strain associated with twist. Thickness stresses, stressNew(*,2) and stressNew(*,3), are assumed to be zero andany values you assign are ignored.
对于梁单元不能使⽤拉伸张量及位移梯度张量。对于空间梁，⽤户必须定义厚度应变、应变增量（*，2）及应变增量
（*，3）。应变增量（*，4）时与扭曲有关的剪应变。厚度应⼒，（*，2）及（*，3）假定为0，并且⽤户分配的任何相关张量都被忽略。
Deformation gradient 位移梯度
The polar decomposition of the deformation gradient is written as , where
and are the right and left symmetric stretch tensors, respectively. The constitutive model is
defined in a corotational coordinate system in which the basis system rotates with the material. All stress and strain tensorquantities are defined with respect to the corotational basis system. The right stretch tensor, , is used. The relative spin tensorrepresents the spin (the antisymmetric
part of the velocity gradient) defined with respect to the corotational basis system.
位移梯度写成，其中及分别为右边及左边的对称拉伸张量。本构模型定义为联合旋转坐标系，在该坐标系中基系随着材料转动。所有的应⼒和应变张量值按照联合旋转坐标系定义。使⽤右边的拉伸位移。相应的旋转张量代表与联合旋转基系相应的转动。
Special considerations for hyperelasticity超弹性的特殊考虑
Hyperelastic constitutive models in VUMAT should be defined in a corotational coordinate system in which the basis systemrotates with the material. This is most effectively accomplished by formulating the hyperelastic constitutive model in terms ofthe stretch tensor, , instead of in terms
of the deformation gradient, . Using the deformation gradient can present some
difficulties because the deformation gradient includes the rotation tensor and the resulting stresses would need to be rotatedback to the corotational basis.
在VUMAT中的超弹性本构模型可以被定义在联合选择坐标系中。这可以通过⽤拉伸张量
表⽰的超弹性本构模型很好的实现，⽽不是使⽤位移剃度来表⽰。使⽤位移剃度可能会带来⼀些困难，因为位移剃度包括旋转张量并且导致应⼒需要选择返回到联合旋转基系。Objective stress rates⽬标应⼒率
The Green-Naghdi stress rate is used when the mechanical behavior of the material is defined using user subroutine
VUMAT. The stress rate obtained with user subroutine VUMAT may differ from that obtained with a built-in ABAQUS materialmodel. For example, most material models used with solid (continuum) elements in ABAQUS/Explicit employ the Jaumannstress rate. This difference in the formulation will cause significant differences in the results only if finite rotation of a materialpoint is accompanied by finite shear. For a discussion of the objective stress rates used in ABAQUS, see “Stress rates,”Section 1.5.3 of the ABAQUS Theory Manual.
在⽤户⼦程序VUMAT中使⽤Green-Naghdi应⼒率来定义材料的⼒学本构关系。通过⽤户⼦程序VUMAT获得的应⼒率可能会与在ABAQUS建⽴的材料模型获得的应⼒率有所不同。例如，在ABAQUS/Explicit中⼤多数实体（连续）单元材料模型使⽤Jaumann应⼒率。只要材料点的有限旋转伴随有限剪切，这种表达⽅式的不同将导致计算结果的明显差异。关于ABAQUS中使⽤的⽬标应⼒率的详细讨论参考“Stress rates,” Section 1.5.3 of the ABAQUS Theory Manual.Material point deletion 材料点删除
Material points that satisfy a user-defined failure criterion can be deleted from the model (see “User-defined mechanicalmaterial behavior,” Section 12.8.1). You must specify the state variable number controlling the element deletion flag whenyou allocate space for the solution-dependent state variables, as explained in “User-defined mechanical material behavior,”Section 12.8.1. The deletion state variable should be set to a value of one or zero in VUMAT. A value of one indicates thatthe material point is active, while a value of zero indicates that ABAQUS/Explicit should delete the material point from themodel by setting the stresses to zero. The structure of the block of material points passed to user subroutine VUMAT remainsunchanged during the analysis; deleted material points are not removed from the block. ABAQUS/Explicit will pass zero
stresses and strain increments for all deleted material points. Once a material point has been flagged as deleted, it cannot bereactivated.
满⾜⽤户定义的破坏准则的材料点可以被从模型中删除（参考“User-defined mechanical material behavior,” Section 12.8.1)。当⽤户给结果依赖状态变量分配空间时，⽤户需要指定控制单元删除标⽰的状态变量号，在“User-defined mechanicalmaterial behavior,” Section 12.8.1
中进⾏详细说明。在VUMAT中删除状态变量可以被赋予1或者0。1表⽰材料点时激活的，0表⽰ABAQUS/Explicit将通过设定应⼒为0删除材料点。在分析过程中传递给⽤户⼦程序VUMAT的材料点结构保持不变；删除的材料点没有从块中移⾛。
ABAQUS/Explicit将传递0应⼒及应变给所有删除的材料点。⼀旦⼀个材料点被标⽰为删除，该材料点将不能够被再次激活。User subroutine interface ⽤户⼦程序subroutine vumat(
C Read only (unmodifiable)variables -1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,5 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only (modifiable) variables -7 stressNew, stateNew, enerInternNew, enerInelasNew )C
include 'vaba_param.inc'C
dimension props(nprops), density(nblock), coordMp(nblock,*),1 charLength(nblock), strainInc(nblock,ndir+nshr),2 relSpinInc(nblock,nshr), tempOld(nblock),3 stretchOld(nblock,ndir+nshr),4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(nblock),8 stretchNew(nblock,ndir+nshr),8 defgradNew(nblock,ndir+nshr+nshr),9 fieldNew(nblock,nfieldv),
1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),2 enerInternNew(nblock), enerInelasNew(nblock),C
character*80 cmnameC
do 100 km = 1,nblockuser coding100 continuereturnend
Variables to be defined 被定义的变量stressNew (nblock, ndir+nshr)
Stress tensor at each material point at the end of the increment.在增量结束时每个材料点的应⼒张量。stateNew (nblock, nstatev)
State variables at each material point at the end of the increment. You define the size of this array by allocating space for it(see “User subroutines: overview,” Section 25.1.1, for more information).
增量结束时每个材料点的状态变量。⽤户通过分配空间来定义该矩阵的⼤⼩（更多的资料参考“User subroutines: overview,”Section 25.1.1）。
Variables that can be updated 可以更新的变量enerInternNew (nblock)
Internal energy per unit mass at each material point at the end of the increment.增量结束时每个材料点单位质量的内能。enerInelasNew (nblock)
Dissipated inelastic energy per unit mass at each material point at the end of the increment.增量结束时每个材料点单位质量的消散的⽆弹性能。Variables passed in for informationnblock
Number of material points to be processed in this call to VUMAT.调⽤VUMAT的材料点号。ndir
Number of direct components in a symmetric tensor.对称张量的直接分量号。nshr
Number of indirect components in a symmetric tensor.对称张量的间接分量号。nstatev
Number of user-defined state variables that are associated with this material type (you define this as described in “Allocatingspace” in “User subroutines: overview,” Section 25.1.1).与材料类型相关的⽤户定义状态变量号。nfieldv
Number of user-defined external field variables.⽤户定义外部场变量号。nprops
User-specified number of user-defined material properties.⽤户定义材料属性的⽤户指定号。lanneal
Flag indicating whether the routine is being called during an annealing process. lanneal=0indicates that the routine is being called during a normal mechanics increment. lanneal=1
indicates that this is an annealing process and you should re-initialize the internal state variables, stateNew, if necessary.ABAQUS/Explicit will automatically set the stresses, stretches, and state to a value of zero during the annealing process.在退⽕处理过程中标⽰程序是否被调⽤。Laneal=0表明程序在正常⼒学增量过程中被调⽤。Laneal＝1表明这是⼀个退⽕过
程，并且如果需要的话⽤户需要重新初始化内部状态变量stateNew。ABAQUS/Explicit将⾃动设置应⼒，延展性及在退⽕过程中0值状态。stepTime
Value of time since the step began.从时间步开始时的的时间totalTime
Value of total time. The time at the beginning of the step is given by totalTime - stepTime.总时间值。时间步开始时的时间定义为totalTime - stepTimedt
Time increment size.时间增量⼤⼩。cmname
User-specified material name, left justified. It is passed in as an upper-case character string.
Some internal material models are given names starting with the “ABQ_” character string. To avoid conflict, you should notuse “ABQ_” as the leading string for cmname.
⽤户指定材料名。按照⼤写字母传⼊。有些内部材料本构以ABQ_字母开头赋名。为了避免冲突，⽤户不能使⽤ABQ_作为cmname的开头字母。coordMp(nblock,*)
Material point coordinates. It is the midplane material point for shell elements and the centroid for beam elements.材料点坐标。对于壳单元为中平⾯材料点，对于梁单元为质⼼。charLength(nblock)
Characteristic element length. This is a typical length of a line across an element. For beams and trusses, it is a characteristiclength along the element axis. For membranes and shells, it is a characteristic length in the reference surface. For
axisymmetric elements, it is a characteristic length in the –plane only. For cohesive elements it is equal to the constitutivethickness.
特征单元长度。它是⼀个穿过单元线的特殊长度。对于梁及桁架，它表⽰沿单元轴的特征长度。对于膜及壳，它表⽰参考⾯上的特征长度。对于对称单元，它仅仅表⽰–平⾯上的特征长度。对于粘聚单元它等于结构厚度。props(nprops)
User-supplied material properties.⽤户指定材料属性density(nblock)
Current density at the material points in the midstep configuration. This value may be inaccurate in problems where the
volumetric strain increment is very small. If an accurate value of the density is required in such cases, the analysis should berun in double precision. This value of the density is not affected by mass scaling.
材料点的当前密度。在体积应变增量⾮常⼩的情况下，该值可能出现精确性问题。在这种情况下如果需要精确的密度值，分析需要在双精度下运⾏。密度值不受质量缩放⽐例的影响strainInc (nblock, ndir+nshr)
Strain increment tensor at each material point.每个材料点的应变增量张量。relSpinInc (nblock, nshr)
Incremental relative rotation vector at each material point defined in the corotational system.
Defined as , where is the antisymmetric part of the velocity gradient, , and
. Stored in 3-D as and in 2-D as .
定义在联合旋转系中的每个材料点的增量相关选择⽮量。定义为，其中为速度梯度，的反对称部分，在三维条件下存储为，在⼆维条件下存储为。tempOld(nblock)
Temperatures at each material point at the beginning of the increment.增量开始时每个材料点的温度。stretchOld (nblock, ndir+nshr)
Stretch tensor, , at each material point at the beginning of the increment defined from the
polar decomposition of the deformation gradient by .
增量开始时每个材料点的，拉伸张量。defgradOld (nblock,ndir+2*nshr)
Deformation gradient tensor at each material point at the beginning of the increment. Stored in
3-D as (, , , , , , , , ) and in 2-D as (, , , ,
).
增量开始时每个材料点的位移梯度张量。在三维状态下存储为(, , , , ,
, , , ) 及在⼆维条件下存储为(, , , , ).fieldOld (nblock, nfieldv)
Values of the user-defined field variables at each material point at the beginning of theincrement.
增量开始时每个材料点的⽤户定义场变量值。stressOld (nblock, ndir+nshr)
Stress tensor at each material point at the beginning of the increment.增量开始时每个材料点的应⼒张量。stateOld (nblock, nstatev)
State variables at each material point at the beginning of the increment.增量开始时每个材料点的状态变量。enerInternOld (nblock)
Internal energy per unit mass at each material point at the beginning of the increment.增量开始时每个材料点的单位质量的内能enerInelasOld (nblock)
Dissipated inelastic energy per unit mass at each material point at the beginning of theincrement.
增量开始时每个材料点的单位质量消散的⽆弹性能。tempNew(nblock)
Temperatures at each material point at the end of the increment.增量开始时每个材料点的温度。stretchNew (nblock, ndir+nshr)
Stretch tensor, , at each material point at the end of the increment defined from the polar
decomposition of the deformation gradient by .
defgradNew (nblock,ndir+2*nshr)
Deformation gradient tensor at each material point at the end of the increment. Stored in 3-D as
(, , , , , , , , ) and in 2-D as (, , , , ).
增量结束时每个材料点的位移剃度张量。三维状态下存储为(, , , , , ,
, , ) ，⼆维状态下存储为(, , , , )。fieldNew (nblock, nfieldv)
Values of the user-defined field variables at each material point at the end of the increment.增量结束时每个材料点的⽤户指定场变量值。
Example: Using more than one user-defined material model例⼦：使⽤多个⽤户定义材料模型
To use more than one user-defined material model, the variable cmname can be tested for different material names insideuser subroutine VUMAT, as illustrated below:
要使⽤多种⽤户材料模型，⽤户⼦程序VUMAT中使⽤cmname值代表不同的材料名，如下所述：if (cmname(1:4) .eq. 'MAT1') thencall VUMAT_MAT1(argument_list)else if (cmname(1:4) .eq. 'MAT2') thencall VUMAT_MAT2(argument_list)end if
VUMAT_MAT1 and VUMAT_MAT2 are the actual user material subroutines containing the constitutive material models foreach material MAT1 and MAT2, respectively. Subroutine VUMAT merely acts as a directory here. The argument list can bethe same as that used in subroutine VUMAT. The material names must be in upper case since cmname is passed in as an
upper-case character string.
VUMAT_MAT1 and VUMAT_MAT2分别为实际的包含MAT1及MAT2材料本构的⽤户材料⼦程序。⼦程序VUMAT在此仅仅作为⼀个地址。
Example: Elastic/plastic material with kinematic hardening例⼦：带有运动硬化的弹性/塑性材料
As a simple example of the coding of subroutine VUMAT, consider the generalized plane strain case for an elastic/plasticmaterial with kinematic hardening. The basic assumptions and definitions of the model are as follows.
作为⼦程序VUMAT的简单例⼦代码，考虑通⽤运动硬化弹性/塑性材料的平⾯应变问题。模型的基本假设及定义如下。Let be the current value of the stress and define to be the deviatoric part of the stress. The
center of the yield surface in deviatoric stress space is given by the tensor , which has initial
values of zero. The stress difference, , is the stress measured from the center of the yield surface
and is given by
设定为应⼒的当前值，并定义应⼒的偏量部分。在偏应⼒空间屈服⾯中⼼由张量给定，其初始值为0。应⼒导数为从屈服⾯中学测量的应⼒，由给出。
The von Mises yield surface is defined as Mises屈服⾯定义为
where is the uniaxial equivalent yield stress. The von Mises yield surface is a cylinder in
deviatoric stress space with a radius of
为单轴等效屈服应⼒，在偏应⼒空间中Mises屈服⾯是⼀个以为半径的圆柱。For the kinematic hardening model, is aconstant. The normal to the Mises yield surface can be
written as
对于运动硬化模型，是⼀个常数。Mises屈服⾯的法线定义为
We decompose the strain rate into an elastic and plastic part using an additive decomposition:
将应变率分解为弹性和塑性部分：
The plastic part of the strain rate is given by a normality condition塑性应变率由给出
where the scalar multiplier must be determined. A scalar measure of equivalent plastic strain rate
is defined by
式中必须指定⽐例乘⼦。等效塑性应变率⽐例量测定义为
The stress rate is assumed to be purely due to the elastic part of the strain rate and is expressed in terms of Hooke's law by
where and are the Lamés constants for the material.
假定应⼒率完全是由弹性应变率引起，并且由Hooke定律表⽰为式中及是材料的拉枚常数。
The evolution law for is given as
where is the slope of the uniaxial yield stress versus plastic strain curve.
的发展准则定义为
式中，单轴屈服应⼒与塑性应变曲线的斜率
During active plastic loading the stress must remain on the yield surface, so that在激活塑性荷载是应⼒必须保持在屈服⾯上，所以
The equivalent plastic strain rate is related to by等效塑性应变率通过式与相关。
The kinematic hardening constitutive model is integrated in a rate form as follows. A trial elastic stress is computed as
where the subscripts and refer to the beginning and end of the increment, respectively.
If the trial stress does not exceed the yield stress, the new stress is set equal to the trial stress. If the yield stress is exceeded,plasticity occurs in the increment. We then write the incremental analogs of the rate equations as
where
From the definition of the normal to the yield surface at the end of the increment, ,
This can be expanded using the incremental equations as
Taking the tensor product of this equation with , using the the yield condition at the end of the
increment, and solving for :
The value for is used in the incremental equations to determine , , and .
This algorithm is often referred to as an elastic predictor, radial return algorithm because the correction to the trial stress
under the active plastic loading condition returns the stress state to the yield surface along the direction defined by the vectorfrom the center of the yield surface to the elastic trial stress. The subroutine would be coded as follows:subroutine vumat(C Read only –只读
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,2 stepTime, totalTime, dt, cmname, coordMp, charLength,3 props, density, strainInc, relSpinInc,4 tempOld, stretchOld, defgradOld, fieldOld,3 stressOld, stateOld, enerInternOld, enerInelasOld,6 tempNew, stretchNew, defgradNew, fieldNew,C Write only –只写
5 stressNew, stateNew, enerInternNew, enerInelasNew )C
include 'vaba_param.inc'C
C J2 Mises Plasticity with kinematic hardening for planeC strain case.平⾯应变状态下带有运动硬化的Mises塑性C Elastic predictor, radial corrector algorithm.C
C The state variables are stored as:
C STATE(*,1) = back stress component 11C STATE(*,2) = back stress component 22C STATE(*,3) = back stress component 33C STATE(*,4) = back stress component 12
C STATE(*,5) = equivalent plastic strain 等效塑性应变CC
C All arrays dimensioned by (*) are not used in this algorithmdimension props(nprops), density(nblock),1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),3 relSpinInc(*), tempOld(*),4 stretchOld(*), defgradOld(*),
5 fieldOld(*), stressOld(nblock,ndir+nshr),6 stateOld(nblock,nstatev), enerInternOld(nblock),7 enerInelasOld(nblock), tempNew(*),8 stretchNew(*), defgradNew(*), fieldNew(*),
9 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),1 enerInternNew(nblock), enerInelasNew(nblock)C
character*80 cmnameC
parameter( zero = 0., one = 1., two = 2., three = 3.,1 third = one/three, half = .5, twoThirds = two/three,2 threeHalfs = 1.5 )C
e = props(1)xnu = props(2)yield = props(3)hard = props(4)C
twomu = e / ( one + xnu )thremu = threeHalfs * twomusixmu = three * twomu
alamda = twomu * ( e - twomu ) / ( sixmu - two * e )
term = one / ( twomu * ( one + hard/thremu ) )con1 = sqrt( twoThirds )C
do 100 i = 1,nblockC
C Trial stress
trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)sig1 = stressOld(i,1) + alamda*trace + twomu*strainInc(i,1)sig2 = stressOld(i,2) + alamda*trace + twomu*strainInc(i,2)sig3 = stressOld(i,3) + alamda*trace + twomu*strainInc(i,3)sig4 = stressOld(i,4) + twomu*strainInc(i,4) CC Trial stress measured from the back stresss1 = sig1 - stateOld(i,1)s2 = sig2 - stateOld(i,2)s3 = sig3 - stateOld(i,3)s4 = sig4 - stateOld(i,4)C
C Deviatoric part of trial stress measured from the back stresssmean = third * ( s1 + s2 + s3 )ds1 = s1 - smeands2 = s2 - smeands3 = s3 - smeanC
C Magnitude of the deviatoric trial stress differencedsmag = sqrt( ds1**2 + ds2**2 + ds3**2 + 2.*s4**2 )C
C Check for yield by determining the factor for plasticity,C zero for elastic, one for yieldradius = con1 * yieldfacyld = zero
if( dsmag - radius .ge. zero ) facyld = oneC
C Add a protective addition factor to prevent a divide by zeroC when dsmag is zero. If dsmag is zero, we will not have exceededC the yield stress and facyld will be zero.
dsmag = dsmag + ( one - facyld )C
C Calculated increment in gamma (this explicitly includes theC time step)diff = dsmag - radiusdgamma = facyld * term * diffC
C Update equivalent plastic straindeqps = con1 * dgamma
stateNew(i,5) = stateOld(i,5) + deqpsC
C Divide dgamma by dsmag so that the deviatoric stresses areC explicitly converted to tensors of unit magnitude in theC following calculationsdgamma = dgamma / dsmagC
C Update back stress
factor = hard * dgamma * twoThirdsstateNew(i,1) = stateOld(i,1) + factor * ds1stateNew(i,2) = stateOld(i,2) + factor * ds2stateNew(i,3) = stateOld(i,3) + factor * ds3stateNew(i,4) = stateOld(i,4) + factor * s4C
C Update the stressfactor = twomu * dgammastressNew(i,1) = sig1 - factor * ds1stressNew(i,2) = sig2 - factor * ds2stressNew(i,3) = sig3 - factor * ds3stressNew(i,4) = sig4 - factor * s4C
C Update the specific internal energy -stressPower = half * (
1 ( stressOld(i,1)+stressNew(i,1) )*strainInc(i,1)1 + ( stressOld(i,2)+stressNew(i,2) )*strainInc(i,2)1 + ( stressOld(i,3)+stressNew(i,3) )*strainInc(i,3)
1 + two*( stressOld(i,4)+stressNew(i,4) )*strainInc(i,4) ) CenerInternNew(i) = enerInternOld(i)1 + stressPower / density(i)C
C Update the dissipated inelastic specific energy -smean = third * ( stressNew(i,1) + stressNew(i,2)1 + stressNew(i,3) )
equivStress = sqrt( threeHalfs *1 ( (stressNew(i,1)-smean)**21 + (stressNew(i,2)-smean)**21 + (stressNew(i,3)-smean)**21 + two * stressNew(i,4)**2 ) )
plasticWorkInc = equivStress * deqpsenerInelasNew(i) = enerInelasOld(i)1 + plasticWorkInc / density(i)100 continueCreturnend
因篇幅问题不能全部显示，请点此查看更多更全内容