SDT-nlsim         Contents     Functions              PDF Index

1.4  Hyper-visco 3D model

• Kinematics
• Invariants, strain rates
• Stress, equations of motion, time integration at a material point
• Hyperelastic potential, polyconvexity
• Compression or volume change behavior
• Deviatoric laws
• Mixed U-P elements small deformation
• Cross check examples
• Mixed U-P elements large deformation
• Viscoelastic implementations

The SDT/OpenFEM implementation of volumes also partially supports hyper viscoelastic laws under large deformation. This section derived from []. This first section summarizes classical results of continuum mechanics, which are needed for later discussion of constitutive laws exhibiting hyperelasticity, viscoelasticity and rate independent hysteresis under large deformation [2, 3, 4].

1.4.1  Kinematics

Considering a three dimensional euclidean space with base e₁,e₂,e₃, a body Ω may be described by a region from this space. Each point P of Ω is called a material point and its coordinates are described by

in the reference configuration of a Lagrangian description of the problem. The body position may evolve with time, and time derivatives give velocity v and acceleration a as

The kinematic description is based on the knowledge of the initial and current positions of any material point. The difference between initial and current positions of material points is the displacement vector {u}. Material behavior is typically described as a function of strains. SDT uses the simplest strain measure, the displacement gradient, given by

where δ_ij denotes the Kronecker delta function. This tensor links the reference and the current configurations, making it a hybrid deformation description. It may be divided in two separate parts: rotation R and dilatation U in a polar decomposition

where R is orthonormal, and U positive definite and symmetric. To obtain rotation independent deformation tensors, quadratic forms of the deformation tensor are used. For a Lagrangian description, the right Cauchy-Green tensor

will be used (the left Cauchy-Green tensor is B=F F^T). The Green Lagrange strain e_ij=1/2(C_ij−δ_ij) is a common alternative. The main difference between the two is that C tends to the identity for nil deformation, while e tends to zero.

1.4.2  Invariants, strain rates

For isotropic materials, properties should not depend of orientation and behavior should not dependent on rotation R. Invariants of the strain tensors are often used in constitutive formulation to achieve rotation independence. The three classical invariants are

where one should note that I₃=det(C)=det²(F) quantifies the square of the volume change at the material point, and J=I₃^1/2 the volume change. The derivatives of invariants with respect to deformations will also be necessary and are classically given by

and

where 1_ijkl=1/2(δ_ikδ_jl+δ_ilδ_jk). For highly incompressible materials, it is usual to describe the potential associated with the deviatoric part of the strain using reduced invariants, defined by

It is useful to remind that for zero deformation F=I, I_i=331 so that it makes sense to use expressions involving Ī₁−3 (as in the Yeoh model). To check values on also reminds that dI1dC=I,dI2dC=2I, dI3dc=I

xxx

Deformation rate tensors are also often needed. The velocity gradient tensor is commonly used and given by

As this gradient is not rotation independent, it is handy to decompose it into strain rate and spin tensors, by taking its symmetric and the skew-symmetric parts,

The strain rate tensor may be used in viscoelastic laws while the spin tensor is usually discussed in finite deformation plasticity models [5].

1.4.3  Stress, equations of motion, time integration at a material point

When computing work of deformable bodies the dual quantities (energy conjugate pair) are stress and strain, as force and displacement are dual for rigid body mechanics.

In the Lagrangian reference frame, the second Piola-Kirchhoff S=J F⁻¹ σ F^−T tensor is used for stress (degenerates to stresses in the original configuration for small deformation but large rotation). The hybrid reference tensor (initial frame for strains to current frame for stresses) is called first Piola-Kirchhoff Π=Jσ F^−T (or Boussinesq, or also nominal stress) and is not symmetric. In linearized conditions, the Eulerian description given by the Cauchy stress σ is used.

Stress and strain energy present different equivalent expressions, forming dual pairs

Material or constitutive models are laws that describe the evolution of stresses for an history of deformations and possibly internal states. They are system models of what happens at a material point.

Given a displacement gradient tensor, stress at timestep n+1 must be an explicit function of the gradient displacement tensor and its internal states at the last timestep n (considering explicit implementation),

The development of material laws described by (1.65) should typically be tested and tuned separately from the global structural application.

With material laws defined, equations of motion at the structure level are derived from the fact that the work of all forces for any kinematically acceptable virtual displacement field w is equal to zero. In other words,

Implicit solution of this equation typically require a tangent stiffness obtained from integration of the tangent material stiffness (1.69).

1.4.4  Hyperelastic potential, polyconvexity

