AbaqusSubroutineDev/docs/AbaqusUserSubroutineManual/AbaqusUserSubroutineManual_023.md

c derivative d(Ni)/d(g)
dshape(1,1) = -(one - h)/four;
dshape(1,2) = (one - h)/four;
dshape(1,3) = (one + h)/four;
dshape(1,4) = -(one + h)/four;
c
c derivative d(Ni)/d(h)
dshape(2,1) = -(one - g)/four;
dshape(2,2) = -(one + g)/four;
dshape(2,3) = (one + g)/four;
dshape(2,4) = (one - g)/four;
c
c compute coordinates at the integration point
c
do k1=1, 3
    coords_ip(k1) = zero
end do
do k1=1, nnode
    do k2=1, mcrd
    coords_ip(k2)=coords_ip(k2)+shape(k1)*coords(k2,k1)
    end do
end do
c
c INTERPOLATE FIELD VARIABLES
c
if (npredf.gt.0) then
    do k1=1, npredf
    predef_loc(k1) = zero
    dpredef_loc(k1) = zero
    do k2=1, nnode
    predef_loc(k1) =
    &    predef_loc(k1)+
    &    (predef(1,k1,k2)-predef(2,k1,k2))*shape(k2)
    dpredef_loc(k1) =
    &    dpredef_loc(k1)+predef(2,k1,k2)*shape(k2)
    end do
    end do
    end if
c
c FORM B-MATRIX
c 

djac = one

do i = 1, ndim
    do j = 1, ndim
    xjac(i,j) = zero
    xjac(i,j) = zero
    end do
end do

do inod = 1, nnode
    do idim = 1, ndim
    do jdim = 1, ndim
    xjac(jdim,idim) = xjac(jdim,idim) +
    dshape(jdim,inod)*coords(idim,inod)
    end do
    end do
end do

djac = xjac(1,1)*xjac(2,2) - xjac(1,2)*xjac(2,1)
    if (djac .gt. zero) then
    ! jacobian is positive - o.k.
    xjac(1,1) = xjac(2,2)/djac
    xjac(2,2) = xjac(1,1)/djac
    xjac(1,2) = -xjac(1,2)/djac
    xjac(2,1) = -xjac(2,1)/djac
    else
    ! negative or zero jacobian
    write(7,*)'WARNING: element',jelem,'has neg.
    Jacobian'
    pnewdt = fourth
    endif

if (pnewdt .lt. pnewdtLocal) pnewdtLocal = pnewdt
    do i = 1, nnode*ndim
    bmat(i) = zero
    end do

do inod = 1, nnode
    do ider = 1, ndim
    do idim = 1, ndim
    irow = idim + (inod - 1)*ndim 

bmat(irow) = bmat(irow) +
1 xjaci(idim,ider)*dshape(ider,inod)
    end do
    end do
    end do
c
c CALCULATE INCREMENTAL STRAINS
c
do i = 1, ntens
    dstran(i) = zero
end do
!
! set deformation gradient to Identity matrix
do k1=1,3
    do k2=1,3
    defGrad(k1,k2) = zero
    end do
    defGrad(k1,k1) = one
end do
c
c COMPUTE INCREMENTAL STRAINS
c
do nodi = 1, nnode
    incr_row = (nodi - 1)*ndof
    do i = 1, ndof
    xdu(i) = du(i + incr_row,1)
    utmp(i) = u(i + incr_row)
    end do
    dNidx = bmat(1 + (nodi-1)*ndim)
    dNidy = bmat(2 + (nodi-1)*ndim)
    dstran(1) = dstran(1) + dNidx*xdu(1)
    dstran(2) = dstran(2) + dNidy*xdu(2)
    dstran(4) = dstran(4) +
1    dNidy*xdu(1) +
2    dNidx*xdu(2)
c deformation gradient 

defGrad(1,1) = defGrad(1,1) + dNidx*utmp(1)
defGrad(1,2) = defGrad(1,2) + dNidy*utmp(1)
defGrad(2,1) = defGrad(2,1) + dNidx*utmp(2)
defGrad(2,2) = defGrad(2,2) + dNidy*utmp(2)
end do

