c user element for uncoupled springs (translation + rotational) c 2 nodes, 2-D c element dof=[u1,v1,rotz1,u2,v2,rotz2] subroutine vuel( * nblock, c to be defined * rhs,amass,dtimeStable, * svars,nsvars, * energy, c * nnode,ndofel, * props,nprops, * jprops,njprops, * coords,ncrd, * u,du,v,a, * jtype,jelem, * time,period,dtimeCur,dtimePrev,kstep,kinc,lflags, * dMassScaleFactor, * predef,npredef, * jdltyp,adlmag) include 'vaba_param.inc' parameter ( zero = 0.d0, half = 0.5d0, one = 1.d0, two=2.d0 ) c operation code parameter ( jMassCalc = 1, * jIntForceAndDtStable = 2, * jExternForce = 3) c flags parameter (iProcedure = 1, * iNlgeom = 2, * iOpCode = 3, * nFlags = 3) c time parameter (iStepTime = 1, * iTotalTime = 2, * nTime = 2) c procedure flags parameter ( jDynExplicit = 17 ) c energies parameter ( iElPd = 1, * iElCd = 2, * iElIe = 3, * iElTs = 4, * iElDd = 5, * iElBv = 6, * iElDe = 7, * iElHe = 8, * iElKe = 9, * iElTh = 10, * iElDmd = 11, * iElDc = 12, * nElEnergy = 12) c predefined variables parameter ( iPredValueNew = 1, * iPredValueOld = 2, * nPred = 2) c indexing in a 3x3 symmetric rotary inertia tensor parameter (iXx = 1, * iYy = 2, * iZz = 3, * iXy = 4, * iYz = 5, * iZx = 6, * nCarTens = 6 ) c indexing in a 3-long vector parameter (iX = 1, * iY = 2, * iZ = 3, * nCar = 3) parameter (factorStable = 0.99d0) **------------------------------------ C dimension rhs(nblock,ndofel), amass(nblock,ndofel,ndofel), * dtimeStable(nblock), * svars(nblock,nsvars), energy(nblock,nElEnergy), * props(nprops), jprops(njprops), * jelem(nblock), time(nTime), lflags(nFlags), * coords(nblock,nnode,ncrd), u(nblock,ndofel), * du(nblock,ndofel), v(nblock,ndofel), a(nblock, ndofel), * predef(nblock, nnode, npredef, nPred), adlmag(nblock) c local variables dimension symm(nCarTens), eigVal(nCar) if (jtype .eq. 1001 .and. * lflags(iProcedure).eq.jDynExplicit) then dStiff = props(1) dDamp = props(2) dMass = props(3) eDampTra = dDamp c properties for rotations-diagonals of a 3x3 dRotStiff_zz = props(4) eRotDamp = props(5) eDampRot = eRotDamp c user specified rotary inertia ri_zz = props(6) if ( lflags(iOpCode).eq.jMassCalc ) then do kblock = 1, nblock c define mass matrix for translations am0 = dMass amass(kblock,1,1) = am0 amass(kblock,2,2) = am0 amass(kblock,4,4) = am0 amass(kblock,5,5) = am0 c define rotary inertia amass(kblock,3,3) = ri_zz amass(kblock,6,6) = ri_zz end do else if ( lflags(iOpCode) .eq. * jIntForceAndDtStable) then do kblock = 1, nblock amElem0 = two*dMass ak = dStiff c undamped stable time increment for translations dtTrialTransl = sqrt(amElem0/ak) c damped stable time increment critDampTransl = two*sqrt(amElem0*ak) csiTra = eDampTra/critDampTransl factDamp = sqrt(one+csiTra*csiTra) - csiTra dtTrialTransl = dtTrialTransl*factDamp*factorStable c undamped stable time increment for rotations c compute the largest eigenvalue of the elastic stiffness matrix symm(iXx) = zero symm(iYy) = zero symm(iZz) = dRotStiff_zz symm(iXy) = zero symm(iYz) = zero symm(iZx) = zero call vuel_uti_eig33Anal( 1, symm, eigVal ) eigKMax=dRotStiff_zz c compute the smallest eigenvalue of the rotary inertia matrix eigRIMin = max(dRotStiff_xx,dRotStiff_yy,dRotStiff_zz) symm(iXx) = zero symm(iYy) = zero symm(iZz) = ri_zz symm(iXy) = zero symm(iYz) = zero symm(iZx) = zero call vuel_uti_eig33Anal( 1, symm, eigVal) eigRIMin=ri_zz c undamped stable time increment for rotations dtTrialRot= sqrt(two*eigRIMin/eigKMax) c damped stable time increment critDampRot = two*sqrt(two*eigRIMin*eigKMax) csiRot = eDampRot/critDampRot factDamp = sqrt(one+csiRot*csiRot) - csiRot dtTrialRot = dtTrialRot*factDamp*factorStable c stable time increment in the element dtimeStable(kblock) = min(dtTrialTransl,dtTrialRot) c internal force and moment calcualtions aTraX = u(kblock,4) - u(kblock,1) aTraY = u(kblock,5) - u(kblock,2) aVelX = v(kblock,4) - v(kblock,1) aVelY = v(kblock,5) - v(kblock,2) fElasTraX = ak*aTraX fElasTraY = ak*aTraY fDampTraX = dDamp*aVelX fDampTraY = dDamp*aVelY forceTraX = fElasTraX + fDampTraX forceTraY = fElasTraY + fDampTraY c Note: this is an oversimplifed kinematics; it does not c match the kinematics of any connection type in the connector c library aRotZ = u(kblock,6) - u(kblock,3) aOmegaZ = v(kblock,6) - v(kblock,3) fElasRotZ = dRotStiff_zz*aRotZ fDampRotZ = eRotDamp*aOmegaZ forceRotZ = fElasRotZ + fDampRotZ c assemble internal load in RHS rhs(kblock,1) = -forceTraX rhs(kblock,2) = -forceTraY rhs(kblock,3) = -forceRotZ rhs(kblock,4) = forceTraX rhs(kblock,5) = forceTraY rhs(kblock,6) = forceRotZ c internal energy calculation c could be done with fewer state vars; used here to c prove concept aTraXOld = svars(kblock,1) aTraYOld = svars(kblock,2) fElasTraXOld = svars(kblock,3) fElasTraYOld = svars(kblock,4) fDampTraXOld = svars(kblock,5) fDampTraYOld = svars(kblock,6) aRotZOld = svars(kblock,7) fElasRotZOld = svars(kblock,8) fDampRotZOld = svars(kblock,9) energy(kblock, iElIe) = energy(kblock, iElIe) + * half*(fElasTraX+fElasTraXOld)*(aTraX-aTraXOld) + * half*(fElasTraY+fElasTraYOld)*(aTraY-aTraYOld) + c * half*(fElasRotZ+fElasRotZOld)*(aRotZ-aRotZOld) energy(kblock, iElDd) = energy(kblock, iElDd) + * half*(fDampTraX+fDampTraXOld)*(aTraX - aTraXOld) + * half*(fDampTraY+fDampTraYOld)*(aTraY - aTraYOld) + c * half*(fDampRotZ+fDampRotZOld)*(aRotZ - aRotZOld) c update state variables svars(kblock,1) =aTraX svars(kblock,2) =aTraY svars(kblock,3) =fElasTraX svars(kblock,4) =fElasTraY svars(kblock,5) =fDampTraX svars(kblock,6) =fDampTraY svars(kblock,7)=aRotZ svars(kblock,8)=fElasRotZ svars(kblock,9)=fDampRotZ end do end if end if c return end subroutine vuel_uti_eig33Anal( ngroup, sMat, eigVal ) * include 'vaba_param.inc' parameter (r_sys_Fuzz = 1.d-16) parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, * three = 3.d0, half = 0.5d0, third = one / three, * pi23 = 2.094395102393195d0, preciz = r_sys_Fuzz * 1.d4 ) c indexing in a 3x3 symmetric rotary inertia tensor parameter (iXx = 1, * iYy = 2, * iZz = 3, * iXy = 4, * iYz = 5, * iZx = 6, * nCarTens = 6 ) parameter (iX = 1, * iY = 2, * iZ = 3, * nCar = 3) * dimension eigVal(ngroup,nCar), sMat(ngroup,nCarTens) * do k = 1, ngroup s11 = sMat(k,iXx) s22 = sMat(k,iYy) s33 = sMat(k,iZz) s12 = sMat(k,iXy) s13 = sMat(k,iZx) s23 = sMat(k,iYz) fac = max(abs(s11), abs(s22), abs(s33)) facs = max(abs(s12), abs(s13), abs(s23)) xnorm = max(fac, facs) s11 = s11/xnorm s22 = s22/xnorm s33 = s33/xnorm s12 = s12/xnorm s13 = s13/xnorm s23 = s23/xnorm sh = third*(s11+s22+s33) s11 = s11 - sh s22 = s22 - sh s33 = s33 - sh fac = max(abs(s11), abs(s22), abs(s33)) facs = facs/xnorm * if( facs .lt. (preciz*fac) ) then eigVal(k,iX) = sMat(k,iXx) eigVal(k,iY) = sMat(k,iYy) eigVal(k,iZ) = sMat(k,iZz) else q = third*((s12**2+s13**2+s23**2)+half*(s11**2+s22**2+s33**2)) fac = two * sqrt(q) if( fac .gt. r_sys_Fuzz ) then ofac = two/fac else ofac = zero end if s11 = ofac*s11 s22 = ofac*s22 s33 = ofac*s33 s12 = ofac*s12 s13 = ofac*s13 s23 = ofac*s23 r = s12*s13*s23 * + half*(s11*s22*s33-s11*s23**2-s22*s13**2-s33*s12**2) if( r .ge. one-r_sys_Fuzz ) then cos1 = -half cos2 = -half cos3 = one else if( r .le. r_sys_Fuzz-one ) then cos1 = -one cos2 = half cos3 = half else ang = third * acos(r) cos1 = cos(ang) cos2 = cos(ang+pi23) cos3 =-cos1-cos2 end if eigVal(k,iX) = (sh + fac*cos1)*xnorm eigVal(k,iY) = (sh + fac*cos2)*xnorm eigVal(k,iZ) = (sh + fac*cos3)*xnorm end if end do * return end