A material is said to be hyperelastic, if its stress only depends on the current strain. It is classically shown that this is verified if stress is obtained as the derivative of a potential ψ that depends on deformation tensors. Derivation from a potential is a sufficient condition to demonstrate that deformations are completely reversible, despite being nonlinear. The non-linearities involved are

• geometric: large deformations is typically involved requiring non-linear strain descriptions
• material: as the operating range can be large, there is no physical reason for materials to respond linearly

For an isotropic material, behavior must be independent of rotation. This can be used to demonstrate that the potential may always be expressed as a function of strain invariants. The stress is thus found as

where the convention of summing repeated indexes is used. Note that the potential is decomposed in deviatoric (section 1.4.6) and compression (section 1.4.5) behavior so that the vector actually used is dWdI=∂ ψ^D/∂ I₁∂ ψ^D/∂ I₂∂ ψ^D/∂ I₃∂ ψ^I/∂ I₃ .

Other consequence of the expression of the potential in terms of the invariants is that it may always be written in the form [3]

with β_i functions of the invariants that may be determined according to the model. These three β functions, therefore can be seen as a non parametric manner of modeling hyperelasticity. Since going through the whole space spanned by I₁, I₂ and I₃ is unpractical, most common approaches are to use either selected order representations as polynomials on the invariants or order independent models that are not based on invariants, as will be described in section 1.4.6.

The only restriction on the energy potential is its polyconvexity [2, 6] which ensures stability. For scalar stress or 0D, the potential is the integral of the force/displacement curve. This function can be a polynomial or any convex function corresponding to the fact that the instantaneous stiffness or slope of the force/displacement curve is positive. In 3D, ensuring polyconvexity for a tensor is not straightforward but corresponds to the positive definite nature of the material stiffness. For an isotropic material, this material stiffness can be expressed as [7, 8]

leading to the expression of the structural geometric stiffness [9]

SDT uses naming conventions dWdI,d2WdI2 for ∂ψ/∂ I_n, ∂² ψ/∂ I_n∂ I_m and expressions of these values for classic materials are given in section 1.4.5 for compression or isotropic volume change and section 1.4.6 for deviatoric or shear behavior.

1.4.5  Compression or volume change behavior

Materials used for bushings are typically incompressible. Compressibility is described by the isotropic part of strains and a specific discussion of constitutive laws is needed. Remaining behavior is associated with deviatoric strains and discussed separately.

To achieve this separation for hyperelastic models, one may divide the potential in a deviatoric and an isotropic part ψ=ψ^D+ψ^I, which means that the same decomposition applies to the stresses, aside from numerical issues. This implies the introduction of two separate models that are superposed in the end.

As the third invariant is the only one that has no deviatoric components, compression models are normally based on J=I₃^1/2. The isotropic potential for linear models normally used is

This potential however allows full compression or volume suppression. For this reason, the Ciarlet-Geymonat [10] potential was developed with a logarithmic penalty term,

The stress tensor resulting from the isotropic part may be computed as σ_p=pδ_ij in an Eulerian reference frame, which may be transferred to a Lagrangian reference frame by S_p=F⁻¹(σ_p)F^−T.

1.4.6  Deviatoric laws

Three main groups of hyperelastic potentials are typically [11] distinguished: invariant, principal strains and physic based potentials.

Invariant based formulations were developed first and are still widely used, because of their simpler formulation and low numerical cost.

The Yeoh model uses polynomials of (Ī₁−3)=(J^−2/3 I₁−3)=(I₃^−1/3 I₁−3), (Ī₂−3)=(J^−2/3 I₂−3)=(I₃^−1/3 I₂−3), and J−1=(I₃^1/2) whose values vanish for zero displacement. And is given by

used in (1.67) for stress computations and (1.69) for Jacobian computations. Note that Abaqus uses a generalization combining polynomials of the form ∑_i,j=1^N C_ij(Ī₁−3)ⁱ(Ī₂−3)^j.

The neo-hookean model only considers C₁₀ (asymptotic elasticity gives G=2C₁₀). The Mooney-Rivlin model adds the second invariant and constant C₀₁ to this formulation and is given by

where asymptotic elasticity gives G=2(C₁₀+C₀₁).

Stability requires ∂² ψ^D/∂ C² to be a positive definite matrix, but the above potentials do not verify stability for some strain. A very simple potential featuring a very good fitting unconditional stability [12, 13] is the Carroll model given by

(1.75)

As the proof for the stability of invariant based models is not straightforward, the use of potentials associated with principal strains λ_i (equivalently, the eigenvalues of the matrix C) ensures the stability, as the matrix C is symmetric (note that these correspond to the square of the singular values of F). The Ogden model [11] is the most used model of this type

This open form ensures that the model is able to fit any curve, but it comes at the cost of a large number of parameters and the need to compute eigenvalues at each step.

