General finite-strain deformation has the following characteristics:

The updated Lagrangian method is applied to simulate geometric nonlinearities (accessed with NLGEOM,ON). Assuming all variables, such as coordinates x_i, displacements u_i, strains ε_ij, stresses σ_ij, velocities v_i, volume V and other material variables have been solved for and are known at time t; one solves for a set of linearized simultaneous equations having displacements (and hydrostatic pressures in the mixed u-P formulation) as primary unknowns to obtain the solution at time t + Δt. These simultaneous equations are derived from the element formulations which are based on the principle of virtual work:

(3–89)

where:

σij = Cauchy stress component

ui = displacement

xi = current coordinate

V = volume of deformed body

S = surface of deformed body on which tractions are prescribed

The internal virtual work can be indicated by:

(3–90)

where:

Element formulations are obtained by differentiating this virtual work expression (Bonet and Wood([237]) and Gadala and Wang([292])). In derivation, only linear differential terms are kept and all higher order terms are ignored so that finally a linear set of equations can be obtained.

In element formulation, material constitutive law is used to create the relation between stress increment and strain increment. The constitutive law only reflects the stress increment due to straining. However, the Cauchy stress is affected by the rigid body rotation and is not objective (not frame invariant). An objective stress is needed, therefore, to be able to be applied in constitutive law. One of these is Jaumann rate of Cauchy stress expressed by McMeeking and Rice([293])

(3–91)

where:

= Jaumann rate of Cauchy stress

= time rate of Cauchy stress

Therefore, the Cauchy stress rate is:

(3–92)

Using the constitutive law, the stress change due to straining can be expressed as:

(3–93)

where:

cijkl = material constitutive tensor

vi = velocity

The Cauchy stress rate can be written as:

(3–94)

Pure displacement formulation only takes displacements or velocities as primary unknown variables. All other quantities such as strains, stresses and state variables in history-dependent material models are derived from the displacements. It is the most widely used formulation and is able to handle most nonlinear deformation problems.

The differentiation of δW:

(3–95)

From Equation 3–94, the stress differentiation can be derived as:

(3–96)

where:

The differentiation of dV is:

(3–97)

where:

Substitution of Equation 3–96 and Equation 3–97 into Equation 3–95 yields:

(3–98)

The third term is unsymmetric and is usually insignificant in most of deformation cases. Hence, it is ignored. The final pure displacement formulation is:

(3–99)

The above equation is a set of linear equations in terms of Du_i or displacement change, which can be solved by standard linear solvers. This formulation is exactly the same as the one published by McMeeking and Rice([293]). The stiffness has two terms: the first one is material stiffness due to straining; the second one is stiffness due to geometric nonlinearity (stress stiffness).

Since no other assumption is made on the size or type of deformation, the formulation can be applied to any deformation problems (small deformation, finite deformation, small deformation with large rotation, stress stiffening, etc.) so it is called a general element formulation.

To achieve higher efficiency, the second term or stress stiffness is included only if requested for analyses with geometric nonlinearities (NLGEOM,ON, or PSTRES,ON), buckling analysis (ANTYPE,BUCKLE), or linear perturbation analysis (PERTURB).

The above pure displacement formulation is computationally efficient. However, the accuracy of any displacement formulation is dependent on Poisson's ratio or the bulk modulus. In such formulations, volumetric strain is determined from derivatives of displacements, which are not as accurately predicted as the displacements themselves. Under nearly incompressible conditions (Poisson's ratio is close to 0.5 or bulk modulus approaches infinity), any small error in the predicted volumetric strain will appear as a large error in the hydrostatic pressure and subsequently in the stresses. This error will, in turn, also affect the displacement prediction since external loads are balanced by the stresses, and may result in displacements very much smaller than they should be for a given mesh--this is called locking-- or, in some cases, in no convergence at all.

Another disadvantage of pure displacement formulation is that it is not to be able to handle fully incompressible deformation, such as fully incompressible hyperelastic materials.

To overcome these difficulties, mixed u-P formulations were developed. In these u-P formulations of the current-technology elements, the hydrostatic pressure or volume change rate is interpolated on the element level and solved on the global level independently in the same way as displacements. The final stiffness matrix has the format of:

(3–100)

where:

Δu = displacement increment

= hydrostatic pressure increment

Since hydrostatic pressure is obtained on a global level instead of being calculated from volumetric strain, the solution accuracy is independent of Poisson's ratio and bulk modulus. Hence, it is more robust for nearly incompressible material. For fully incompressible material, mixed u-P formulation has to be employed in order to get solutions.

The pressure degrees of freedom are brought to global level by using internal or external nodes. The internal nodes are different from the regular (external) nodes in the following aspects:

The interpolation function of pressure is determined according to the order of elements. To remedy the locking problem, they are one order less than the interpolation function of strains or stresses. For most current-technology elements, the hydrostatic pressure degrees of freedom are introduced by the internal nodes. The number of pressure degrees of freedom, number of internal nodes, and interpolation functions are shown in Table 3.1: Interpolation Functions of Hydrostatic Pressure of Current-Technology Elements.

Table 3.1: Interpolation Functions of Hydrostatic Pressure of Current-Technology Elements

Element	KEYOPT(6)	Internal nodes		Functions
PLANE182 selective reduced integration and uniform reduced integration	1	1	1
PLANE182 Enhanced strain formulation	1	2	3
PLANE183
SOLID185 selective reduced integration and uniform reduced integration	1	1	1
SOLID185 Enhanced strain formulation	1	2	4
SOLID186 Uniform reduced integration and full integration
SOLID187	1	1	1
SOLID187	2	2	4
SOLSH190	1	2	4
SOLID272	1	KEYOPT(2) / 3	KEYOPT(2)	on r-z plane and Fourier interpolation in the circumferential (θ) direction
SOLID273	1	KEYOPT(2)	KEYOPT(2) x 3	on r-z plane and Fourier interpolation in the circumferential (θ) direction

In Table 3.1: Interpolation Functions of Hydrostatic Pressure of Current-Technology Elements, , , , , and are the pressure degrees of freedom at internal node i. s, t, and r are the natural coordinates.

For SOLID285, the hydrostatic pressure degrees of freedom are introduced by extra degrees of freedom (HDSP) at each node. The total number of pressures and interpolation function of hydrostatic pressure are shown in Table 3.2: Interpolation Functions of Hydrostatic Pressure for SOLID285.

Table 3.2: Interpolation Functions of Hydrostatic Pressure for SOLID285

Element		Functions
285	4

, , , and are the pressure degrees of freedom at each element node i. s, t, and r are the natural coordinates.

This formulation is for nearly incompressible materials other than hyperelastic materials. For these materials, the volumetric constraint equations or volumetric compatibility can be defined as (see Bathe([2]) for details):

(3–101)

where:

K = bulk modulus

P can also be defined as:

(3–102)

In mixed formulation, stress is updated and reported by:

(3–103)

where:

so that the internal virtual work Equation 3–90 can be expressed as:

(3–104)

Introduce the constraint Equation 3–101 by Lagrangian multiplier , the augmented internal virtual work is:

(3–105)

Substitute Equation 3–103 into above; it is obtained:

(3–106)

where:

Take differentiation of Equation 3–105, ignore all higher terms of Du_i and D than linear term, the final formulation can be expressed as:

(3–107)

This is a linear set of equations of Du_i and D (displacement and hydrostatic pressure changes). In the final mixed u-P formulation, the third term is the stress stiffness and is included only if requested (NLGEOM,ON or PSTRES,ON). The rest of the terms are based on the material stiffness. The first term is from material constitutive law directly or from straining; the second term is because of the stress modification (Equation 3–103); the fourth and fifth terms are the extra rows and columns in stiffness matrix due to the introduction of the extra degree of freedom: pressure, that is, K_uP, K_Pu and K_PP as in Equation 3–100.

The stress stiffness in the above formulation is the same as the one in pure displacement formulation. All other terms exist even for small deformation and are the same as the one derived by Bathe([2]) for small deformation problems.

It is worthwhile to indicate that in the mixed formulation of the higher order elements (PLANE183 , SOLID186 and SOLID187 with KEYOPT(6) = 1), elastic strain only relates to the stress in the element on an averaged basis, rather than pointwise. The reason is that the stress is updated by Equation 3–103 and pressure is interpolated independently in an element with a function which is one order lower than the function for volumetric strain. For lower order elements (PLANE182, SOLID185), this problem is eliminated since either B-bar technology or uniform reduced integration is used; volumetric strain is constant within an element, which is consistent with the constant pressure interpolation functions (see Table 3.1: Interpolation Functions of Hydrostatic Pressure of Current-Technology Elements). In addition, this problem will not arise in element SOLID187 with linear interpolation function of (KEYOPT(6) = 2). This is because the order of interpolation function of is the same as the one for volumetric strain. In other words, the number of degrees of freedom in one element is large enough to make consistent with the volumetric strain at each integration point. Therefore, when mixed formulation of element SOLID187 is used with nearly incompressible material, the linear interpolation function of or KEYOPT(6) = 2 is recommended.

A special formulation is necessary for fully incompressible hyperelastic material because the volume constraint equation is different and hydrostatic pressure can not be obtained from material constitutive law. Instead, it must be calculated separately. For these types of materials, the stress is updated by:

(3–108)

where:

= deviatoric component of Cauchy stress tensor

