KOKO\Mimi
15 KiB
if (cmname(1:4) .eq. 'MAT1') then
call VUANISOHYPER_STRAIN1(argument_list)
else if (cmname(1:4) .eq. 'MAT2') then
call VUANISOHYPER_STRAIN2(argument_list)
end if
VUANISOHYPER_STRAIN1 and VUANISOHYPER_STRAIN2 are the actual subroutines containing the anisotropic hyperelastic models for each material MAT1 and MAT2, respectively. Subroutine VUANISOHYPER_STRAIN merely acts as a directory here. The argument list can be the same as that used in subroutine VUANISOHYPER_STRAIN. The material names must be in uppercase characters since cmname is passed in as an uppercase character string.
Example: Orthotropic Saint-Venant Kirchhoff model
As a simple example of the coding of subroutine VUANISOHYPER_STRAIN, consider the generalization to anisotropic hyperelasticity of the Saint-Venant Kirchhoff model. The strain energy function of the Saint-Venant Kirchhoff model can be expressed as a quadratic function of the Green strain tensor, \varepsilon ^ { G } , as
U (\boldsymbol {\varepsilon} ^ {G}) = \frac {1}{2} \boldsymbol {\varepsilon} ^ {G}: \mathbf {D}: \boldsymbol {\varepsilon} ^ {G},
where is the fourth-order elasticity tensor. The derivatives of the strain energy function with respect to the Green strain are given as
\frac {\partial U}{\partial \varepsilon^ {G}} = \mathbf {D}: \varepsilon^ {G},
\frac {\partial^ {2} U}{\partial \varepsilon^ {G} \partial \varepsilon^ {G}} = \mathbf {D}.
However, subroutine VUANISOHYPER_STRAIN must return the derivatives of the strain energy function with respect to the modified Green strain tensor, \overline { { \varepsilon } } ^ { G } , and the volume ratio, J, which can be accomplished easily using the following relationship between \varepsilon ^ { G } , \overline { { \varepsilon } } ^ { G } , and :
\boldsymbol {\varepsilon} ^ {G} = J ^ {\frac {2}{3}} \overline {{\boldsymbol {\varepsilon}}} ^ {G} + \frac {1}{2} (J ^ {\frac {2}{3}} - 1) \mathbf {I},
where is the second-order identity tensor. Thus, using the chain rule we find
\frac {\partial U}{\partial \overline {{\varepsilon}} ^ {G}} = J ^ {\frac {2}{3}} \frac {\partial U}{\partial \varepsilon^ {G}},
\frac {\partial U}{\partial J} = \frac {\partial \varepsilon^ {G}}{\partial J}: \frac {\partial U}{\partial \varepsilon^ {G}},
\frac {\partial^ {2} U}{\partial \overline {{\varepsilon}} ^ {G} \partial \overline {{\varepsilon}} ^ {G}} = J ^ {\frac {4}{3}} \frac {\partial^ {2} U}{\partial \varepsilon^ {G} \partial \varepsilon^ {G}},
\frac {\partial^ {2} U}{\partial J ^ {2}} = \frac {\partial^ {2} \varepsilon^ {G}}{\partial J ^ {2}}: \frac {\partial U}{\partial \varepsilon^ {G}} + \frac {\partial \varepsilon^ {G}}{\partial J}: \frac {\partial^ {2} U}{\partial \varepsilon^ {G} \partial \varepsilon^ {G}}: \frac {\partial \varepsilon^ {G}}{\partial J},
\frac {\partial^ {2} U}{\partial \overline {{\varepsilon}} ^ {G} \partial J} = \frac {2}{3 J} J ^ {\frac {2}{3}} \frac {\partial U}{\partial \varepsilon^ {G}} + J ^ {\frac {2}{3}} \frac {\partial^ {2} U}{\partial \varepsilon^ {G} \partial \varepsilon^ {G}}: \frac {\partial \varepsilon^ {G}}{\partial J},
where
\frac {\partial \pmb {\varepsilon} ^ {G}}{\partial J} = \frac {2}{3 J} J ^ {\frac {2}{3}} (\overline {{\pmb {\varepsilon}}} ^ {G} + \frac {1}{2} \mathbf {I}) = \frac {2}{3 J} (\pmb {\varepsilon} ^ {G} + \frac {1}{2} \mathbf {I})
and
\frac {\partial^ {2} \varepsilon^ {G}}{\partial J ^ {2}} = - \frac {1}{3 J} \frac {\partial \varepsilon^ {G}}{\partial J}.
In this example an auxiliary function is used to facilitate indexing into a fourth-order symmetric tensor. The subroutine would be coded as follows:
subroutine vuanisohyper_strain (
C Read only -
* nblock,
* jElem, kIntPt, kLayer, kSecPt,
* cmname,
* ndir, nshr, nstatev, nfieldv, nprops,
* props, tempOld, tempNew, fieldOld, fieldNew,
* stateOld, ebar, detu,
C Write only -
* uDev, duDe, duDj,
* d2uDeDe, d2uDjDj, d2uDeDj,
* stateNew )
C
include 'vaba_param.inc'
C
dimension props(nprops),
* tempOld(nblock),
* fieldOld(nblock,nfieldv),
* stateOld(nblock,nstatev),
* tempNew(nblock),
* fieldNew(nblock,nfieldv),
* ebar(nblock,ndir+nshr), detu(nblock),
* uDev(nblock), duDe(nblock,ndir+nshr), duDj(nblock),
* d2uDeDe(nblock, *), d2uDjDj(nblock),
* d2uDeDj(nblock,ndir+nshr),
* stateNew(nblock,nstatev)
character*80 cmname
parameter( half = 0.5d0, one = 1.d0, two = 2.d0,
* third = 1.d0/3.d0, twoths = 2.d0/3.d0, four = 4.d0,
* dinv = 0.d0 )
Orthotropic Saint-Venant Kirchhoff strain energy function
(3D)
D1111 = props(1)
D1122 = props(2)
D2222 = props(3)
D1133 = props(4)
D2233 = props(5)
D3333 = props(6)
D1212 = props(7)
D1313 = props(8)
D2323 = props(9)
do k = 1, nblock
d2UdE11dE11 = D1111
d2UdE11dE22 = D1122
d2UdE11dE33 = D1133
d2UdE22dE11 = d2UdE11dE22
d2UdE22dE22 = D2222
d2UdE22dE33 = D2233
d2UdE33dE11 = d2UdE11dE33
d2UdE33dE22 = d2UdE22dE33
d2UdE33dE33 = D3333
d2UdE12dE12 = D1212
d2UdE13dE13 = D1313
d2UdE23dE23 = D2323
xpow = exp ( log(detu(k)) * twoths )
detuInv = one / detu(k)
C
C
E11 = xpow * ebar(k,1) + half * ( xpow - one )
E22 = xpow * ebar(k,2) + half * ( xpow - one )
E33 = xpow * ebar(k,3) + half * ( xpow - one )
E12 = xpow * ebar(k,4)
E23 = xpow * ebar(k,5)
E13 = xpow * ebar(k,6)
term1 = twothds * detuInv
dE11Dj = term1 * (E11 + half)
dE22Dj = term1 * (E22 + half)
dE33Dj = term1 * (E33 + half)
dE12Dj = term1 * E12
dE23Dj = term1 * E23
dE13Dj = term1 * E13
term2 = - third * detuInv
d2E11DjDj = term2 * dE11Dj
d2E22DjDj = term2 * dE22Dj
d2E33DjDj = term2 * dE33Dj
d2E12DjDj = term2 * dE12Dj
d2E23DjDj = term2 * dE23Dj
d2E13DjDj = term2 * dE13Dj
C
dUdE11 = d2UdE11dE11 * E11
* + d2UdE11dE22 * E22
* + d2UdE11dE33 * E33
dUdE22 = d2UdE22dE11 * E11
* + d2UdE22dE22 * E22
* + d2UdE22dE33 * E33
dUdE33 = d2UdE33dE11 * E11
* + d2UdE33dE22 * E22
* + d2UdE33dE33 * E33
dUdE12 = two * d2UdE12dE12 * E12
dUdE23 = two * d2UdE23dE23 * E23
dUdE13 = two * d2UdE13dE13 * E13
U = half * ( E11*dUdE11 + E22*dU
* + E12*dUdE12 + E13*dUdE13 +
uDev(k) = U
duDe(k,1) = xpow * dUdE11
duDe(k,2) = xpow * dUdE22
duDe(k,3) = xpow * dUdE33
C
C
duDe(k,4) = xpow * dUdE12
duDe(k,5) = xpow * dUdE23
duDe(k,6) = xpow * dUdE13
C
xpow2 = xpow * xpow
C Only update nonzero components
d2uDeDe(k,indx(1,1)) = xpow2 * d2UdE11dE11
d2uDeDe(k,indx(1,2)) = xpow2 * d2UdE11dE22
d2uDeDe(k,indx(2,2)) = xpow2 * d2UdE22dE22
d2uDeDe(k,indx(1,3)) = xpow2 * d2UdE11dE33
d2uDeDe(k,indx(2,3)) = xpow2 * d2UdE22dE33
d2uDeDe(k,indx(3,3)) = xpow2 * d2UdE33dE33
d2uDeDe(k,indx(4,4)) = xpow2 * d2UdE12dE12
d2uDeDe(k,indx(5,5)) = xpow2 * d2UdE23dE23
d2uDeDe(k,indx(6,6)) = xpow2 * d2UdE13dE13
C
duDj(k) = dUdE11*dE11Dj + dUdE22*dE22Dj + dUdE33*dE33Dj
* + two * (dUdE12*dE12Dj + dUdE13*dE13Dj
* + dUdE23*dE23Dj)
d2uDjDj(k) = dUdE11*d2E11DjDj + dUdE22*d2E22DjDj
* + dUdE33*d2E33DjDj
* + two * (dUdE12*d2E12DjDj + dUdE13*d2E13DjDj
* + dUdE23*d2E23DjDj)
* + d2UdE11dE11 * dE11Dj * dE11Dj
* + d2UdE22dE22 * dE22Dj * dE22Dj
* + d2UdE33dE33 * dE33Dj * dE33Dj
* + two * (d2UdE11dE22 * dE11Dj * dE22Dj
* + d2UdE11dE33 * dE11Dj * dE33Dj
* + d2UdE22dE33 * dE22Dj * dE33Dj)
* + four * (d2UdE12dE12 * dE12Dj * dE12Dj
* d2UdE13dE13 * dE13Dj * dE13Dj
* d2UdE23dE23 * dE23Dj * dE23Dj)
C
d2uDeDj(k,1) = xpow * (term1 * dUdE11
* + d2UdE11dE11 * dE11Dj
* + d2UdE11dE22 * dE22Dj
* + d2UdE11dE33 * dE33Dj)
d2uDeDj(k,2) = xpow * (term1 * dUdE22
* + d2UdE22dE11 * dE11Dj
* + d2UdE22dE22 * dE22Dj
* + d2UdE22dE33 * dE33Dj)
d2uDeDj(k,3) = xpow * (term1 * dUdE33
* + d2UdE33dE11 * dE11Dj
* + d2UdE33dE22 * dE22Dj
* + d2UdE33dE33 * dE33Dj )
d2uDeDj(k,4) = xpow * ( term1 * dUdE12
* + two * d2UdE12dE12 * dE12Dj )
d2uDeDj(k,5) = xpow * ( term1 * dUdE23
* + two * d2UdE23dE23 * dE23Dj )
d2uDeDj(k,6) = xpow * ( term1 * dUdE13
* + two * d2UdE13dE13 * dE13Dj )
end do
C
return
end
C
integer function index( i, j )
C
include 'vaba_param.inc'
C
C
Function to map index from Square to Triangular storage
C of symmetric matrix
C
ii = min(i,j)
jj = max(i,j)
C
index = ii + jj*(jj-1)/2
C
return
end
1.2.12 VUCHARLENGTH: User subroutine to define characteristic element length at a material point.
Product: Abaqus/Explicit
References
• *CHARACTERISTIC LENGTH
• “VUCHARLENGTH,” Section 4.1.32 of the Abaqus Verification Guide
Overview
User subroutine VUCHARLENGTH:
• is called at all material points of elements for which the material definition includes a user-defined characteristic element length and the constitutive model requires a characteristic length;
• allows the definition of characteristic element length at a material point as a function of element topology, nodal and material point coordinates, and material orientation;
• can use field variables that are passed in; and
• can use solution-dependent state variables that are passed in.
Defining characteristic element length
The characteristic element length defined in user subroutine VUCHARLENGTH is used by Abaqus in regularization schemes needed to mitigate mesh dependency in constitutive models that include strain-softening, such as damage models (“Damage evolution and element removal for ductile metals,” Section 24.2.3 of the Abaqus Analysis User’s Guide), and concrete (“Concrete smeared cracking,” Section 23.6.1 of the Abaqus Analysis User’s Guide). It could be used with built-in Abaqus material models as well as user subroutine–based material models.
The characteristic element length coming in user subroutine VUCHARLENGTH has a default value based on the geometric mean. This default value is a typical length of a line across an element for a first-order element and is half of the same typical length for a second-order element. For trusses the default value is a characteristic length along the element axis. For membranes and shells the default value is a characteristic length in the reference surface. For axisymmetric elements the default value is a characteristic length in the r–z plane only.
Inside user subroutine VUCHARLENGTH you can redefine the value of the characteristic element length based on the element topology and geometry. The characteristic element length defined in user subroutine VUCHARLENGTH at a particular material point is passed to other user subroutines that are called at the same material point, such as user subroutines VFABRIC, VUMAT, VUSDFLD, and VUEOS. The characteristic element length calculated in user subroutine VUCHARLENGTH at a particular material point is also used in built-in Abaqus material models that require characteristic length and are called at the same material point.
Array of element type and geometric properties
jElType contains information about the element type and the geometry. jElType(1) provides information about the shape of the element.
| jElType(1) | Shape |
| 1 | line |
| 2 | triangle |
| 3 | quadrilateral |
| 4 | tetrahedron |
| 5 | wedge |
| 6 | hexahedron |
jElType(2) provides information about the element’s dimensionality.
| jElType (2) | Space |
| 1 | 2D and plane strain |
| 2 | 3D |
| 3 | axisymmetric |
| 4 | plane stress |
jElType(3) provides information about the section of the element.
| jElType(3) | Section |
| 1 | solid |
| 2 | shell |
| 3 | truss |
| 4 | membrane |
Elements
User subroutine VUCHARLENGTH can be used with membrane elements; shell elements; truss elements; and plane stress, plane strain, axisymmetric, and three-dimensional solid elements.
Special consideration for 8-node continuum shell elements
For 8-node hexahedron continuum shell elements (SC8R and SC8RT), the order of nodes in the array of nodal coordinates (coordNode) passed to user subroutine VUCHARLENGTH depends on the stacking
direction. For a single element nodes 1–4 correspond to the bottom face and nodes 5–6 correspond to the top face.
User subroutine interface
subroutine vucharlength(
c Read only variables-
1 nblock, nfieldv, nprops, ncomp, ndim, nnode, nstatev,
2 kSecPt, kLayer, kIntPt, jElType, jElem,
3 totalTime, stepTime, dt,
4 cmname, coordMp, coordNode, direct, T, props,
5 field, stateOld,
c Write only variables-
6 charLength )
c
include 'vaba_param.inc'
c
dimension jElType(3), jElem(nblock), coordMp(nblock,ndim),
1 coordNode(nblock, nnode, ndim),
2 direct(nblock,3,3), T(nblock,3,3), props(nprops),
3 stateOld(nblock, nstatev), charLength(nblock, ncomp),
4 field(nblock, nfieldv)
c
character*80 cmname
c
do 100 k = 1, nblock
user coding to define charLength(nblock, ncomp)
100 continue
c
return
end
Variable to be defined
charLength(nblock,ncomp)
Characteristic element length.
Variables passed in for information
nblock
Number of material points to be processed in this call to user subroutine VUCHARLENGTH.
nfieldv
Number of user-defined external field variables.
nprops
User-specified number of user-defined material properties.
ncomp
User-specified number of components of characteristic element length. If ncomp is greater than 1, only the first component of characteristic element length would be used in the built-in Abaqus material models. However, all the components could be used in the above-mentioned user subroutines.
ndim
Number of coordinate directions: 2 for two-dimensional models and 3 for three-dimensional models.
nnode
Number of nodes of the element.
nstatev
Number of user-defined state variables that are associated with this material type (define this as described in “Allocating space” in “User subroutines: overview,” Section 18.1.1 of the Abaqus Analysis User’s Guide).
kSecPt
Section point number within the current layer.
kLayer
Layer number (for composite shells).
kIntPt
Integration point number.
jElType(3)
Array containing information about the element type and geometry.
jElem(nblock)
Array of element numbers.
totalTime
Value of total time. The time at the beginning of the step is given by totalTime-stepTime.
stepTime
Value of time since the step began.
dt
Time increment size.