c
c CALL CONSTITUTIVE ROUTINE
c
isvinc= (kintk-1)*nsvint ! integration point increment
c
c prepare arrays for entry into material routines
c
do i = 1, nsvint
statevLocal(i)=svars(i+isvinc)
end do
c
c state variables
c
!DEC$ NOVECTOR
do k1=1,ntens
stran(k1) = statevLocal(k1)
stress(k1) = zero
end do
c
do i=1, ntens
!DEC$ NOVECTOR
do j=1, ntens
ddsdde(i,j) = zero
end do
ddsdde(i,j) = one
enddo
c
c compute characteristic element length
c
celent = sqrt(djac*dble(ninpt))
dvmat = djac*thickness
c
dvdv0 = one
call material_lib_mech(materialllib,stress,ddsdde, 

1 stran,dstran,kintk,dvdv0,dvmat,defGrad,
2 predef_loc,dpredef_loc,npredf,celent,coords_ip)
do k1=1,ntens
statevLocal(k1) = stran(k1) + dstran(k1)
end do
isvinc= (kintk-1)*nsvint ! integration point increment
update element state variables
do i = 1, nsvint
    svars(i+isvinc)=statevLocal(i)
end do
form stiffness matrix and internal force vector
dNjdx = zero
dNjdy = zero
do i = 1, ndof*nnode
    force(i) = zero
    do j = 1, ndof*nnode
    stiff(j,i) = zero
    end do
end do
dvol= wght(kintk)*djac
do nodj = 1, nnode
incr_col = (nodj - 1)*ndof
dNjdx = bmat(1+(nodj-1)*ndim)
dNjdy = bmat(2+(nodj-1)*ndim)
force_p(1) = dNjdx*stress(1) + dNjdy*stress(4)
force_p(2) = dNjdy*stress(2) + dNjdx*stress(4)
do jdof = 1, ndof
jcol = jdof + incr_col
force(jcol) = force(jcol) +
& force_p(jdof)*dvol 

end do
do nodi = 1, nnode
incr_row = (nodi -1)*ndof
dNidx = bmat(1+(nodi-1)*ndim)
dNidy = bmat(2+(nodi-1)*ndim)
stiff_p(1,1) = dNidx*ddsdde(1,1)*dNjdx
& + dNidy*ddsdde(4,4)*dNjdy
& + dNidx*ddsdde(1,4)*dNjdy
& + dNidy*ddsdde(4,1)*dNjdx
stiff_p(1,2) = dNidx*ddsdde(1,2)*dNjdy
& + dNidy*ddsdde(4,4)*dNjdx
& + dNidx*ddsdde(1,4)*dNjdx
& + dNidy*ddsdde(4,2)*dNjdy
stiff_p(2,1) = dNidy*ddsdde(2,1)*dNjdx
& + dNidx*ddsdde(4,4)*dNjdy
& + dNidy*ddsdde(2,4)*dNjdy
& + dNidx*ddsdde(4,1)*dNjdx
stiff_p(2,2) = dNidy*ddsdde(2,2)*dNjdy
& + dNidx*ddsdde(4,4)*dNjdx
& + dNidy*ddsdde(2,4)*dNjdx
& + dNidx*ddsdde(4,2)*dNjdy
do jdof = 1, ndof
icol = jdof + incr_col
do idof = 1, ndof
irow = idof + incr_row
stiff(irow,icol) = stiff(irow,icol) +
& stiff_p(idof,jdof)*dvol
end do
end do
end do
end do
end do
c
c assemble rhs and lhs
c 

do k1=1, ndof*nnode
    rhs(k1, 1) = rhs(k1, 1) - force(k1)
    do k2=1, ndof*nnode
    amatrix(k1, k2) = amatrix(k1, k2) + stiff(k1,k2)
    end do
    end do
    end do    ! end loop on material integration points
    pnewdt = pnewdtLocal
c
999 continue
c
return
end 

1.1.30 UEXPAN: User subroutine to define incremental thermal strains.

Product: Abaqus/Standard

References

• “Thermal expansion,” Section 26.1.2 of the Abaqus Analysis User’s Guide
• *EXPANSION
• “UEXPAN,” Section 4.1.16 of the Abaqus Verification Guide

Overview

User subroutine UEXPAN:

• can be used to define incremental thermal strains as functions of temperature, predefined field variables, and state variables;
• is intended for models in which the thermal strains depend on temperature and/or predefined field variables in complex ways or depend on state variables, which can be used and updated in this routine;
• is called at all integration points of elements for which the material or gasket behavior definition contains user-subroutine-defined thermal expansion; and
• is called twice per material point in each iteration during coupled temperature-displacement and coupled thermal-electrical-structural analyses.

User subroutine interface

SUBROUTINE UEXPAN(EXPAN, DEXPANDT, TEMP, TIME, DTIME, PREDEF, 1 DPRED, STATEV, CMNAME, NSTATV, NOEL)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
C
DIMENSION EXPAN(*), DEXPANDT(*), TEMP(2), TIME(2), PREDEF(*), 1 DPRED(*), STATEV(NSTATV)
user coding to define EXPAN, DEXPANDT and update
STATEV if necessary.
RETURN
END 

EXPAN(*)

Increments of thermal strain. The number of values to be defined and the order in which they are arranged depend on the type of thermal expansion being defined.

• For isotropic expansion give the isotropic thermal strain increment as the first and only component of the matrix.
• For orthotropic expansion give \Delta \epsilon _ { 1 1 } ^ { t h } , \Delta \epsilon _ { 2 2 } ^ { t h } , and \Delta \epsilon _ { 3 3 } ^ { t h } as the first, second, and third components of the matrix, respectively.
• For anisotropic expansion give \Delta \epsilon _ { 1 1 } ^ { t h } , \Delta \epsilon _ { 2 2 } ^ { t h } , \Delta \epsilon _ { 3 3 } ^ { t h } , \Delta \epsilon _ { 1 2 } ^ { t h } , \Delta \epsilon _ { 1 3 } ^ { t h } , and \Delta \epsilon _ { 2 3 } ^ { t h } . Direct components are stored first, followed by shear components in the order presented here. For plane stress only three components of the matrix are needed; give \Delta \epsilon _ { 1 1 } ^ { t h } , \Delta \epsilon _ { 2 2 } ^ { t h } , and \Delta \epsilon _ { 1 2 } ^ { t h } , as the first, second, and third components, respectively.

DEXPANDT(*)

Variation of thermal strains with respect to temperature, \partial \epsilon ^ { t h } / \partial \theta . The number of values and the order in which they are arranged depend on the type of thermal expansion being defined.

• For isotropic expansion give the variation of the isotropic thermal strain with respect to temperature as the first and only component of the matrix.
• For orthotropic expansion give \partial \epsilon _ { 1 1 } ^ { t h } / \partial \theta , \partial \epsilon _ { 2 2 } ^ { t h } / \partial \theta , and \partial \epsilon _ { 3 3 } ^ { t h } / \partial \theta as the first, second, and third components of the matrix, respectively.
• For anisotropic expansion give \partial \epsilon _ { 1 1 } ^ { t h } / \partial \theta , \partial \epsilon _ { 2 2 } ^ { t h } / \partial \theta , \partial \epsilon _ { 3 3 } ^ { t h } / \partial \theta , \partial \epsilon _ { 1 2 } ^ { t h } / \partial \theta , \partial \epsilon _ { 1 3 } ^ { t h } / \partial \theta , and \partial \epsilon _ { 2 3 } ^ { t h } / \partial \theta . Direct components are stored first, followed by shear components in the order presented here. For plane stress only three components of the matrix are needed; give \partial \epsilon _ { 1 1 } ^ { t h } / \partial \theta , \partial \epsilon _ { 2 2 } ^ { t h } / \partial \theta , and \partial \epsilon _ { 1 2 } ^ { t h } / \partial \theta , as the first, second, and third components, respectively.

Variable that can be updated

STATEV(NSTATV)

Array containing the user-defined solution-dependent state variables at this point. Except for coupled temperature-displacement and coupled thermal-electrical-structural analyses, these are supplied as values at the start of the increment and can be updated to their values at the end of the increment. For coupled temperature-displacement and coupled thermal-electrical-structural analyses, UEXPAN is called twice per material point per iteration. In the first call for a given material point and iteration, the values supplied are those at the start of the increment and can be updated. In the second call for the same material point and iteration, the values supplied are those returned from the first call, and they can be updated again to their values at the end of the increment.

User subroutine UEXPAN allows for the incremental thermal strains to be only weakly dependent on the state variables. The Jacobian terms arising from the derivatives of the thermal strains with respect to the state variables are not taken into account.