2023年6月23日发(作者：)

1.1.41 UMATUser subroutine to define a material's mechanical t: Abaqus/Standard

Warning: The use of this subroutine generally requires considerable expertise. You arecautioned that the implementation of any realistic constitutive model requires extensivedevelopment and testing. Initial testing on a single-element model with prescribedtraction loading is strongly nces“User-defined mechanical material behavior,” Section 26.7.1 of the Abaqus Analysis User's Guide“User-defined thermal material behavior,” Section 26.7.2 of the Abaqus Analysis User's Guide*USER MATERIAL“SDVINI,” Section 4.1.11 of the Abaqus Verification Guide“UMAT and UHYPER,” Section 4.1.21 of the Abaqus Verification GuideOverviewUser subroutine UMAT:can be used to define the mechanical constitutive behavior of a material;will be called at all material calculation points of elements for which the material definition includes auser-defined material behavior;can be used with any procedure that includes mechanical behavior;can use solution-dependent state variables;must update the stresses and solution-dependent state variables to their values at the end of theincrement for which it is called;must provide the material Jacobian matrix, , for the mechanical constitutive model;can be used in conjunction with user subroutine USDFLD to redefine any field variables before they arepassed in; andis described further in “User-defined mechanical material behavior,” Section 26.7.1 of the AbaqusAnalysis User's e of stress and strain componentsIn the stress and strain arrays and in the matrices DDSDDE, DDSDDT, and DRPLDE, direct components are storedfirst, followed by shear components. There are NDI direct and NSHR engineering shear components. The orderof the components is defined in “Conventions,” Section 1.2.2 of the Abaqus Analysis User's Guide. Since thenumber of active stress and strain components varies between element types, the routine must be coded toprovide for all element types with which it will be ng local orientationsIf a local orientation (“Orientations,” Section 2.2.5 of the Abaqus Analysis User's Guide) is used at the samepoint as user subroutine UMAT, the stress and strain components will be in the local orientation; and, in the caseof finite-strain analysis, the basis system in which stress and strain components are stored rotates with ityYou should ensure that the integration scheme coded in this routine is stable—no direct provision is made toinclude a stability limit in the time stepping scheme based on the calculations in gence 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 Newtonscheme is to be achieved. In most cases the accuracy of this definition is the most important factor governingthe convergence rate. Since nonsymmetric equation solution is as much as four times as expensive as thecorresponding symmetric system, if the constitutive Jacobian (DDSDDE) is only slightly nonsymmetric (forexample, a frictional material with a small friction angle), it may be less expensive computationally to use asymmetric approximation and accept a slower convergence incorrect definition of the material Jacobian affects only the convergence rate; the results (if obtained) l considerations for various element typesThere are several special considerations that need to be bility of 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) . For fully integrated first-order isoparametricelements (4-node quadrilaterals in two dimensions and 8-node hexahedra in three dimensions) the selectivelyreduced integration technique is used (also known as the technique). Thus, a modified deformation gradientis passed into user subroutine UMAT. For more details, see “Solid isoparametric quadrilaterals and hexahedra,”Section 3.2.4 of the Abaqus Theory and shells that calculate transverse shear energyIf user subroutine UMAT is used to describe the material of beams or shells that calculate transverse shearenergy, you must specify the transverse shear stiffness as part of the beam or shell section definition to definethe transverse shear behavior. See “Shell section behavior,” Section 29.6.4 of the Abaqus Analysis User'sGuide, and “Choosing a beam element,” Section 29.3.3 of the Abaqus Analysis User's Guide, for informationon specifying this -section beam elementsWhen user subroutine UMAT is used to describe the material response of beams with open sections (forexample, an I-section), the torsional stiffness is obtained aswhere J is the torsional constant, A is the section area, k is a shear factor, and is the user-specifiedtransverse shear stiffness (see “Transverse shear stiffness definition” in “Choosing a beam element,” Section29.3.3 of the Abaqus Analysis User's Guide).Elements with hourglassing modesIf this capability is used to describe the material of elements with hourglassing modes, you must define thehourglass stiffness factor for hourglass control based on the total stiffness approach as part of the elementsection definition. The hourglass stiffness factor is not required for enhanced hourglass control, but you candefine a scaling factor for the stiffness associated with the drill degree of freedom (rotation about the surfacenormal). See “Section controls,” Section 27.1.4 of the Abaqus Analysis User's Guide, for information onspecifying the stiffness -soil interaction elementsThe constitutive behavior of the pipe-soil interaction elements (see “Pipe-soil interaction elements,” Section32.12.1 of the Abaqus Analysis User's Guide) is defined by the force per unit length caused by relativedisplacement between two edges of the element. The relative-displacements are available as “strains” (STRANand DSTRAN). The corresponding forces per unit length must be defined in the STRESS array. The Jacobianmatrix defines the variation of force per unit length with respect to relative two-dimensional elements two in-plane components of “stress” and “strain” exist (NTENS=NDI=2, andNSHR=0). For three-dimensional elements three components of “stress” and “strain” exist (NTENS=NDI=3, andNSHR=0).Large volume changes with geometric nonlinearityIf the material model allows large volume changes and geometric nonlinearity is considered, the exact definitionof the consistent Jacobian should be used to ensure rapid convergence. These conditions are most commonlyencountered 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; inthe latter case, rate-form constitutive laws are generally total-form constitutive laws, the exact consistent Jacobian is defined through the variation in Kirchhoffstress:Here, J is the determinant of the deformation gradient, is the Cauchy stress,

