* * Example user subroutine VUINTER to model a uniform thickness interface * bonded to both surfaces. Decisions about which slave nodes are in contact * are made based on the initial configuration for the step in which the * contact pair is introduced. The interface material has "uniaxial" plasticity * in the normal direction with linear hardening. The shear behavior remains * elastic. Membrane straining of the interface does not affect the stress * transmitted across the interface. Simple heat transfer across the interface * is modeled using a conductance that is independent of gap distance or * pressure. Heat generation at the interface is not considered. * subroutine vuinter( sfd, scd, spd, svd, * stress, fluxSlv, fluxMst, sed, statev, * kStep, kInc, nFacNod, nSlvNod, nMstNod, nSurfDir, * nDir, nUSdv, nProps, NumTemp, NumExfv, numDefTfv, * jSlvUid, jMstUid, jConMstid, timStep, timGlb, * dTimCur, surfInt, surfSlv, surfMst, * rdisp, drdisp, drot, stiffDflt, condDflt, * shape, coordSlv, coordMst, alocaldir, props, * areaSlv, tempSlv, dtempSlv, exfvSlv, dexfvSlv, * tempMst, dtempMst, exfvMst, dexfvMst ) include 'vaba_param.inc' character*80 surfInt, surfSlv, surfMst dimension props(nProps), statev(nUSdv,nSlvNod), * drot(2,2,nSlvNod), sed(nSlvNod), sfd(nSlvNod), * scd(nSlvNod), spd(nSlvNod), svd(nSlvNod), * rdisp(nDir,nSlvNod), drdisp(nDir,nSlvNod), * stress(nDir,nSlvNod), fluxSlv(nSlvNod), * fluxMst(nSlvNod), areaSlv(nSlvNod), * stiffDflt(nSlvNod), condDflt(nSlvNod), * alocaldir(nDir,nDir,nSlvNod), shape(nFacNod,nSlvNod), * coordSlv(nDir,nSlvNod), coordMst(nDir,nMstNod), * jSlvUid(nSlvNod), jMstUid(nMstNod), * jConMstid(nFacNod,nSlvNod), tempSlv(nSlvNod), * dtempSlv(nSlvNod), exfvSlv(NumExfv,nSlvNod), * dexfvSlv(NumExfv,nSlvNod), tempMst(nSlvNod), * dtempMst(nSlvNod), exfvMst(NumExfv,nSlvNod), * dexfvMst(NumExfv,nSlvNod) * Indices to user-defined properties (nprops=7): parameter ( i_prp_GapCutOff = 1, * i_prp_YoungsModulus = 2, * i_prp_PoissonsRatio = 3, * i_prp_InitYield = 4, * i_prp_HardenMod = 5, * i_prp_ThickInter = 6, * i_prp_IfcCond = 7 ) * Descriptions: * i_prp_GapCutOff: cut-off init. gap dist. above which slave nodes * are not bonded. * (The rest are material properties of the interface.) * i_prp_YoungsModulus: E * i_prp_PoissonsRatio: Poisson's Ratio * i_prp_InitYield: initial yield stress * i_prp_HardenMod: hardening modulus * i_prp_ThickInter: interface thickness * i_prp_IfcCond: interface conductivity * Indices to user-defined state variables per slave node (nUSdv=3): parameter ( i_usv_CompletedInit = 1, * i_usv_BondStatus = 2, * i_usv_CurYield = 3 ) * Descriptions: * i_usv_CompletedInit: whether initializations have occured * i_usv_BondStatus: whether a slave node is bonded * i_usv_CurYield: current yield stress * Indices to the stress array: parameter ( i_str_S11 = 1, i_str_S12 = 2, i_str_S13 = 3 ) * i_str_S11: normal stress * i_str_S12: shear traction in first tangent direction * i_str_S13: shear traction in second tangent direction parameter ( zero = 0.d0, one = 1.d0, half = 0.5d0 ) * if ( statev(i_usv_CompletedInit,1) .eq. zero ) then * Note that state variables are initialized to zero by default outside * of this subroutine. Reintitialize some of them the first time VUINTER * is called for a contact pair. do kSlv = 1, nSlvNod statev(i_usv_CompletedInit,kSlv) = one gapInit = -rDisp(1,kSlv) if ( gapInit .le. props(i_prp_GapCutOff) ) then statev(i_usv_BondStatus,kSlv) = one statev(i_usv_CurYield,kSlv) = props(i_prp_InitYield) end if end do end if * Compute shear modulus. xMu = half * props(i_prp_YoungsModulus) / * (one + props(i_prp_PoissonsRatio)) if ( nDir .eq. 2 ) then * 2D: do kSlv = 1, nSlvNod if ( statev(i_usv_BondStatus,kSlv) .eq. one ) then * Compute nominal strain increments. dE11 = drDisp(1,kSlv) / props(i_prp_ThickInter) dE12 = -drDisp(2,kSlv) / props(i_prp_ThickInter) * Compute elastic trial stress. stress(i_str_S11,kSlv) = stress(i_str_S11,kSlv) * + props(i_prp_YoungsModulus)*dE11 stress(i_str_S12,kSlv) = stress(i_str_S12,kSlv) + xMu*dE12 * Determine how much to scale S11 due to any yielding. sTrial = abs(stress(i_str_S11,kSlv)) if ( sTrial .gt. statev(i_usv_CurYield,kSlv) ) then * Plastic strain increment. dEpl = ( sTrial - statev(i_usv_CurYield,kSlv) ) * / (props(i_prp_YoungsModulus)+props(i_prp_HardenMod)) * Change in yield stress. statev(i_usv_CurYield,kSlv) =statev(i_usv_CurYield,kSlv) * + props(i_prp_HardenMod) * dEpl * Corrected normal stress. stress(i_str_S11,kSlv)=sign(statev(i_usv_CurYield,kSlv), * stress(i_str_S11,kSlv)) end if * Compute heat fluxes. if( NumTemp .gt. 0 ) then if( numDefTfv .eq. nSlvNod ) then fluxSlv(kSlv) = props(i_prp_IfcCond)*areaSlv(kSlv) * * (tempMst(kSlv) - tempSlv(kSlv)) * / props(i_prp_ThickInter) else fluxSlv(kSlv) = props(i_prp_IfcCond)*areaSlv(kSlv) * * (tempMst(1) - tempSlv(kSlv)) * / props(i_prp_ThickInter) end if fluxMst(kSlv) = -fluxSlv(kSlv) end if else * Zero stress and heat flux for unbonded nodes. stress(i_str_S11,kSlv) = zero stress(i_str_S12,kSlv) = zero if( NumTemp .gt. 0 ) then fluxSlv(kSlv) = zero fluxMst(kSlv) = zero end if end if end do else * 3D: do kSlv = 1, nSlvNod if ( statev(i_usv_BondStatus,kSlv) .eq. one ) then * Compute nominal strain increments. dE11 = drDisp(1,kSlv) / props(i_prp_ThickInter) dE12 = -drDisp(2,kSlv) / props(i_prp_ThickInter) dE13 = -drDisp(3,kSlv) / props(i_prp_ThickInter) * Compute elastic trial stress. stress(i_str_S11,kSlv) = stress(i_str_S11,kSlv) * + props(i_prp_YoungsModulus)*dE11 stress(i_str_S12,kSlv) = stress(i_str_S12,kSlv) + xMu*dE12 stress(i_str_S13,kSlv) = stress(i_str_S13,kSlv) + xMu*dE13 * Determine how much to scale S11 due to any yielding. sTrial = abs(stress(i_str_S11,kSlv)) if ( sTrial .gt. statev(i_usv_CurYield,kSlv) ) then * Plastic strain increment. dEpl = ( sTrial - statev(i_usv_CurYield,kSlv) ) * / (props(i_prp_YoungsModulus)+props(i_prp_HardenMod)) * Change in yield stress. statev(i_usv_CurYield,kSlv) =statev(i_usv_CurYield,kSlv) * + props(i_prp_HardenMod) * dEpl * Corrected normal stress. stress(i_str_S11,kSlv)=sign(statev(i_usv_CurYield,kSlv), * stress(i_str_S11,kSlv) ) end if * Compute heat fluxes. if( NumTemp .gt. 0 ) then if( numDefTfv .eq. nSlvNod ) then fluxSlv(kSlv) = props(i_prp_IfcCond)*areaSlv(kSlv) * * (tempMst(kSlv) - tempSlv(kSlv)) * / props(i_prp_ThickInter) else fluxSlv(kSlv) = props(i_prp_IfcCond)*areaSlv(kSlv) * * (tempMst(1) - tempSlv(kSlv)) * / props(i_prp_ThickInter) end if fluxMst(kSlv) = -fluxSlv(kSlv) end if else * Zero stress and heat flux for unbonded nodes. stress(i_str_S11,kSlv) = zero stress(i_str_S12,kSlv) = zero stress(i_str_S13,kSlv) = zero if( NumTemp .gt. 0 ) then fluxSlv(kSlv) = zero fluxMst(kSlv) = zero end if end if end do end if return end