* * Example user subroutine VUINTER to model a uniform thickness interface * bonded to both surfaces. Decisions about which secondary 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, fluxSec, fluxMain, sed, statev, * kStep, kInc, nFacNod, nSecNod, nMainNod, nSurfDir, * nDir, nUSdv, nProps, NumTemp, NumExfv, numDefTfv, * jSecUid, jMainUid, jConMainid, timStep, timGlb, * dTimCur, surfInt, surfSec, surfMain, * rdisp, drdisp, drot, stiffDflt, condDflt, * shape, coordSec, coordMain, alocaldir, props, * areaSec, tempSec, dtempSec, exfvSec, dexfvSec, * tempMain, dtempMain, exfvMain, dexfvMain ) include 'vaba_param.inc' character*80 surfInt, surfSec, surfMain dimension props(nProps), statev(nUSdv,nSecNod), * drot(2,2,nSecNod), sed(nSecNod), sfd(nSecNod), * scd(nSecNod), spd(nSecNod), svd(nSecNod), * rdisp(nDir,nSecNod), drdisp(nDir,nSecNod), * stress(nDir,nSecNod), fluxSec(nSecNod), * fluxMain(nSecNod), areaSec(nSecNod), * stiffDflt(nSecNod), condDflt(nSecNod), * alocaldir(nDir,nDir,nSecNod), shape(nFacNod,nSecNod), * coordSec(nDir,nSecNod), coordMain(nDir,nMainNod), * jSecUid(nSecNod), jMainUid(nMainNod), * jConMainid(nFacNod,nSecNod), tempSec(nSecNod), * dtempSec(nSecNod), exfvSec(NumExfv,nSecNod), * dexfvSec(NumExfv,nSecNod), tempMain(nSecNod), * dtempMain(nSecNod), exfvMain(NumExfv,nSecNod), * dexfvMain(NumExfv,nSecNod) * 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 secondary 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 secondary 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 secondary 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 kSec = 1, nSecNod statev(i_usv_CompletedInit,kSec) = one gapInit = -rDisp(1,kSec) if ( gapInit .le. props(i_prp_GapCutOff) ) then statev(i_usv_BondStatus,kSec) = one statev(i_usv_CurYield,kSec) = 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 kSec = 1, nSecNod if ( statev(i_usv_BondStatus,kSec) .eq. one ) then * Compute nominal strain increments. dE11 = drDisp(1,kSec) / props(i_prp_ThickInter) dE12 = -drDisp(2,kSec) / props(i_prp_ThickInter) * Compute elastic trial stress. stress(i_str_S11,kSec) = stress(i_str_S11,kSec) * + props(i_prp_YoungsModulus)*dE11 stress(i_str_S12,kSec) = stress(i_str_S12,kSec) + xMu*dE12 * Determine how much to scale S11 due to any yielding. sTrial = abs(stress(i_str_S11,kSec)) if ( sTrial .gt. statev(i_usv_CurYield,kSec) ) then * Plastic strain increment. dEpl = ( sTrial - statev(i_usv_CurYield,kSec) ) * / (props(i_prp_YoungsModulus)+props(i_prp_HardenMod)) * Change in yield stress. statev(i_usv_CurYield,kSec) =statev(i_usv_CurYield,kSec) * + props(i_prp_HardenMod) * dEpl * Corrected normal stress. stress(i_str_S11,kSec)=sign(statev(i_usv_CurYield,kSec), * stress(i_str_S11,kSec)) end if * Compute heat fluxes. if( NumTemp .gt. 0 ) then if( numDefTfv .eq. nSecNod ) then fluxSec(kSec) = props(i_prp_IfcCond)*areaSec(kSec) * * (tempMain(kSec) - tempSec(kSec)) * / props(i_prp_ThickInter) else fluxSec(kSec) = props(i_prp_IfcCond)*areaSec(kSec) * * (tempMain(1) - tempSec(kSec)) * / props(i_prp_ThickInter) end if fluxMain(kSec) = -fluxSec(kSec) end if else * Zero stress and heat flux for unbonded nodes. stress(i_str_S11,kSec) = zero stress(i_str_S12,kSec) = zero if( NumTemp .gt. 0 ) then fluxSec(kSec) = zero fluxMain(kSec) = zero end if end if end do else * 3D: do kSec = 1, nSecNod if ( statev(i_usv_BondStatus,kSec) .eq. one ) then * Compute nominal strain increments. dE11 = drDisp(1,kSec) / props(i_prp_ThickInter) dE12 = -drDisp(2,kSec) / props(i_prp_ThickInter) dE13 = -drDisp(3,kSec) / props(i_prp_ThickInter) * Compute elastic trial stress. stress(i_str_S11,kSec) = stress(i_str_S11,kSec) * + props(i_prp_YoungsModulus)*dE11 stress(i_str_S12,kSec) = stress(i_str_S12,kSec) + xMu*dE12 stress(i_str_S13,kSec) = stress(i_str_S13,kSec) + xMu*dE13 * Determine how much to scale S11 due to any yielding. sTrial = abs(stress(i_str_S11,kSec)) if ( sTrial .gt. statev(i_usv_CurYield,kSec) ) then * Plastic strain increment. dEpl = ( sTrial - statev(i_usv_CurYield,kSec) ) * / (props(i_prp_YoungsModulus)+props(i_prp_HardenMod)) * Change in yield stress. statev(i_usv_CurYield,kSec) =statev(i_usv_CurYield,kSec) * + props(i_prp_HardenMod) * dEpl * Corrected normal stress. stress(i_str_S11,kSec)=sign(statev(i_usv_CurYield,kSec), * stress(i_str_S11,kSec) ) end if * Compute heat fluxes. if( NumTemp .gt. 0 ) then if( numDefTfv .eq. nSecNod ) then fluxSec(kSec) = props(i_prp_IfcCond)*areaSec(kSec) * * (tempMain(kSec) - tempSec(kSec)) * / props(i_prp_ThickInter) else fluxSec(kSec) = props(i_prp_IfcCond)*areaSec(kSec) * * (tempMain(1) - tempSec(kSec)) * / props(i_prp_ThickInter) end if fluxMain(kSec) = -fluxSec(kSec) end if else * Zero stress and heat flux for unbonded nodes. stress(i_str_S11,kSec) = zero stress(i_str_S12,kSec) = zero stress(i_str_S13,kSec) = zero if( NumTemp .gt. 0 ) then fluxSec(kSec) = zero fluxMain(kSec) = zero end if end if end do end if return end