deformation, and is the virtual spin tensor, defined as is the virtual rate ofandFor rate-form constitutive laws, the exact consistent Jacobian is given byUse with incompressible elastic materialsFor user-defined incompressible elastic materials, user subroutine UHYPER should be used rather than usersubroutine UMAT. In UMAT incompressible materials must be modeled via a penalty method; that is, you mustensure that a finite bulk modulus is used. The bulk modulus should be large enough to model incompressibilitysufficiently but small enough to avoid loss of precision. As a general guideline, the bulk modulus should beabout – times the shear modulus. The tangent bulk modulus can be calculated fromIf a hybrid element is used with user subroutine UMAT, Abaqus/Standard will replace the pressure stresscalculated from your definition of STRESS with that derived from the Lagrange multiplier and will modify theJacobian incompressible pressure-sensitive materials the element choice is particularly important when using usersubroutine UMAT. In particular, first-order wedge elements should be avoided. For these elements the

technique is not used to alter the deformation gradient that is passed into user subroutine UMAT, which increasesthe risk of volumetric ents for which only the Jacobian can be definedAbaqus/Standard passes zero strain increments into user subroutine UMAT to start the first increment of all thesteps and all increments of steps for which you have suppressed extrapolation (see “Defining an analysis,”Section 6.1.2 of the Abaqus Analysis User's Guide). 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 stressinvariants for a stress tensor and calculating principal values and directions for stress or strain tensors. Theseutility routines are discussed in detail in “Obtaining stress invariants, principal stress/strain values and directions,and rotating tensors in an Abaqus/Standard analysis,” Section 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,KSTEP,KINC)C INCLUDE 'ABA_'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) user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT

RETURN ENDVariables to be definedIn all situationsDDSDDE(NTENS,NTENS)Jacobian matrix of the constitutive model, , where are the stress increments and arethe strain increments. DDSDDE(I,J) defines the change in the Ith stress component at the end of the timeincrement caused by an infinitesimal perturbation of the Jth component of the strain increment you invoke the unsymmetric equation solution capability for the user-defined material,Abaqus/Standard will use only the symmetric part of DDSDDE. The symmetric part of the matrix iscalculated by taking one half the sum of the matrix and its (NTENS)This array is passed in as the stress tensor at the beginning of the increment and must be updated in thisroutine to be the stress tensor at the end of the increment. If you specified initial stresses (“Initial conditionsin Abaqus/Standard and Abaqus/Explicit,” Section 34.2.1 of the Abaqus Analysis User's Guide), thisarray will contain the initial stresses at the start of the analysis. The size of this array depends on the valueof NTENS as defined below. In finite-strain problems the stress tensor has already been rotated to accountfor rigid body motion in the increment before UMAT is called, so that only the corotational part of the stressintegration should be done in UMAT. The measure of stress used is “true” (Cauchy) (NSTATV)An array containing the solution-dependent state variables. These are passed in as the values at thebeginning of the increment unless they are updated in user subroutines USDFLD or UEXPAN, in which casethe updated values are passed in. In all cases STATEV must be returned as the values at the end of theincrement. The size of the array is defined as described in “Allocating space” in “User subroutines:overview,” Section 18.1.1 of the Abaqus Analysis User's finite-strain problems any vector-valued or tensor-valued state variables must be rotated to account forrigid body motion of the material, in addition to any update in the values associated with constitutivebehavior. The rotation increment matrix, DROT, is provided for this , SPD, SCDSpecific elastic strain energy, plastic dissipation, and “creep” dissipation, respectively. These are passed inas the values at the start of the increment and should be updated to the corresponding specific energyvalues at the end of the increment. They have no effect on the solution, except that they are used forenergy in a fully coupled thermal-stress or a coupled thermal-electrical-structural analysisRPLVolumetric heat generation per unit time at the end of the increment caused by mechanical working of (NTENS)Variation of the stress increments with respect to the (NTENS)Variation of RPL with respect to the strain Variation of RPL with respect to the in a geostatic stress procedure or a coupled pore fluid diffusion/stress analysis for porepressure cohesive elementsRPLRPL is used to indicate whether or not a cohesive element is open to the tangential flow of pore fluid. SetRPL equal to 0 if there is no tangential flow; otherwise, assign a nonzero value to RPL if an element is opened, a cohesive element will remain open to the fluid le that can be updatedPNEWDTRatio of suggested new time increment to the time increment being used (DTIME, see discussion later inthis section). This variable allows you to provide input to the automatic time incrementation algorithms inAbaqus/Standard (if automatic time incrementation is chosen). For a quasi-static procedure the automatictime stepping that Abaqus/Standard uses, which is based on techniques for integrating standard creeplaws (see “Quasi-static analysis,” Section 6.2.5 of the Abaqus Analysis User's Guide), cannot becontrolled from within the UMAT is set to a large value before each call to PNEWDT is redefined to be less than 1.0, Abaqus/Standard must abandon the time increment andattempt it again with a smaller time increment. The suggested new time increment provided to theautomatic time integration algorithms is PNEWDT × DTIME, where the PNEWDT used is the minimum valuefor all calls to user subroutines that allow redefinition of PNEWDT for this PNEWDT is given a value that is greater than 1.0 for all calls to user subroutines for this iteration and theincrement converges in this iteration, Abaqus/Standard may increase the time increment. The suggestednew time increment provided to the automatic time integration algorithms is PNEWDT × DTIME, where thePNEWDT used is the minimum value for all calls to user subroutines for this automatic time incrementation is not selected in the analysis procedure, values of PNEWDT that aregreater than 1.0 will be ignored and values of PNEWDT that are less than 1.0 will cause the job to les passed in for informationSTRAN(NTENS)An array containing the total strains at the beginning of the increment. If thermal expansion is included inthe same material definition, the strains passed into UMAT are the mechanical strains only (that is, thethermal strains computed based upon the thermal expansion coefficient have been subtracted from thetotal strains). These strains are available for output as the “elastic” finite-strain problems the strain components have been rotated to account for rigid body motion in theincrement before UMAT is called and are approximations to logarithmic (NTENS)Array of strain increments. If thermal expansion is included in the same material definition, these are themechanical strain increments (the total strain increments minus the thermal strain increments).TIME(1)Value of step time at the beginning of the current increment or (2)Value of total time at the beginning of the current ime mperature at the start of the ncrement of Array of interpolated values of predefined field variables at this point at the start of the increment, basedon the values read in at the rray of increments of predefined field User-defined material name, left justified. Some internal material models are given names starting with the“ABQ_” character string. To avoid conflict, you should not use “ABQ_” as the leading string for ber of direct stress components at this mber of engineering shear stress components at this ize of the stress or strain component array (NDI + NSHR).NSTATVNumber of solution-dependent state variables that are associated with this material type (defined asdescribed in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus AnalysisUser's Guide).PROPS(NPROPS)User-specified array of material constants associated with this user User-defined number of material constants associated with this user An array containing the coordinates of this point. These are the current coordinates if geometricnonlinearity is accounted for during the step (see “Defining an analysis,” Section 6.1.2 of the AbaqusAnalysis User's Guide); otherwise, the array contains the original coordinates of the (3,3)Rotation increment matrix. This matrix represents the increment of rigid body rotation of the basis systemin which the components of stress (STRESS) and strain (STRAN) are stored. It is provided so that vector-or tensor-valued state variables can be rotated appropriately in this subroutine: stress and straincomponents are already rotated by this amount before UMAT is called. This matrix is passed in as a unitmatrix for small-displacement analysis and for large-displacement analysis if the basis system for thematerial point rotates with the material (as in a shell element or when a local orientation is used).CELENTCharacteristic element length, which is a typical length of a line across an element for a first-order element;it is half of the same typical length for a second-order element. For beams and trusses it is a characteristiclength along the element axis. For membranes and shells it is a characteristic length in the referencesurface. For axisymmetric elements it is a characteristic length in the plane only. For cohesiveelements it is equal to the constitutive 0(3,3)Array containing the deformation gradient at the beginning of the increment. If a local orientation is definedat the material point, the deformation gradient components are expressed in the local coordinate systemdefined by the orientation at the beginning of the increment. For a discussion regarding the availability ofthe deformation gradient for various element types, see “Availability of deformation gradient.”DFGRD1(3,3)Array containing the deformation gradient at the end of the increment. If a local orientation is defined at thematerial point, the deformation gradient components are expressed in the local coordinate system definedby the orientation. This array is set to the identity matrix if nonlinear geometric effects are not included inthe step definition associated with this increment. For a discussion regarding the availability of thedeformation gradient for various element types, see “Availability of deformation gradient.”NOELElement egration point ayer number (for composite shells and layered solids).KSPTSection point number within the current tep crement e: 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 differentmaterial 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 IFUMAT_MAT1 and UMAT_MAT2 are the actual user material subroutines containing the constitutive material modelsfor each material MAT1 and MAT2, respectively. Subroutine UMAT merely acts as a directory here. The argumentlist may be the same as that used in subroutine e: Simple linear viscoelastic materialAs a simple example of the coding of user subroutine UMAT, consider the linear, viscoelastic model shown inFigure 1.1.41–1. Although this is not a very useful model for real materials, it serves to illustrate how to codethe 1.1.41–1 Simple linear viscoelastic behavior of the one-dimensional model shown in the figure iswhere and are the time rates of change of stress and strain. This can be generalized for small straining of anisotropic solid asandwhereand , , , , and are material constants ( and are the Lamé constants).A simple, stable integration operator for this equation is the central difference operator:where f is some function, is its value at the beginning of the increment,

the increment, and is the time ng this to the rate constitutive equations above gives is the change in the function overandso that the Jacobian matrix has the termsandThe total change in specific energy in an increment for this material iswhile the change in specific elastic strain energy iswhere D is the elasticity matrix:No state variables are needed for this material, so the allocation of space for them is not necessary. In a morerealistic case a set of parallel models of this type might be used, and the stress components in each model mightbe stored as state our simple case a user material definition can be used to read in the five constants in the order , , , ,and so thatThe 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,KSTEP,KINC)C INCLUDE 'ABA_'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) DIMENSION DSTRES(6),D(3,3)CC EVALUATE NEW STRESS TENSORC EV = 0. DEV = 0. DO K1=1,NDI EV = EV + STRAN(K1) DEV = DEV + DSTRAN(K1) END DOC TERM1 = .5*DTIME + PROPS(5) TERM1I = 1./TERM1 TERM2 = (.5*DTIME*PROPS(1)+PROPS(3))*TERM1I*DEV TERM3 = (DTIME*PROPS(2)+2.*PROPS(4))*TERM1IC 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 DOC 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 DOCC CREATE NEW JACOBIANC 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 DOC DO K1=1,NDI DDSDDE(K1,K1) = TERM2 END DOC 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 DOCC TOTAL CHANGE IN SPECIFIC ENERGYC TDE = 0. DO K1=1,NTENS TDE = TDE + (STRESS(K1)-.5*DSTRES(K1))*DSTRAN(K1) END DOCC CHANGE IN SPECIFIC ELASTIC STRAIN ENERGYC 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