So called physic based potentials are based on assumptions of how the elastomer behaves at a microscopic level. The main origin of this type of model is the chain model [11], which assumes that the elasticity of networked chains is due to the entropic changes, thus it may be determined by the number of possible states. This chain model gave origin to the widely used Arruda-Boyce model, which assumes that the chains are distributed in a form that 8 chains are attached to the center and to the vertices of a cube. More recently the assumption that the chain is restricted to a tube gave origin to the tube model, and finally the assumption that these chains are continually distributed in a sphere gave origin to the sphere models which are very good at capturing the behavior with few parameters.

m_hyper implements basic checks of material properties as shown in the examples below. This allows verification of how κ controls how close ν is to the fully incompressible value 0.5.

m_hyper('urn','PadC{2.5264,0,0,120,3,rho1n,tyMoon,unTM}');
m_hyper('urn','PadA{2.5264,-0.9177,0.4711,1200,3,f .35,g .5688,rho1n,tyYeoh,unTM}')

1.4.7  Mixed U-P elements small deformation

For nearly incompressible materials one uses mixed elements that interpolate displacement and pressure using different shape functions (see [14] section 12.5, [15] or Abaqus theory manual Section 3.2.3 Hybrid incompressible solid element formulation). These are called mixed U-P elements. In SDT switch to this formulation is done for rule 20002 (hexa20 for u, hexa8 for p/J with 8 Gauss point, similarly tetra10/tetra4 with 4 Gauss point, penta15/penta6 xxx, see) For small deformation, the mixed formulation uses a generalized strain of the form

The constitutive law is then split in deviatoric and pressure/isotropic/isochore parts using the expression of the deviatoric part of the strain {є_dev} = [D₀ − 2/3 {m} {m}^T] {є}, where {m}^T=111000 and D₀=diag222111, leading to

the virtual work associated with s_p is assumed to be equal to zero, so that one has the weak set of equations 0=∫_ΩN_p (p/K−{m}^T {є})= [K_qp^T]{q}+K_pp{p} (which is avoids numerical problems by considering a p interpolation that is coarser than the displacement interpolation). When doing transients, the equations of motion would be of the form

where the augmentation by pressure DOF is not a dynamic equation but a set of constraints.

As a reminder, in the isotropic elasticity case, one has K=E/(3(1−2ν)), G=E/(2(1+ν)) and

For cross checking of large deformation, one considers the zero deformation case. One then has I₁=I₂=3, I₃=1, thus ∂ ψ^I/∂ I₃ is vanishing and ∂² ψ^I/∂ I₃²=κ/4, which applying (1.69) is consistent with (1.80). For the Yeoh model one has the isochore part

and the deviatoric contribution

leading to 2C₁₀=G (as also said in [16] table 2.2).

1.4.8  Cross check examples

Incompressible cube with edge length L is considered as a first example. The static problem is

The incompressibility condition yields є_yy=є_zz=−є_xx/2. Assuming linear displacements, boundary conditions leads to the system

Hence, we deduce that

Nearly compressible cube is a second example. The static problem is

The constitutive law can be rewritten as

Assuming linear displacements, the deformation reads

so that

Boundary conditions lead to the following system

Recalling that

the system becomes

So, we deduce

Note that E=3G for incompressible materials so that this solution tends to the one obtained in the incompressible case.

1.4.9  Mixed U-P elements large deformation

In the large deformation case, one considers a development similar to the two field case described in [15]. The generalized strain combines the displacement gradient and the internal state g related to volume change/pressure. The second Piola-Kirchhoff stress is given by

and g related to volume change by considering the weak relation

where the potential is given by (1.71) where ∂ ψ^I∂ I₃=κ (1−I₃^−1/2) or (1.72) where ∂ ψ^I∂ I₃=κ2 (I₃^−1/2−I₃⁻¹ ). As a result the generalized stress components are obtained using

For time integration, considering a bulk viscosity xxx.

1.4.10  Viscoelastic implementations

Since [17, 18] demonstrated that the stress relaxation with saturation was the most appropriate, this is the currently implemented model on the derivatives of the potential dWdI=∂ ψ^D/∂ I₁∂ ψ^D/∂ I₂∂ψ^D/∂ I₃∂ ψ^I/∂ I₃ . The law (1.26) is used on the deviatoric components of ∂ψ^D/∂ I_i.

For the compression part ∂ ψ^I/∂ I₃, which is much stiffer for many elastomers, a viscous model ∂ψ^I = ∂ψ₀^I+(κ_v/κ)∂ψ₀^I is used leading to the increment equation

note that the have a viscosity dominate compression behavior at a target frequency f_v in Hz you should use k_v=κ/(2π f_v).