c user element for uncoupled springs (translation + rotational) c 2 nodes, 3-D element with dof arranged in different non-standard way: c element dof=[u1,u2,v1,v2,w1,w2,rotx1,rotx2,roty1,roty2,rotz1,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, * pt001 = 1.d-3) 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, * nCarSymm = 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), * dMassScaleFactor(nblock), * predef(nblock, nnode, npredef, nPred), adlmag(nblock) c local variables dimension symm(nCarSymm), 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_xx = props(4) dRotStiff_yy = props(5) dRotStiff_zz = props(6) eRotDamp = props(7) eDampRot = eRotDamp c user specified rotary inertia (anisotropic) ri_xx = props(9) ri_yy = props(10) ri_zz = props(11) ri_xy = props(12) ri_yz = props(13) ri_xz = props(14) 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,3,3) = am0 amass(kblock,4,4) = am0 amass(kblock,5,5) = am0 amass(kblock,6,6) = am0 c define rotary inertia amass(kblock,7,7) = ri_xx amass(kblock,9,9) = ri_yy amass(kblock,11,11) = ri_zz amass(kblock,7,9) = ri_xy amass(kblock,9,11) = ri_yz amass(kblock,7,11) = ri_xz amass(kblock,9,7) = amass(kblock,7,9) amass(kblock,11,9) = amass(kblock,9,11) amass(kblock,11,7) = amass(kblock,7,11) amass(kblock,8,8) = ri_xx amass(kblock,10,10) = ri_yy amass(kblock,12,12) = ri_zz amass(kblock,8,10) = ri_xy amass(kblock,10,12) = ri_yz amass(kblock,8,12) = ri_xz amass(kblock,10,8) = amass(kblock,8,10) amass(kblock,12,10) = amass(kblock,10,12) amass(kblock,12,8) = amass(kblock,8,12) end do else if ( lflags(iOpCode) .eq. jDtStableOnly ) then else if ( lflags(iOpCode) .eq. jIntForceOnly )then 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) = dRotStiff_xx symm(iYy) = dRotStiff_yy symm(iZz) = dRotStiff_zz symm(iXy) = zero symm(iYz) = zero symm(iZx) = zero call vuel_uti_eig33Anal( 1, symm, eigVal ) eigKMax=max(eigVal(iX),eigVal(iY), * eigVal(iZ)) c compute the smallest eigenvalue of the rotary inertia matrix eigRIMin = max(dRotStiff_xx,dRotStiff_yy,dRotStiff_zz) symm(iXx) = ri_xx symm(iYy) = ri_yy symm(iZz) = ri_zz symm(iXy) = ri_xy symm(iYz) = ri_yz symm(iZx) = ri_xz call vuel_uti_eig33Anal( 1, symm, eigVal ) eigRIMin=min(eigVal(iX),eigVal(iY), * eigVal(iZ)) 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,2) - u(kblock,1) aTraY = u(kblock,4) - u(kblock,3) aTraZ = u(kblock,6) - u(kblock,5) aVelX = v(kblock,2) - v(kblock,1) aVelY = v(kblock,4) - v(kblock,3) aVelZ = v(kblock,5) - v(kblock,5) fElasTraX = ak*aTraX fElasTraY = ak*aTraY fElasTraZ = ak*aTraZ fDampTraX = dDamp*aVelX fDampTraY = dDamp*aVelY fDampTraZ = dDamp*aVelZ forceTraX = fElasTraX + fDampTraX forceTraY = fElasTraY + fDampTraY forceTraZ = fElasTraZ + fDampTraZ c Note: this is an oversimplifed kinematics; it does not c match the kinematics of any connection type in the connector c library aRotX = u(kblock,8) - u(kblock,7) aRotY = u(kblock,10) - u(kblock,9) aRotZ = u(kblock,12) - u(kblock,11) aOmegaX = v(kblock,8) - v(kblock,7) aOmegaY = v(kblock,10) - v(kblock,9) aOmegaZ = v(kblock,12) - v(kblock,11) fElasRotX = dRotStiff_xx*aRotX fElasRotY = dRotStiff_yy*aRotY fElasRotZ = dRotStiff_zz*aRotZ fDampRotX = eRotDamp*aOmegaX fDampRotY = eRotDamp*aOmegaY fDampRotZ = eRotDamp*aOmegaZ forceRotX = fElasRotX + fDampRotX forceRotY = fElasRotY + fDampRotY forceRotZ = fElasRotZ + fDampRotZ c assemble internal load in RHS rhs(kblock,1) = -forceTraX rhs(kblock,2) = forceTraX rhs(kblock,3) = -forceTraY rhs(kblock,4) = forceTraY rhs(kblock,5) = -forceTraZ rhs(kblock,6) = forceTraZ rhs(kblock,7) = -forceRotX rhs(kblock,8) = forceRotX rhs(kblock,9) = -forceRotY rhs(kblock,10) = forceRotY rhs(kblock,11) =-forceRotZ rhs(kblock,12) = 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) aTraZOld = svars(kblock,3) fElasTraXOld = svars(kblock,4) fElasTraYOld = svars(kblock,5) fElasTraZOld = svars(kblock,6) fDampTraXOld = svars(kblock,7) fDampTraYOld = svars(kblock,8) fDampTraZOld = svars(kblock,9) aRotXOld = svars(kblock,10) aRotYOld = svars(kblock,11) aRotZOld = svars(kblock,12) fElasRotXOld = svars(kblock,13) fElasRotYOld = svars(kblock,14) fElasRotZOld = svars(kblock,15) fDampRotXOld = svars(kblock,16) fDampRotYOld = svars(kblock,17) fDampRotZOld = svars(kblock,18) energy(kblock, iElIe) = energy(kblock, iElIe) + * half*(fElasTraX+fElasTraXOld)*(aTraX-aTraXOld) + * half*(fElasTraY+fElasTraYOld)*(aTraY-aTraYOld) + * half*(fElasTraZ+fElasTraZOld)*(aTraZ-aTraZOld) + c * half*(fElasRotX+fElasRotXOld)*(aRotX-aRotXOld) + * half*(fElasRotY+fElasRotYOld)*(aRotY-aRotYOld) + * half*(fElasRotZ+fElasRotZOld)*(aRotZ-aRotZOld) energy(kblock, iElDd) = energy(kblock, iElDd) + * half*(fDampTraX+fDampTraXOld)*(aTraX - aTraXOld) + * half*(fDampTraY+fDampTraYOld)*(aTraY - aTraYOld) + * half*(fDampTraZ+fDampTraZOld)*(aTraZ - aTraZOld) + c * half*(fDampRotX+fDampRotXOld)*(aRotX - aRotXOld) + * half*(fDampRotY+fDampRotYOld)*(aRotY - aRotYOld) + * half*(fDampRotZ+fDampRotZOld)*(aRotZ - aRotZOld) c update state variables svars(kblock,1) =aTraX svars(kblock,2) =aTraY svars(kblock,3) =aTraZ svars(kblock,4) =fElasTraX svars(kblock,5) =fElasTraY svars(kblock,6) =fElasTraZ svars(kblock,7) =fDampTraX svars(kblock,8) =fDampTraY svars(kblock,9) =fDampTraZ svars(kblock,10)=aRotX svars(kblock,11)=aRotY svars(kblock,12)=aRotZ svars(kblock,13)=fElasRotX svars(kblock,14)=fElasRotY svars(kblock,15)=fElasRotZ svars(kblock,16)=fDampRotX svars(kblock,17)=fDampRotY svars(kblock,18)=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, * nCarSymm = 6 ) parameter (iX = 1, * iY = 2, * iZ = 3, * nCar = 3) * dimension eigVal(ngroup,nCar), sMat(ngroup,nCarSymm) * 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