ProductsAbaqus/Standard Storage of stress and strain componentsIn the stress and strain arrays and in the matrices DDSDDE, DDSDDT, and DRPLDE, direct components are stored first, followed by shear components. There are NDI direct and NSHR engineering shear components. The order of the components is defined in Conventions. Since the number of active stress and strain components varies between element types, the routine must be coded to provide for all element types with which it will be used. Defining local orientationsIf a local orientation (Orientations) is used at the same point as user subroutine UMAT, the stress and strain components will be in the local orientation; and, in the case of finite-strain analysis, the basis system in which stress and strain components are stored rotates with the material. StabilityYou should ensure that the integration scheme coded in this routine is stable—no direct provision is made to include a stability limit in the time stepping scheme based on the calculations in UMAT. Convergence rateDDSDDE and—for coupled temperature-displacement and coupled thermal-electrical-structural analyses—DDSDDT, DRPLDE, and DRPLDT must be defined accurately if rapid convergence of the overall Newton scheme is to be achieved. In most cases the accuracy of this definition is the most important factor governing the convergence rate. Since nonsymmetric equation solution is as much as four times as expensive as the corresponding symmetric system, if the constitutive Jacobian (DDSDDE) is only slightly nonsymmetric (for example, a frictional material with a small friction angle), it may be less expensive computationally to use a symmetric approximation and accept a slower convergence rate. An incorrect definition of the material Jacobian affects only the convergence rate; the results (if obtained) are unaffected. Viscoelastic behavior in frequency domainThe constitutive Jacobian (DDSDDE) must provide both the stiffness (storage modulus) and damping (loss modulus) for modeling frequency domain viscoelastic behavior. Special considerations for various element typesThere are several special considerations that need to be noted. Deformation gradientThe deformation gradient is available for solid (continuum) elements, membranes, and finite-strain shells (S3/S3R, S4, S4R, SAXs, and SAXAs). It is not available for beams or small-strain shells. It is stored as a 3 × 3 matrix with component equivalence DFGRD0(I,J)$\iff {F}_{IJ}$. For fully integrated first-order isoparametric elements (4-node quadrilaterals in two dimensions and 8-node hexahedra in three dimensions) the selectively reduced integration technique is used (also known as the $\overline{B}$ technique). Thus, a modified deformation gradient $$\overline{\mathbf{F}}=\mathbf{F}{\left(\frac{\overline{J}}{J}\right)}^{\frac{1}{n}}$$
is passed into user subroutine UMAT. For more details, see Solid isoparametric quadrilaterals and hexahedra. The deformation gradient, which is passed to the user subroutine, is computed with respect to the initial configuration. If a local orientation is not specified, the components of the deformation gradient are expressed in the global coordinate system. If a local orientation is used, the components of the same deformation gradient are expressed in the local coordinate system; in the case of finite-strain analysis, the basis system rotates with the material. Beams and shells that calculate transverse shear energyIf user subroutine UMAT is used to describe the material of beams or shells that calculate transverse shear energy, you must specify the transverse shear stiffness as part of the beam or shell section definition to define the transverse shear behavior. See Shell section behavior and Choosing a beam element for information on specifying this stiffness. Open-section beam elementsWhen user subroutine UMAT is used to describe the material response of beams with open sections (for example, an I-section), the torsional stiffness is obtained as $$GJ=\frac{\left({K}_{13}+{K}_{23}\right)J}{2kA},$$
where J is the torsional constant, A is the section area, k is a shear factor, and ${K}_{\alpha 3}$ is the user-specified transverse shear stiffness (see Transverse shear stiffness definition). Elements with hourglassing modesIf this capability is used to describe the material of elements with hourglassing modes, you must define the hourglass stiffness factor for hourglass control based on the total stiffness approach as part of the element section definition. The hourglass stiffness factor is not required for enhanced hourglass control, but you can define a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surface normal). See Section controls for information on specifying the stiffness factor. Pipe-soil interaction elementsThe constitutive behavior of the pipe-soil interaction elements (see Pipe-soil interaction elements) is defined by the force per unit length caused by relative displacement between two edges of the element. The relative-displacements are available as “strains” (STRAN and DSTRAN). The corresponding forces per unit length must be defined in the STRESS array. The Jacobian matrix defines the variation of force per unit length with respect to relative displacement. For two-dimensional elements two in-plane components of “stress” and “strain” exist (NTENS=NDI=2, and NSHR=0). For three-dimensional elements three components of “stress” and “strain” exist (NTENS=NDI=3, and NSHR=0). Large volume changes with geometric nonlinearityIf the material model allows large volume changes and geometric nonlinearity is considered, the exact definition of the consistent Jacobian should be used to ensure rapid convergence. These conditions are most commonly encountered when considering either large elastic strains or pressure-dependent plasticity. In the former case, total-form constitutive equations relating the Cauchy stress to the deformation gradient are commonly used; in the latter case, rate-form constitutive laws are generally used. For total-form constitutive laws, the exact consistent Jacobian $\mathbf{C}$ is defined through the variation in Kirchhoff stress: $$\delta \left(J\mathit{\sigma}\right)=J(\mathbf{C}:\delta \mathbf{D}+\delta \mathbf{W}\cdot \mathit{\sigma}-\mathit{\sigma}\cdot \delta \mathbf{W})$$
Here, J is the determinant of the deformation gradient, $\mathit{\sigma}$ is the Cauchy stress, $\delta \mathbf{D}$ is the virtual rate of deformation, and $\delta \mathbf{W}$ is the virtual spin tensor, defined as $$\delta \mathbf{D}\stackrel{\mathrm{def}}{=}\mathrm{sym}\left(\delta \mathbf{F}\cdot {\mathbf{F}}^{-1}\right)$$
and $$\delta \mathbf{W}\stackrel{\mathrm{def}}{=}\mathrm{asym}\left(\delta \mathbf{F}\cdot {\mathbf{F}}^{-1}\right).$$
For rate-form constitutive laws, the exact consistent Jacobian is given by $$\mathbf{C}=\frac{1}{J}\frac{\partial \mathrm{\Delta}\left(J\mathit{\sigma}\right)}{\partial \mathrm{\Delta}\mathit{\epsilon}}.$$ Use with almost incompressible or fully incompressible elastic materialsFor user-defined almost incompressible or incompressible elastic materials, a few different options are available depending on whether hybrid or nonhybrid elements are used. For all cases the first option should be to use user subroutine UHYPER instead of user subroutine UMAT when it is possible to do so. In user subroutine UMAT incompressible materials can be modeled via a penalty method; that is, you ensure that a finite bulk modulus is used. The bulk modulus should be large enough to model incompressibility sufficiently but small enough to avoid loss of precision. As a general guideline, the bulk modulus should be about ${10}^{4}$–${10}^{6}$ times the shear modulus. The tangent bulk modulus ${K}^{t}$ can be calculated from $${K}^{t}\mathit{}=\mathit{}\frac{1}{9}\underset{\mathrm{I}=1}{\sum ^{3}}\underset{\mathrm{J}=1}{\sum ^{3}}\mathrm{DDSDDE}(\mathrm{I},\mathrm{J}).$$
If a hybrid element is used with user subroutine UMAT, Abaqus/Standard, by default, replaces the pressure stress calculated from your definition of STRESS with that derived from the Lagrange multiplier and modifies the Jacobian appropriately (Hybrid incompressible solid element formulation). This approach is suitable for material models that use an incremental formulation (for example, metal plasticity) but is not consistent with a total formulation that is commonly used for hyperelastic materials. In the latter situation, the default formulation may lead to convergence problems. Such convergence problems may be observed, for example, when an almost incompressible nonlinear elastic user material is subjected to large deformations. Abaqus/Standard provides an alternate total formulation when user materials are used with hybrid elements (see User-defined mechanical material behavior). This formulation is consistent with the native almost incompressible formulation used by Abaqus for hyperelastic materials (Hyperelastic material behavior) and works better than the default formulation for such cases. The total hybrid formulation assumes that the response of the material can be written as the sum of its deviatoric and volumetric parts and that these parts are decoupled from each other. In particular, the volumetric response is assumed to be defined in terms of a strain energy potential, $U\left(\widehat{J}\right)$, which is a function of an alternate variable, $\widehat{J}$ in place of the actual volume change $J$. The alternate variable is made available inside user subroutine UMAT by extending the STRESS array beyond NTENS, with the NTENS+1 entry providing read access to $\widehat{J}$. You must define the hydrostatic part of the stress tensor as $\widehat{p}=-\frac{\partial U}{\partial \widehat{J}}$. The formulation also requires the additional derivatives $\widehat{K}=J\frac{{\partial}^{2}U}{\partial {\widehat{J}}^{2}}$ and $\frac{\partial \widehat{K}}{\partial \widehat{J}}=J\frac{{\partial}^{3}U}{\partial {\widehat{J}}^{3}}$. You must define these additional derivatives inside user subroutine UMAT as the NTENS+1 and the NTENS+2 entry, respectively, of the STRESS array. In addition, the bulk modulus of the material (contributes toward the material Jacobian matrix, DDSDDE) must be defined as $\widehat{K}$ Abaqus/Standard also provides a fully incompressible user material formulation for use with hybrid elements to define a fully incompressible user material response. This formulation is consistent with the native formulation used by Abaqus for incompressible hyperelastic materials and assumes that the deviatoric stress can be derived from a strain energy potential function. You need define only the deviatoric stress and Jacobian to define a fully incompressible material response through user subroutine UMAT. For incompressible pressure-sensitive materials the element choice is particularly important when using user subroutine UMAT. In particular, first-order wedge elements should be avoided. For these elements the $\overline{B}$ technique is not used to alter the deformation gradient that is passed into user subroutine UMAT, which increases the risk of volumetric locking. Increments for which only the Jacobian can be definedAbaqus/Standard passes zero strain increments into user subroutine UMAT to start the first increment of all the steps and all increments of steps for which you have suppressed extrapolation (see Defining an analysis). In this case you can define only the Jacobian (DDSDDE). Utility routinesSeveral utility routines may help in coding user subroutine UMAT. Their functions include determining stress invariants for a stress tensor and calculating principal values and directions for stress or strain tensors. These utility routines are discussed in detail in Obtaining stress invariants, principal stress/strain values and directions, and rotating tensors in an Abaqus/Standard analysis. User subroutine interfaceSUBROUTINE UMAT(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,JSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME 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),DFGRD0(3,3),DFGRD1(3,3), 4 JSTEP(4) user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT RETURN END Variables to be defined
Variables that can be updated
Variables passed in for information
Example: Using more than one user-defined mechanical material modelTo use more than one user-defined mechanical material model, the variable CMNAME can be tested for different material names inside user subroutine UMAT as illustrated below: IF (CMNAME(1:4) .EQ. 'MAT1') THEN CALL UMAT_MAT1(argument_list) ELSE IF(CMNAME(1:4) .EQ. 'MAT2') THEN CALL UMAT_MAT2(argument_list) END IF UMAT_MAT1 and UMAT_MAT2 are the actual user material subroutines containing the constitutive material models for each material MAT1 and MAT2, respectively. Subroutine UMAT merely acts as a directory here. The argument list may be the same as that used in subroutine UMAT. Example: Simple linear viscoelastic materialAs a simple example of the coding of user subroutine UMAT, consider the linear, viscoelastic model shown in Figure 1. Although this is not a very useful model for real materials, it serves to illustrate how to code the routine. Figure 1. Simple linear viscoelastic model.
The behavior of the one-dimensional model shown in the figure is $$\sigma +\frac{{\mu}_{1}}{\left({E}_{1}+{E}_{2}\right)}\dot{\sigma}=\frac{{\mu}_{1}}{\left(1+{E}_{1}/{E}_{2}\right)}\dot{\epsilon}+\frac{1}{\left(1/{E}_{1}+1/{E}_{2}\right)}\epsilon ,$$
where $\dot{\sigma}$ and $\dot{\epsilon}$ are the time rates of change of stress and strain. This can be generalized for small straining of an isotropic solid as $${\sigma}_{xx}+\stackrel{~}{\nu}{\dot{\sigma}}_{xx}=\lambda {\epsilon}_{V}+2\mu {\epsilon}_{xx}+\stackrel{~}{\lambda}{\dot{\epsilon}}_{V}+2\stackrel{~}{\mu}{\dot{\epsilon}}_{xx},\text{etc.,}$$
and $${\sigma}_{xy}+\stackrel{~}{\nu}{\dot{\sigma}}_{xy}=\mu {\gamma}_{xy}+\stackrel{~}{\mu}{\dot{\gamma}}_{xy},\text{etc.,}$$
where $${\epsilon}_{V}={\epsilon}_{xx}+{\epsilon}_{yy}+{\epsilon}_{zz},$$
and $\stackrel{~}{\nu}$, $\lambda $, $\mu $, $\stackrel{~}{\lambda}$, and $\stackrel{~}{\mu}$ are material constants ($\lambda $ and $\mu $ are the Lamé constants). A simple, stable integration operator for this equation is the central difference operator: $$\begin{array}{cc}\hfill {\dot{f}}_{t+\frac{1}{2}\mathrm{\Delta}t}& =\frac{\mathrm{\Delta}f}{\mathrm{\Delta}t},\hfill \\ \hfill {f}_{t+\frac{1}{2}\mathrm{\Delta}t}& ={f}_{t}+\frac{\mathrm{\Delta}f}{2},\hfill \end{array}$$
where f is some function, ${f}_{t}$ is its value at the beginning of the increment, $\mathrm{\Delta}f$ is the change in the function over the increment, and $\mathrm{\Delta}t$ is the time increment. Applying this to the rate constitutive equations above gives $$\left(\frac{\mathrm{\Delta}t}{2}+\stackrel{~}{\nu}\right)\mathrm{\Delta}{\sigma}_{xx}=\left(\mathrm{\Delta}t\frac{\lambda}{2}+\stackrel{~}{\lambda}\right)\mathrm{\Delta}{\epsilon}_{V}+\left(\mathrm{\Delta}t\mu +2\stackrel{~}{\mu}\right)\mathrm{\Delta}{\epsilon}_{xx}+\mathrm{\Delta}t{\left(\lambda {\epsilon}_{V}+2\mu {\epsilon}_{xx}-{\sigma}_{xx}\right)}_{t},\text{etc.,}$$
and $$\left(\frac{\mathrm{\Delta}t}{2}+\stackrel{~}{\nu}\right)\mathrm{\Delta}{\sigma}_{xy}=\left(\mathrm{\Delta}t\frac{\mu}{2}+\stackrel{~}{\mu}\right)\mathrm{\Delta}{\gamma}_{xy}+\mathrm{\Delta}t{\left(\mu {\gamma}_{xy}-{\sigma}_{xy}\right)}_{t},\text{etc.,}$$
so that the Jacobian matrix has the terms $$\frac{\partial \mathrm{\Delta}{\sigma}_{xx}}{\partial \mathrm{\Delta}{\epsilon}_{xx}}=\frac{1}{\left(\frac{\mathrm{\Delta}t}{2}+\stackrel{~}{\nu}\right)}\left[\mathrm{\Delta}t\left(\frac{\lambda}{2}+\mu \right)+\stackrel{~}{\lambda}+2\stackrel{~}{\mu}\right],$$
$$\frac{\partial \mathrm{\Delta}{\sigma}_{xx}}{\partial \mathrm{\Delta}{\epsilon}_{yy}}=\frac{1}{\left(\frac{\mathrm{\Delta}t}{2}+\stackrel{~}{\nu}\right)}\left[\mathrm{\Delta}t\frac{\lambda}{2}+\stackrel{~}{\lambda}\right],$$
and $$\frac{\partial \mathrm{\Delta}{\sigma}_{xy}}{\partial \mathrm{\Delta}{\gamma}_{xy}}=\frac{1}{\left(\frac{\mathrm{\Delta}t}{2}+\stackrel{~}{\nu}\right)}\left[\mathrm{\Delta}t\frac{\mu}{2}+\stackrel{~}{\mu}\right].$$
The total change in specific energy in an increment for this material is $$\left({\sigma}_{ij}+\frac{1}{2}\mathrm{\Delta}{\sigma}_{ij}\right)\mathrm{\Delta}{\epsilon}_{ij},$$
while the change in specific elastic strain energy is $$\left({\epsilon}_{ij}+\frac{1}{2}\mathrm{\Delta}{\epsilon}_{ij}\right){D}_{ijkl}\mathrm{\Delta}{\epsilon}_{kl},$$
where D is the elasticity matrix: $$\left[\begin{array}{cccccc}\hfill \lambda +2\mu \hfill & \hfill \lambda \hfill & \hfill \lambda \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \lambda \hfill & \hfill \lambda +2\mu \hfill & \hfill \lambda \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \lambda \hfill & \hfill \lambda \hfill & \hfill \lambda +2\mu \hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mu \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mu \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mu \hfill \end{array}\right].$$
No state variables are needed for this material, so the allocation of space for them is not necessary. In a more realistic case a set of parallel models of this type might be used, and the stress components in each model might be stored as state variables. For our simple case a user material definition can be used to read in the five constants in the order $\lambda $, $\mu $, $\stackrel{~}{\lambda}$, $\stackrel{~}{\mu}$, and $\stackrel{~}{\nu}$ so that $$\mathtt{\text{PROPS(1)}}=\lambda ,$$
$$\mathtt{\text{PROPS(2)}}=\mu ,$$
$$\mathtt{\text{PROPS(3)}}=\stackrel{~}{\lambda},$$
$$\mathtt{\text{PROPS(4)}}=\stackrel{~}{\mu},$$
$$\mathtt{\text{PROPS(5)}}=\stackrel{~}{\nu}.$$
The routine can then be coded as follows: SUBROUTINE UMAT(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,JSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C 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), 5 JSTEP(4) DIMENSION DSTRES(6),D(3,3) C C EVALUATE NEW STRESS TENSOR C EV = 0. DEV = 0. DO K1=1,NDI EV = EV + STRAN(K1) DEV = DEV + DSTRAN(K1) END DO C TERM1 = .5*DTIME + PROPS(5) TERM1I = 1./TERM1 TERM2 = (.5*DTIME*PROPS(1)+PROPS(3))*TERM1I*DEV TERM3 = (DTIME*PROPS(2)+2.*PROPS(4))*TERM1I C DO K1=1,NDI DSTRES(K1) = TERM2+TERM3*DSTRAN(K1) 1 +DTIME*TERM1I*(PROPS(1)*EV 2 +2.*PROPS(2)*STRAN(K1)-STRESS(K1)) STRESS(K1) = STRESS(K1) + DSTRES(K1) END DO C TERM2 = (.5*DTIME*PROPS(2) + PROPS(4))*TERM1I I1 = NDI DO K1=1,NSHR I1 = I1+1 DSTRES(I1) = TERM2*DSTRAN(I1)+ 1 DTIME*TERM1I*(PROPS(2)*STRAN(I1)-STRESS(I1)) STRESS(I1) = STRESS(I1)+DSTRES(I1) END DO C C CREATE NEW JACOBIAN C TERM2 = (DTIME*(.5*PROPS(1)+PROPS(2))+PROPS(3)+ 1 2.*PROPS(4))*TERM1I TERM3 = (.5*DTIME*PROPS(1)+PROPS(3))*TERM1I DO K1=1,NTENS DO K2=1,NTENS DDSDDE(K2,K1) = 0. END DO END DO C DO K1=1,NDI DDSDDE(K1,K1) = TERM2 END DO C DO K1=2,NDI N2 = K1−1 DO K2=1,N2 DDSDDE(K2,K1) = TERM3 DDSDDE(K1,K2) = TERM3 END DO END DO TERM2 = (.5*DTIME*PROPS(2)+PROPS(4))*TERM1I I1 = NDI DO K1=1,NSHR I1 = I1+1 DDSDDE(I1,I1) = TERM2 END DO C C TOTAL CHANGE IN SPECIFIC ENERGY C TDE = 0. DO K1=1,NTENS TDE = TDE + (STRESS(K1)-.5*DSTRES(K1))*DSTRAN(K1) END DO C C CHANGE IN SPECIFIC ELASTIC STRAIN ENERGY C TERM1 = PROPS(1) + 2.*PROPS(2) DO K1=1,NDI D(K1,K1) = TERM1 END DO DO K1=2,NDI N2 = K1-1 DO K2=1,N2 D(K1,K2) = PROPS(1) D(K2,K1) = PROPS(1) END DO END DO DEE = 0. DO K1=1,NDI TERM1 = 0. TERM2 = 0. DO K2=1,NDI TERM1 = TERM1 + D(K1,K2)*STRAN(K2) TERM2 = TERM2 + D(K1,K2)*DSTRAN(K2) END DO DEE = DEE + (TERM1+.5*TERM2)*DSTRAN(K1) END DO I1 = NDI DO K1=1,NSHR I1 = I1+1 DEE = DEE + PROPS(2)*(STRAN(I1)+.5*DSTRAN(I1))*DSTRAN(I1) END DO SSE = SSE + DEE SCD = SCD + TDE − DEE RETURN END |