User Subroutine Interface
SUBROUTINE VUEL(nblock,rhs,amass,dtimeStable,svars,nsvars,
1 energy,
2 nnode,ndofel,props,nprops,jprops,njprops,
3 coords,mcrd,u,du,v,a,
4 jtype,jElem,
5 time,period,dtimeCur,dtimePrev,kstep,kinc,
6 lflags,
7 dMassScaleFactor,
8 predef,npredef,
9 jdltyp, adlmag)
C
include 'vaba_param.inc'
C operational code keys
parameter ( jMassCalc = 1,
* jIntForceAndDtStable = 2,
* jExternForce = 3)
C flag indices
parameter (iProcedure = 1,
* iNlgeom = 2,
* iOpCode = 3,
* nFlags = 3)
C energy array indices
parameter ( iElPd = 1,
* iElCd = 2,
* iElIe = 3,
* iElTs = 4,
* iElDd = 5,
* iElBv = 6,
* iElDe = 7,
* iElHe = 8,
* iUnused = 9,
* iElTh = 10,
* iElDmd = 11,
* iElDc = 12,
* nElEnergy = 12)
C predefined variables indices
parameter ( iPredValueNew = 1,
* iPredValueOld = 2,
* nPred = 2)
C time indices
parameter (iStepTime = 1,
* iTotalTime = 2,
* nTime = 2)
dimension rhs(nblock,ndofel),amass(nblock,ndofel,ndofel),
1 dtimeStable(nblock),
2 svars(nblock,nsvars),energy(nblock,nElEnergy),
3 props(nprops),jprops(njprops),
4 jElem(nblock),time(nTime),lflags(nFlags),
5 coords(nblock,nnode,mcrd),
6 u(nblock,ndofel), du(nblock,ndofel),
7 v(nblock,ndofel), a(nblock, ndofel),
8 dMassScaleFactor(nblock),
9 predef(nblock,nnode,npredef,nPred),
* adlmag(nblock)
do kblock = 1,nblock
user coding to define rhs, amass, dtimeStable, svars and energy
end do
RETURN
END
Variables to Be Defined
Some of
the following arrays depend on the value of the
lflags array.
- rhs
-
An array containing the contributions of each element to the right-hand-side
vector of the overall system of equations. Depending on the settings of the
lflags array, it contains either the internal
force from the element or the external load calculated from the specified
distributed loads.
- amass
-
An array containing the contribution of each element to the mass matrix of
the overall system of equations.
All nonzero entries in amass should be
defined. Moreover, the mass matrix must be symmetric. There are several other
requirements that apply depending on the active degrees of freedom specified.
These requirements are explained in detail below.
- dtimeStable
-
A scalar value defining, for each element, the upper limit of the time
increment for stability considerations. This would be the maximum time
increment to be used in the subsequent increment for this element to be stable
(to satisfy the Courant condition). This value depends strongly on the element
formulation, and it is important that is computed appropriately.
- svars
-
An array containing the values of the solution-dependent state variables
associated with each element. The number of such variables is
nsvars (see below). You define the meaning of
these variables.
This array is passed into
VUEL containing the values of these variables at the start of
the current increment. In most cases they should be updated to be the values at
the end of the increment. In rare cases such an update is not required.
- energy
-
The array energy contains the values of the energy quantities
associated with each element. The values in this array when VUEL is called are the element
energy quantities at the start of the current increment. They should be updated to the
correct values at the end of the current increment; otherwise, plots of the energy
balance for the entire model will not be accurate. You cannot modify or update the
kinetic energy. Depending on the element formulation, many of these energies could be
zero at all times. The entries in the array are as follows:
| energy(nblock,iElPd
)
|
Plastic dissipation.
|
| energy(nblock,iElCd)
|
Creep dissipation.
|
| energy(nblock,iElIe)
|
Internal energy.
|
| energy(nblock,iElTs)
|
Transverse shear energy.
|
| energy(nblock,iElDd)
|
Material damping dissipation.
|
| energy(nblock,iElBv)
|
Bulk viscosity dissipation.
|
| energy(nblock,iElDe)
|
Drill energy.
|
| energy(nblock,iElHe)
|
Hourglass energy.
|
| energy(nblock,iElTh)
|
Heat energy.
|
| energy(nblock,iElDmd)
|
Damage dissipation.
|
| energy(nblock,iElDc )
|
Distortion control energy.
|
Variables Passed in for Information
- Arrays:
- props
-
A floating point array containing the nprops
real property values defined for use with each element processed.
nprops is the user-specified number of real
property values. See
Defining the Element Properties.
- jprops
-
An integer array containing the njprops
integer property values defined for use with each element processed.
njprops is the user-specified number of integer
property values. See
Defining the Element Properties.
- coords
-
An array containing the original coordinates of the nodes of the element.
coords(kblock,k1,k2) is the
k2th coordinate of the
k1th node of the
kblock element.
- u, du, v,
a
-
Arrays containing the basic solution variables (displacements, rotations,
temperatures, pressures, depending on the degree of freedom) at the nodes of
the element. Values are provided as follows, illustrated below for the
k1th degree of freedom of the
kblock element:
| u(kblock,k1)
|
Total value of the variables (such as
displacements or rotations) at the end of the current increment.
|
| du(kblock,k1)
|
Incremental values of the variables in
the current increment.
|
| v(kblock,k1)
|
Time rate of change of the variables
(velocities, rates of rotation) at the midpoint of the increment.
|
| a(kblock,k1)
|
Accelerations of the variables at the
end of the current increment.
|
- jElem
-
jElem(kblock) contains the element number for
the kblock element.
- adlmag
-
adlmag(kblock) is the total load magnitude of
the load type jdltyp (integer identifying the
load number for distributed load type Un) distributed load at the end of the current increment for
distributed loads of type Un.
- predef
-
An array containing the values of predefined field variables, such as
temperature in an uncoupled stress/displacement analysis, at the nodes of the
element (Predefined Fields).
The second index, k2, indicates the local
node number on the kblock element. The third
index, k3, indicates the variable: the
temperature is stored if the index is 1, and the predefined field variables are
stored if the indices are greater than or equal to 2. The fourth index of the
array, k4, is either 1 or 2, with 1 indicating
the value of the field variable at the end of the increment and 2 indicating
the value of the field variable at the beginning of the increment.
| predef(kblock,k2,1,k4)
|
Temperature.
|
| predef(kblock,k2,2,k4)
|
First predefined field variable.
|
| predef(kblock,k2,3,k4)
|
Second predefined field variable.
|
| Etc.
|
Any other predefined field variable.
|
| predef(kblock,k2,k3,k4)
|
Value of the
(k3–1)th predefined field variable at the
k2th node of the element at the beginning or the
end of the increment.
|
| predef(kblock,k2,k3,1)
|
Values of the variables at the end of
the current increment.
|
| predef(kblock,k2,k3,2)
|
Values of the variables at the
beginning of the current increment.
|
- lflags
-
An array containing the flags that define the current solution procedure and
requirements for element calculations.
| lflags(iProcedure)
|
Defines the procedure type. See
Results File
for the key used for each procedure.
|
| lflags(iNlgeom)=0
|
Small-displacement analysis.
|
| lflags(iNlgeom)=1
|
Large-displacement analysis (nonlinear
geometric effects included in the step; see
General and Perturbation Procedures).
|
| lflags(iOpCode)=jMassCalc
|
Define the mass matrix
amass in the beginning of the analysis.
|
| lflags(iOpCode)=jIntForceAndDtStable
|
Define the element internal force.
Define the stable time increment as well.
|
| lflags(iOpCode)=jExternForce
|
Define the distributed load effect on the external force
associated with the element.
|
- dMassScaleFactor
-
An array containing the mass scale factors for each element.
- time(iStepTime)
-
Current value of step time.
- time(iTotalTime)
-
Current value of total time.
- Scalar parameters:
- nblock
-
Number of user elements to be processed in this call to
VUEL.
- dtimeCur
-
Current time increment.
- dtimePrev
-
Previous time increment.
- period
-
Time period of the current step.
- ndofel
-
Number of degrees of freedom in the elements processed.
- nsvars
-
User-defined number of solution-dependent state variables associated with
the element (Defining the Number of Solution-Dependent Variables That Must Be Stored within the Element).
- nprops
-
User-defined number of real property values associated with the elements
processed (Defining the Element Properties).
- njprops
-
User-defined number of integer property values associated with the elements
processed (Defining the Element Properties).
- mcrd
-
mcrd is defined as the maximum of the
user-defined maximum number of coordinates needed at any node point (Defining the Maximum Number of Coordinates Needed at Any Nodal Point)
and the value of the largest active degree of freedom of the user element that
is less than or equal to 3. For example, if you specify that the maximum number
of coordinates is 1 and the active degrees of freedom of the user element are
2, 3, and 6 mcrd will be 3. If you specify that
the maximum number of coordinates is 2 and the active degree of freedom of the
user element is 11, mcrd will be 2.
- nnode
-
User-defined number of nodes on the elements (Defining the Number of Nodes Associated with the Element).
- jtype
-
Integer defining the element type. This is the user-defined integer value
n in element type VUn (Assigning an Element Type Key to a User-Defined Element).
- kstep
-
Current step number.
- kinc
-
Current increment number.
- npredef
-
Number of predefined field variables, including temperature. For user
elements
Abaqus/Explicit
uses one value for each field variable per node.
VUEL Conventions
The solution variables (displacement, velocity, etc.) are arranged on a
node/degree of freedom basis. The degrees of freedom of the first node are
first, followed by the degrees of freedom of the second node, etc. The degrees
of freedom that will be updated automatically in
Abaqus/Explicit
are: 1–3 (displacements), 4–6 (rotations), 8 (pressure), and 11 (temperature).
Depending on the procedure type (see below), only some of the degrees of
freedom listed above will be updated. Other degrees of freedom will not be
updated by the time integration procedure in
Abaqus/Explicit
and, hence, should not be used.
The mass matrix defined in user subroutine
VUEL must be symmetric. In addition, the following requirements
apply:
You must be using appropriate lumping techniques to provide a mass matrix
that follows these requirements. For the rotational degrees of freedom at a
particular node in three-dimensional analyses, you can specify a fully
populated symmetric 3 × 3 inertia tensor.
Usage with General Nonlinear Procedures
The following illustrates the use in explicit dynamic procedures:
Direct-Integration Explicit Dynamic Analysis (lflags(iProcedure)=17)
-
Automatic updates for degrees of freedom 1–6, 8, and 11.
-
The governing equations are as described in
Explicit Dynamic Analysis.
-
Coding for the operational code
lflags(iOpCode)=jExternForceis optional.
Transient Fully Coupled Thermal-Stress Analysis (lflags(iProcedure)=74)
Example: Structural User Element
A structural user element has been created to demonstrate the usage of
subroutine
VUEL. These user-defined elements are applied in a number of
analyses. The following excerpt is from the verification problem that invokes
the structural user element in an explicit dynamic procedure:
USER ELEMENT, NODES=2, TYPE=VU1, PROPERTIES=4, COORDINATES=3,
VARIABLES=12
1, 2, 3
ELEMENT, TYPE=VU1
101, 101, 102
ELGEN, ELSET=VUTRUSS
101, 5
UEL PROPERTY, ELSET=VUTRUSS
0.002, 2.1E11, 0.3, 7200.
The user element consists of two nodes that are assumed to lie parallel to
the x-axis. The element behaves similarly to a linear
truss element. The supplied element properties are the cross-sectional area,
Young's modulus, Poisson's ratio, and density, respectively.
The next excerpt shows the listing of the subroutine. The user subroutine
has been coded for use in an explicit dynamic analysis. The names of the
verification input files associated with the subroutine and these procedures
can be found in
VUEL.
subroutine vuel(
* nblock,
* rhs,amass,dtimeStable,
* svars,nsvars,
* energy,
* 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'
c operation code
parameter ( jMassCalc = 1,
* jIntForceAndDtStable = 4)
c flags
parameter (iProcedure = 1,
* iNlgeom = 2,
* iOpCode = 3,
* nFlags = 3)
c procedure flags
parameter ( jDynExplicit = 17 )
c time
parameter (iStepTime = 1,
* iTotalTime = 2,
* nTime = 2)
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)
parameter (factorStable = 0.99d0)
parameter (zero = 0.d0, half = 0.5d0, one = 1.d0, two=2.d0)
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), l(nFlags),
* coords(nblock,nnode,ncrd), u(nblock,ndofel),
* du(nblock,ndofel), v(nblock,ndofel), a(nblock, ndofel),
* predef(nblock, nnode, npredef, nPred), adlmag(nblock),
* dMassScaleFactor(nblock)
c Notes:
c Define only nonzero entries; the arrays to be defined
c have been zeroed out just before this call
if (jtype .eq. 1001 .and.
* lflags(iProcedure).eq.jDynExplicit) then
area0 = props(1)
eMod = props(2)
anu = props(3)
rho = props(4)
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
ak = area*eMod/alen
c stable time increment for translations
dtimeStable(kblock) = factorStable*sqrt(amElem0/ak)
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
|