c user element for uncoupled springs (translation + rotational)
c 2 nodes, 3-D
c element dof=[u1,v1,w1,rotx1,roty1,rotz1,u2,v2,w2,rotx2,roty2,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,
* ndload,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),
* dMassScaleFactor(nblock),
* 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_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,7,7) = am0
amass(kblock,8,8) = am0
amass(kblock,9,9) = am0
c define rotary inertia
amass(kblock,4,4) = ri_xx
amass(kblock,5,5) = ri_yy
amass(kblock,6,6) = ri_zz
amass(kblock,4,5) = ri_xy
amass(kblock,5,6) = ri_yz
amass(kblock,4,6) = ri_xz
amass(kblock,5,4) = amass(kblock,4,5)
amass(kblock,6,5) = amass(kblock,5,6)
amass(kblock,6,4) = amass(kblock,4,6)
amass(kblock,10,10) = ri_xx
amass(kblock,11,11) = ri_yy
amass(kblock,12,12) = ri_zz
amass(kblock,10,11) = ri_xy
amass(kblock,11,12) = ri_yz
amass(kblock,10,12) = ri_xz
amass(kblock,11,10) = amass(kblock,10,11)
amass(kblock,12,11) = amass(kblock,11,12)
amass(kblock,12,10) = amass(kblock,10,12)
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) = 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,7) - u(kblock,1)
aTraY = u(kblock,8) - u(kblock,2)
aTraZ = u(kblock,9) - u(kblock,3)
aVelX = v(kblock,7) - v(kblock,1)
aVelY = v(kblock,8) - v(kblock,2)
aVelZ = v(kblock,9) - v(kblock,3)
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,10) - u(kblock,4)
aRotY = u(kblock,11) - u(kblock,5)
aRotZ = u(kblock,12) - u(kblock,6)
aOmegaX = v(kblock,10) - v(kblock,4)
aOmegaY = v(kblock,11) - v(kblock,5)
aOmegaZ = v(kblock,12) - v(kblock,6)
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) = -forceTraY
rhs(kblock,3) = -forceTraZ
rhs(kblock,4) = -forceRotX
rhs(kblock,5) = -forceRotY
rhs(kblock,6) = -forceRotZ
rhs(kblock,7) = forceTraX
rhs(kblock,8) = forceTraY
rhs(kblock,9) = forceTraZ
rhs(kblock,10)= forceRotX
rhs(kblock,11)= forceRotY
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,
* 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