Q. I was wondering if it might be possible to change the code so that more than
48 material parameters can be read into LS-DYNA.
A. There is no direct way of getting all your parameters into the *MAT_USER card.
An indirect way is to define a load curve with your parameter values, e.g.
*DEFINE_CURVE
10
1., parameter 1
2., parameter 2
...
Then you can access this load curve from within your material routine.
Q. How do we get volumetric strain?
A. The volumetric strain may be obtained from the variable dfe
(=V/V0=1+eps_{kk}) in common block aux18loc, and is also computed in usrsld.
Q. Is there a common block that contains the following information?
1) user-specified size of cm array (lmc)
2) user-specified size of hsv array (nhv)
3) flag that identifies an analysis as implicit static or dynamic analysis
A. Regarding the quantities lmc and nhv (did you mean nhsv?) - if you are looking
for the parameters that are entered on the *SECTION cards, then these are
available from the cmusr array as
lmc=nint(cmusr(6,mxe))
nhsv=nint(cmusr(7,mxe))
e.g. see subroutines usrshl and usrsld in dyn21b.F. The cm array that you
refer to is in fact in the cmusr array for element properties, not material
properties - see the code samples in dyn21b.F.
Note that lmc cannot be greater than 40, and the maximum value for nhsv
depends on NHISVUE, defined in nhisparm.inc.
However, if you are looking for the corresponding parameters for the materials
properties, rather than the element properties, then the material properties
are indeed available in the cm array, and the size of that is hard-coded to be
48. If you want to input more material properties, the best way is to put them
in a curve definition, to which the material subroutines and refer and read
off the properties. The material history variables are limited by NHISVAR,
also defined in nhisparm.inc.
The implicit flag is in ascntl(11):
c ascntl(11) = IMASS, flag for implicit time integration
c .lt.0.0: time driven dynamic, load curve
c 0.0: static
c 1.0: dynamic
c 2.0: dynamic analysis using modal superposition
where ascntl is defined in implicit1.inc
Q. I have a question about parameters for coordinates of center of mass for rigid parts
in user subroutines. A customer would like to use them with user defined friction control.
Could you please let me know how they can be treated in usrfrc()?
A. The rigid body coordinates should be available in the array a(n54),
where n54 is defined in bk07.inc. If we call this array rbcor, then its
dimensions are
rbcor(3,*)
and the center of mass of part n (internal numbering) are in
xmc = rbcor(1,n)
ymc = rbcor(2,n)
zmc = rbcor(3,n)
Q. We need to get the element formulation. For years we are using a "line"
that was advised to get that info in the usermat. Here it is:
mxe=ia(lc1s+5*(nnm1+lft-1))
and then
elform=ia(n1+nmmat+mxe-1)
(n1 is from BK07 common, nmat from BK00, lc1s from BK13).
This does not work for solids.
A. lc1s is for shell elements. For solid elements, try
mxe = ia(lc1h+9*(nnm1+lft-1))
Qf: What does lc1* stand for what does it contain?
Af: These contain the element information, i.e. element number and nodal
connectivities. The suffixes are h: hex (or solids), s: shells, b: beams, t:
thick shells.
Q. How do we access the initial coordinates of nodes in an user defined
material subroutine?
A. I presume you are familiar with how to access current coordinates, element
nodes etc. The initial coordinates are stored in the array r_mem(dm_rots),
accessed as real*8 x0(3,*).
Q. Do you have an example for a discrete beam umat?
A. The file dbeam_umat.tgz has an example.
Q. Is it possible to create user materials routines for SPH?
A. You can write a user material for etype.eq.'solid' and that should be
usable for SPH.
Q. What is the bulk modulus and shear modulus for anisotropic materials?
A. The bulk modulus is typically calculated as max(C11,C22,C33) for
anisotropic materials, and the shear modulus as max(C44,C55,C66).
Q. How do we change the bulk modulus from within umat, and should I change it
if there is damage?
You can change the sound speed modulus from with your material. If you look,
e.g. in subroutine urmathn, you will see that is sets variables ss(lft) in
common block soundloc and cb(lft) in common block aux35loc to the P-wave
modulus used for calculating the sound speed. You can modify these as
necessary at each time step by putting these common blocks in your umat
routine. Note that while ss(lft) only is assigned, i.e. only for the first
element in the group lft-llt, the variable cb(i) is assigned individually for
each element. It is the cb values that are used to calculate the time step for
each element, so you can assign the appropriate value to cb in your routine.
Whether or not you want to change the time step based on the damaged moduli is
a decision that you have to make. Some materials in LS-DYNA do it while some
don't.
Q. Is there a way to gain access to the coordinates of the integration points
in the user material subroutine?
A. The shape function values are available in common block bk32:
common/bk32/h(8,9)
where the first index is the node number and the second is the integration point number.
You can add this common block to your subroutine and then combine these with
the nodal coordinate data to compute the coordinates of the integration
points. Please see one of the user element routines, e.g. usrsld for how to
obtain the nodal coordinate data.
The integration point number is available as parameter ipt in the call to
usrmat, and is passed to both urmathn and urmats - you can change the
parameter list for your material to include this information. It will take on
values 1-8 for a fully-integrated element.
Q. How do we locate the coordinates for integration points in shells that have
integration points along the thickness of the element?
A. For fully-integrated shells, e.g. type 16, which has 4 integration points
in the plane, the loops run as follows:
ipt = 0
! in-plane loop
do np=1,4
! through-thickness
do ipt_thk=1,nip
ipt=ipt+1
call material_routine
So ipt is the index to all the integration points, and ipt is the
through-thickness integration point, and both are available in urmats. For
reduced-integration elements with one point in the plane, ipt=ipt_thk.
The in-plane location of the points are at the standard 2x2 Gauss locations,
and the through-thickness locations are available in tcoor in urmats, and the
standard locations are given in the quadrature point tables under
SECTION_SHELL in the manual.
Q. It appears tcoor gives the thickness location for the parent element
e.g. 0,+-0.577 for 3 (Gauss) integration points through the thickness . Is
there a way to get the thickness or coordinates of the integration points
through the thickness of the actual element.
A. To convert the Gauss isoparametric locations into physical locations, you
need two pieces of information: (1) the thickness of the shell, and (2) where
the mid-plane of the shell lies relative to the nodal plane, i.e. the NLOC
value.
To get the thickness see the question "How can I get shell thickness from
within a umat routine?" in
http://ftp.lstc.com/anonymous/outgoing/jday/faq/user_defined_materials.faq
Getting the NLOC value is a little trickier in my experience. I think the most
reliable way would be to enter it as another parameter in the user material -
I understand it would mean repeating the same information in two places, but
perhaps you could make it a parameter in the input file if you have a non-zero
NLOC.
Q. a. Inside a UMAT subroutine, how do I access to the integration points at
different layers of a multi-layer shell element? I need to update the stress
differently on each layer.
b. In case of a fully-integrated shell, how could I figure out the location
of 4 integration points on a given layer of the shell with respect to the 4
shell element nodes? The reason I need to know this is that in the
composite unit cell I am looking at, different integration points will have
different local elastic properties and I need to capture that.
A. In the usrmat routine, you have:
nip: total number of integration points
(in-plane and through-thickness)
ipt: index of current integration point, from 1 through nip
ipt_thk: index of current through-thickness integration point for
fully-integrated elements (4 in-plane integration points)
=ipt for single-point integration
tcoor: normalized coordinate of through-thickness integration point
For fully-integrated elements, you should be able to figure out the index of
the current in-plane integration point from ipt and ipt_thk.
The in-plane Gauss points for the fully-integrated elements are arranged in
the standard way as (-,-), (+,-), (+,+), (-,+), so you should be able to
compute the coordinates of the Gauss points as necessary.
Q. With setting ICOMP = 1 under *SECTION_SHELL, I have defined different
angles for each integration point of the shell element: B1 B2 B3 B4 B5 B6 B7
B8 in the input simulation file. My question now is how can I access B1 B2 B3
B4 B5 B6 B7 B8 angles inside the UMAT subroutine?
A. You can get them from a(lc2+nint(cm(48))), where they are stored in order,
and lc2 is in bk13.inc
If these will be used to transform stresses etc., then you should keep in
mind that (a) the transformations should consider the material angles with
respect to the current element system, and (b) LS-DYNA has already calculated
them for you. If you look in urmats, hsvs(i,ind_q1) and hsvs(i,ind_q2) have
the cosine and negative sine of the relevant angle, and these are used to
transform the stresses and strains before you get to the umat routines.
Q. Can I get the area of a beam element?
A. For the type 2 Belytschko beam and for type 3 truss element, the
cross-section area is accessed in the array a(nb05) [nb05 is in
bk05.inc] as follows -
refer to a(nb05) as fibl(5,*)
then area is fibl(1,nnm1+i), where i is from lft to llt.
Q. Where can I get the following information:
1) Implicit step/increment number
2) Implicit iteration number
3) Implicit analysis type (value of IMASS in *CONTROL_IMPLICIT_DYNAMICS)
A. These are available in
1. Implicit time-step counter is in isolvr(104) - the array is given in
implicit1.inc
2. Non-linear iteration number is in scalar variable itr, available in common
block imloci:
common /imloci/ l20,ite,itr
3. Analysis type is available in ascntl(11), with the array ascntl being
accessible through implicit1.inc
ascntl(11) = IMASS, flag for implicit time integration
.lt.0.0: time driven dynamic, load curve
.eq.0.0: static
1.0: dynamic
2.0: dynamic analysis using modal superposition
Q. This involves the use of beam/truss elements during an implicit analysis
with nonlinear materials (namely the user materials). Regardless of how I
configure the test problem or the properties of the mateiral, LS-DYNA is
calling the user material with a strain increment value equal to 2.0.
A. It seems that the user material subroutine for the truss was not returning
the appropriate stiffness entry for the main element routine to calculate the
element stiffness.
To resolve this, insert the following line in the subroutine urmatt in
dyn21.F:
do 10 i=lft,llt
einc(i)=d1(i)*sig1(i)
cb(i)=ss(lft)
wxxdt(i)=ym ! ---- insert this line
10 continue
The variable labelled wxxdt here is being used to hold the material stiffness
for calculating the element stiffness in the main routine.
Q. How do I define a tangent stiffness for user-defined cohesive elements?
A. Dev r102784 (10/29/2015) onwards allows defining the tangent stiffness in
the umatXXc routines.
For older versions, here's how you can do it:
1. Include the following two common blocks in your subroutine along with their
threadprivate/TASKCOMMON declarations:
common/vect8loc/dsave(nlq,6,6)
common/bki03iloc/lnodim(nlq,48),ndofpn,nnpke,melemt,imlft,imllt,
& is17loc,is18loc,imp_mxe
c$omp threadprivate (/vect8loc/)
c$omp threadprivate (/bki03iloc/)
You will find these in subroutine usrsld in dyn21b.f and in other places.
2. In your code, check for is17loc.eq.1, and if true, populate the dsave array
in vect8loc with the directional stiffnesses on the diagonal and with zeros in
the off-diagonal elements.
if (is17loc.eq.1) then
do i=lft,llt
dsave(i,1,1)=stifr(i)
dsave(i,2,1)=0.
dsave(i,3,1)=0.
dsave(i,1,2)=0.
dsave(i,2,2)=stifs(i)
dsave(i,3,2)=0.
dsave(i,1,3)=0.
dsave(i,2,3)=0.
dsave(i,3,3)=stift(i)
enddo
endif
where stifr, stifs, and stift are the stiffnesses in the two in-plane
directions r and s and the normal direction t, respectively, and have been
defined appropriately.
Q. How do I access history variables in cohesive user materials?
A. The aux array contains the history variables. This array is accessed as
do i=lft,llt
aux(i,k) = (k-th history variable)
enddo
This aux array starts at the 6th position in the larger history variables
sequence, after the six stresses. This is in contrast to the history variables
for regular solid elements, where they start in the 7th position after the six
stresses and the equivalent plastic strain.
Q. How do I arrange the history variables when using both a user-defined
material and a user-defined EOS?
A. If you use an user material with an user EOS, you should set the NHV for
the material to be 4 larger than what you actually need, and start your
material history variables (in hsv or hsvs) only from the 5th position. The
first 4 positions are reserved for EOS variables such as internal energy, bulk
viscosity and current volume.
The hist array passed to the user EOS routines is located after the material
history, so you can use that directly.