```txt 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 ``` ```txt 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 ``` ```fortran 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 ``` ```fortran 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, ``` ```txt 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 ``` ```txt 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 ``` ```fortran 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 ```txt 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.