The deviatoric component of deformation tensor defined by the e_ij term of Equation 3–89 can be expressed as:

(3–109)

The internal virtual work (Equation 3–90) can be shown using and :

(3–110)

The volume constraint is the incompressible condition. For a fully incompressible hyperelastic material, it can be as defined by Sussman and Bathe([125]), Bonet and Wood([237]), Crisfield([294]):

(3–111)

where:

= determinant of deformation gradient tensor

Xi = original coordinate

Vo = original volume

As in the mixed u-P formulation I (u-P Formulation I), the constraint Equation 3–111 was introduced to the internal virtual work by the Lagrangian multiplier . Then, differentiating the augmented internal virtual work, the final formulation is obtained.

This formulation is similar to the formulation for nearly incompressible materials, i.e. Equation 3–107. The only major difference is that [K_PP] = [0] in this formulation. This is because material in this formulation is fully incompressible.

When material behavior is almost incompressible, the pure displacement formulation may be applicable. The bulk modulus of material, however, is usually very large and thus often results in an ill-conditioned matrix. To avoid this problem, a special mixed u-P formulation is therefore introduced. The almost incompressible material usually has small volume changes at all material integration points. A new variable is introduced to quantify this small volume change, and the constraint equation

(3–112)

is enforced by the Lagrangian multiplier method, similar to the u-P Formulation II as shown by Equation 3–110 and Equation 3–111.

The final set of linear equations of mixed formulations (see Equation 3–100) can be grouped into two:

(3–113)

(3–114)

Equation 3–113 are the equilibrium equations and Equation 3–114 are the volumetric constraint equations. The total number of active equilibrium equations on a global level (indicated by N_d) is the total number of displacement degrees of freedom without any prescribed displacement boundary condition. The total number of volumetric constraint equations (indicated by N_p) is the total number of pressure degrees of freedom in all mixed u-P elements. The optimal ratio of N_d/N_p is 2 for 2D elements and 3 for 3D elements. When N_d/N_p is too small, the system may have too many constraint equations which may result in a severe locking problem. On the other hand, when N_d/N_p is too large, the system may have too few constraint equations which may result in too much deformation and loss of accuracy.

When N_d/N_p < 1, the system has more volumetric constraint equations than equilibrium equations, thus the system is over-constrained. In this case, if the u-P formulation I is used, the system equations will be very ill-conditioned so that it is hard to keep accuracy of solution and may cause divergence. If the u-P formulation II is used, the system equation will be singular because [K_PP] = [0] in this formulation so that the system is not solvable. Therefore, over-constrained models should be avoided as described in the Element Reference.

Volumetric constraint is incorporated into the final equations as extra conditions. A check is made at the element level for elements with internal nodes for pressure degrees of freedom and at degrees of freedom (HDSP) at global level for SOLID285 to see if the constraint equations are satisfied. The number of elements in which constraint equations have not been satisfied is reported for current-technology elements if the check is done at element level.

For u-P formulation I, the volumetric constraint is met if:

(3–115)

and for u-P formulation II, the volumetric constraint is met if:

(3–116)

and for u-P formulation III, the volumetric constraint is met if:

(3–117)

where:

When element formulations (or simply formulations) are mentioned in this reference, they refer to traditional forward solving, where input consists of an undeformed geometry/configuration, and the solution is the deformed geometry/configuration and the stresses and strains on it.

The term inverse formulation refers to formulations used for inverse solving, where input consists of a deformed geometry/configuration under loading, and the solution is the undeformed geometry/configuration and the stresses and strains on the deformed input geometry/configuration. Such formulations are therefore called inverse element formulations (or simply inverse formulations).

In this section, the term forward formulation appears (instead of simply formulation) to clearly differentiate between forward and inverse formulations.

For more information, see Nonlinear Static Analysis with Inverse Solving in the Structural Analysis Guide and Inverse-Solving Analysis of a Rotor Fan Blade with Disk in the Technology Showcase: Example Problems.

 

3.6.8.1. Displacement-Based Inverse Formulation

Displacement-based inverse formulation starts with the principal of virtual work (Equation 3–89 and Equation 3–90) used in forward formulation.

The diffentiation of :

(3–118)

is similar to Equation 3–96 but without the third term, as volume is the volume of the input configuration and is fixed during solution, as are the virtual strains (Govindjeet and Mihalic[[433]] and [[434]], and Lu and Li[[435]]).

The stress differentiation is the same as Equation 3–96, but the displacements for calculating and are the inverse displacements (rather than as in forward formulation). Here, and are coordinates in the input configuration and the current configuration, respectively.

The final displacement-based formulation is the same as Equation 3–99, but with:

Therefore, the stiffness matrices are unsymmetric in general, requiring the unsymmetric linear solver.