c user element for truss along the X global axis
c the Young's modulus is made a function of field variable
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)
C
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 3-long vector
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 VUEL subroutine for a truss element along the global X-axis only
c superimposed with rotations at both ends
c Limitations:
c It will not work if motions along another direction other than X
c Nodes must be noncoincident
c Notes:
c Define only nonzero entries; the arrays to be defined have
c been zeroed out just before this call
c for jInternForceOnly, one can use state vars to
c compute initial forces at time zero (prestress)
if (jtype .eq. 2 .and.
* lflags(iProcedure).eq.jDynExplicit) then
area0 = props(1)
eMod1 = props(2)
fv1 = props(3)
eMod2 = props(4)
fv2 = props(5)
anu = props(6)
rho = props(7)
eDampTra = zero
amassFact0 = half*area0*rho
if ( lflags(iOpCode).eq.jMassCalc ) then
do kblock = 1, nblock
c use original distance to compute mass
alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
* alenZ0*alenZ0)
am0 = amassFact0*alen0
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
end do
else if ( lflags(iOpCode) .eq.
* jIntForceAndDtStable) then
do kblock = 1, nblock
alenX0 = (coords(kblock,2,1) - coords(kblock,1,1))
alenY0 = (coords(kblock,2,2) - coords(kblock,1,2))
alenZ0 = (coords(kblock,2,3) - coords(kblock,1,3))
alen0 = sqrt(alenX0*alenX0 + alenY0*alenY0 +
* alenZ0*alenZ0)
vol0 = area0*alen0
amElem0 = two*amassFact0*alen0
alenX = alenX0
* + (u(kblock,4) - u(kblock,1))
alenY = alenY0
* + (u(kblock,5) - u(kblock,2))
alenZ = alenZ0
* + (u(kblock,6) - u(kblock,3))
alen = sqrt(alenX*alenX + alenY*alenY + alenZ*alenZ)
area = vol0/alen
c linear interpolation to determine the value of E at the
c current value of field variable. It the f.v. are outside
c the range specified, E is constant
fv_node1=predef(kblock,1,2,iPredValueNew)
fv_node2=predef(kblock,2,2,iPredValueNew)
fv=half*(fv_node1+fv_node2)
if(fv.lt.fv1) then
eMod=eMod1
else if(fv.gt.fv2) then
eMod=eMod2
else
eMod=eMod1+(eMod2-eMod1)*(fv-fv1)/(fv2-fv1)
end if
ak = area*eMod/alen
c undamped stable time increment for translations
dtTrialTransl = sqrt(amElem0/ak)
c damped stable time increment; since eDampTra=0, the
c stable time increment does not change because of damping
critDampTransl = two*sqrt(amElem0*ak)
csiTra = eDampTra/critDampTransl
factDamp = sqrt(one+csiTra*csiTra) - csiTra
dtTrialTransl = dtTrialTransl*factDamp*factorStable
dtimeStable(kblock) = dtTrialTransl
c force = E * logarithmic strain *current area
strainLog = log(alen/alen0)
fElasTra = eMod*strainLog*area
forceTra = fElasTra
c assemble internal load in RHS
rhs(kblock,1) = -forceTra
rhs(kblock,4) = forceTra
c internal energy calculation
alenOld = svars(kblock,1)
fElasTraOld = svars(kblock,2)
energy(kblock, iElIe) = energy(kblock, iElIe) +
* half*(fElasTra+fElasTraOld)*(alen - alenOld)
c update state variables
svars(kblock,1) = alen
svars(kblock,2) = fElasTra
end do
end if
end if
c
return
end