Finite Element Formulation

Finite Element Approximation

In the finite element method, the motion

x (X, t)

is approximated by:(1)

x_{i} (X, t) = Φ_{I} (X) x_{i I} (t)

Where,

$Φ_{I} (X)$: Interpolating shape functions
$x_{i I}$: Position vector of node $I$

Summation over repeated indices is implied. In the case of lower indices, summation is over the number of space dimensions. For upper case indices, summation is over the number of nodes. The nodes in the sum depend on the type of entity considered. When the volume is considered, the summation is over all the nodes in the domain. When an element is considered, the sum is over the nodes of the element.

Similarly, nodal displacements are defined using Material and Spatial Coordinates, Equation 2 at nodes:(2)

u_{i I} (X, t) = x_{i I} (t) - X_{i I}

The displacement field is:(3)

u_{i} (X, t) = Φ_{I} (X) u_{i I} (t)

The velocities are obtained by taking the material time derivative of the displacement giving:(4)

v_{i} (X, t) = \frac{\partial u_{i} (X, t)}{\partial t} = Φ_{I} (X) v_{i I} (t)

It is worth pointing out that the velocity is a material time derivative of displacements, that is, the partial derivative with respect to time with the material coordinate fixed.

Finally, accelerations are similarly given by the material time derivative of velocities:(5)

{\dot{v}}_{i} (X, t) = Φ_{I} (X) {\dot{v}}_{i I} (t)

Emphasis is placed on the fact that shapes functions are functions of the material coordinates whatever the updated or the total Lagrangian formulation is used. All the dependency in the finite element approximation of the motion is taken into account in the values of the nodal variables.

From Kinematic Description, Equation 8, the velocity gradient is given by:(6)

L_{i j} = \frac{\partial v_{i}}{\partial x_{j}} = v_{i I} \frac{\partial Φ_{I}}{\partial x_{j}} = v_{i I} Φ_{I, j}

and the rate of deformation (Kinematic Description, Equation 14) by:(7)

D_{i j} = \frac{1}{2} (L_{i j} + L_{j i}) = \frac{1}{2} (v_{i I} Φ_{I, j} + v_{j I} Φ_{I, i})

Similarly, the test functions are approximated as:(8)

δ v_{i} (X) = Φ_{I} (X) δ v_{i I}

Where, $δ v_{i I}$ are the virtual nodal velocities.

The test functions are next substituted into the principle of virtual power (Virtual Power Principle, Equation 5) giving:(9)

δ v_{i I} \int_{Ω} \frac{\partial Φ_{I}}{\partial x_{j}} σ_{j i} d Ω - δ v_{i I} \int_{Ω} Φ_{I} ρ b_{i} d Ω - δ v_{i I} \int_{Γ_{σ}} Φ_{I} τ_{i} d Γ + δ v_{i I} \int_{Ω} Φ_{I} ρ {\dot{v}}_{i} d Ω = 0

The virtual velocities must be kinematically admissible, that is, satisfy boundary conditions on

Γ_{u}

, the part of the boundary where kinematical conditions are specified. Using the arbitrariness of the virtual nodal velocities everywhere except on

Γ_{u}

, the weak form of the momentum equation is:(10)

\int_{Ω} \frac{\partial Φ_{I}}{\partial x_{j}} σ_{j i} d Ω - \int_{Ω} Φ_{I} ρ b_{i} d Ω - \int_{Γ_{σ}} Φ_{I} τ_{i} d Γ + \int_{Ω} Φ_{I} ρ {\dot{v}}_{i} d Ω = 0

with $Γ_{σ}$ the part of the boundary where traction loads are imposed.

Internal and External Nodal Forces

As in Virtual Power Term Names, you define the nodal forces corresponding to each term in the virtual power equation.

The internal nodal forces are defined by:(11)

δ P^{int} = δ v_{i I} f_{i I}^{int} = \int_{Ω} \frac{\partial δ v_{i}}{\partial x_{j}} σ_{j i} d Ω = δ v_{i I} \int_{Ω} \frac{\partial Φ_{I}}{\partial x_{j}} σ_{j i} d Ω

The stress is the true (Cauchy) stress.(12)

f_{i I}^{int} = \int_{Ω} σ_{j i} (\frac{\partial Φ_{I}}{\partial x_{j}}) d Ω

These nodal forces are called internal because they represent the stresses in the body. The expression applies to both the complete mesh or to any element. It is pointed out that derivatives are taken with respect to spatial coordinates and that integration is taken over the current deformed configuration.

The external forces are similarly defined in terms of the virtual external power:(13)

δ P^{e x t} = δ v_{i I} f_{i I}^{e x t} = δ v_{i l} \int_{Ω} Φ_{I} ρ b_{i} d Ω + δ v_{i l} \int_{Γ_{σ}} Φ_{I} τ_{i} d Γ

so that external forces are given by:(14)

f_{i I}^{e x t} = \int_{Ω} Φ_{I} ρ b_{i} d Ω + \int_{Γ_{σ}} Φ_{I} τ_{i} d Γ

Mass Matrix and Inertial Forces

The inertial body forces are defined by:(15)

δ P^{i n e r t} = δ v_{i I} f_{i I}^{i n e r t} = δ v_{i I} \int_{Ω} Φ_{I} ρ {\dot{v}}_{i} d Ω = δ v_{i I} \int_{Ω} Φ_{I} ρ {\dot{v}}_{i} d Γ

so that the inertia forces are given by:(16)

f_{i I}^{i n e r t} = \int_{Ω} Φ_{I} ρ {\dot{v}}_{i} d Ω

or using the Finite Element Approximation, Equation 5 for the accelerations:(17)

f_{i I}^{i n e r t} = \int_{Ω} ρ Φ_{I} Φ_{J} d Ω {\dot{v}}_{i J}

It is usual to define the inertial nodal forces as the product of a mass matrix and the nodal accelerations. Defining the mass matrix as:(18)

M_{i j I J} = δ_{i j} \int_{Ω} ρ Φ_{I} Φ_{J} d Ω

the inertial forces are given by:(19)

f_{i I}^{i n e r t} = M_{i j I J} {\dot{v}}_{j J}

Discrete Equations

Using the definitions of the internal and external forces, as well as the definition of the inertial forces, it is possible to write the weak form of the virtual power principle as:(20)

δ v_{i I} (M_{i j I J} {\dot{v}}_{j J} + f_{i I}^{int} - f_{i I}^{e x t}) = 0

or taking into account the arbitrariness of the virtual velocities:(21)

M_{i j I J} {\dot{v}}_{j J} + f_{i I}^{int} = f_{i I}^{e x t}

Equation of Motion for Translational Velocities

Discrete Equations, Equation 21 is written in matrix notation as:(22)

M \frac{d v}{d t} = f^{e x t} - f^{int}

This is Newton's equation, where the mass matrix is:(23)

M = \int_{Ω} ρ Φ^{T} Φ d Ω

Radioss uses a lumped mass approach, that is, each node represents a discrete mass of zero size. This creates a diagonal mass matrix

M

, eliminating, as you will see in Dynamic Analysis, the need to solve simultaneous equations for the solution of nodal accelerations.(24)

f^{e x t} = \int_{Γ} Φ^{T} τ d Γ + \int_{Ω} ρ Φ^{T} b d Ω

is the externally applied load vector, and the internal force vector is:(25)

f_{i I}^{int} = \int_{Ω} σ_{i j} \frac{\partial Φ_{I}}{\partial x_{j}} d Ω

Adding to the internal and external forces the anti-hourglass force vector and the contact force vector which will be described in the following sections, you obtain the overall equation of motion:(26)

M \frac{d v}{d t} = f^{e x t} - f^{int} + f^{h g r} + f^{c o n t}

Equation of Motion for Angular Velocities

Shell, beam and rigid body theory introduces nodal rotational degrees of freedom. The equations of motion for rotational degrees of freedom are complicated if written in the global reference frame. They are much simpler if written for each node in the principal reference frame attached to the node. The resulting equations are the standard Euler equations. They are completely analogous to Newton's law governing translational degrees of freedom and are stated as:(27)

l_{1} α_{1} + (I_{3} - I_{2}) ω_{2} ω_{3} = m_{1}^{e x t} - m_{1}^{int}

(28)

l_{2} α_{2} + (I_{1} - I_{3}) ω_{1} ω_{3} = m_{2}^{e x t} - m_{2}^{int}

(29)

l_{3} α_{3} + (I_{2} - I_{1}) ω_{1} ω_{2} = m_{3}^{e x t} - m_{3}^{int}

Where,

$I_{1}, I_{2}, I_{3}$: Principal moments of inertia about the x, y and z axes respectively
$α 1, α 2, α 3$: Angular accelerations expressed in the principal reference frame
$ω 1, ω 2, ω 3$: Angular velocities
$m_{1}^{e x t}, m_{2}^{e x t}, m_{3}^{e x t}$: Principal externally applied moments
$m_{1}^{int}, m_{2}^{int}, m_{3}^{int}$: Principal internal moments

The equation of motion for rotational degrees of freedom is thus very similar to that for translational degrees of freedom. In matrix notation and in the nodal principal reference frame:(30)

I \frac{d ω}{d t} = M^{e x t} - M^{int} + F (ω)

The vector function $F (ω)$ is computed for a value of $ω$ at $t - δ t / 2$ . Equation 30 is used for rigid body motion.

For shell, beam and spring using a spherical inertia, the equation of motion becomes:(31)

I \frac{d ω}{d t} = M^{e x t} - M^{int} + M^{h g r}

Where,

$I = \sum_{e l e m e n t s} I_{e}$: Diagonal inertia matrix
$M^{e x t} = \sum_{e l e m e n t s} m^{e x t}$: Externally applied moment vector
$M^{int} = \sum_{e l e m e n t s} m^{int}$: Internal moment vector
$M^{h g r} = \sum_{s h e l l s} m^{h g r}$: Anti-hourglass shell moment vector

Element Coordinates

Finite elements are usually developed with shape functions expressed in terms of an intrinsic coordinates system $ξ, η, ζ$ . It is shown hereafter that expressing the shape functions in terms of intrinsic coordinates is equivalent to using material coordinates.

When an element is treated in terms of intrinsic coordinates, we are concerned with three domains that correspond to this element:

The domain in the intrinsic coordinates system
The current element domain
The initial reference element domain

$ξ$ is associated with the direction $u$

$η$ is associated with the direction $v$

$ς$ is associated with the direction $w$

The motion in each element can thus be described by the composition of three maps (the reasoning is described only for the direction

u

):

The map from the intrinsic coordinates system to the initial configuration:
(32) $X (ξ)$
The map from the intrinsic coordinates system to the current configuration:
(33) $x (ξ, t)$
The map from the initial to the current configuration:
(34) $x = ϕ (X, t)$

So, it is possible to approximate the motion in an element by:(35)

x_{i} (ξ, t) = Φ_{I} (ξ) x_{i I} (t)

Shape functions

Φ_{I} (ξ)

have no dimensions. They simply relate coordinates in the physical world to the intrinsic coordinates system. Writing Equation 35 at

t

=0, you obtain:(36)

x_{i} (ξ) = x_{i} (ξ, 0) = Φ_{I} (ξ) x_{i I} (0) = Φ_{I} (ξ) x_{i I}

So, it can be seen from the last equation that the material coordinates system and the intrinsic coordinates system are invariant in a Lagrangian element. As a result, as intrinsic coordinates are time invariant and it is possible to write displacements, velocities and accelerations in terms of intrinsic coordinates (one coordinate system, the two other coordinates have similar shape functions):(37)

u_{i} (ζ, t) = Φ_{I} (ξ) u_{i I} (t)

(38)

{\dot{u}}_{i} (ζ, t) = Φ_{I} (ξ) {\dot{u}}_{i I} (t)

(39)

{\dot{v}}_{i} (ζ, t) = Φ_{I} (ξ) {\dot{v}}_{i I} (t)

Isoparametric elements use the same shape functions for the interpolation of $x, u, \dot{u}$ and $\dot{v}$ .

Integration and Nodal Forces

In practice, integrals over the current domain in the definition of the internal nodal forces (Equation of Motion for Translational Velocities, Equation 25), of the external nodal forces (Equation of Motion for Translational Velocities, Equation 24) and of the mass matrix have to be transformed into integrals over the domain in the intrinsic coordinate system $Δ$ .

Using Vicinity Transformation, Equation 5, integrals on the current configuration are related to integrals over the reference configuration and over the domain in the intrinsic coordinate system by:(40)

\int_{Ω} g (x) d Ω = \int_{Ω^{0}} g (x) | F | d Ω_{0} = \int_{Δ} g (ξ) | F_{ξ} | d Δ

and(41)

\int_{Ω_{0}} g (X) d Ω_{0} = \int_{Δ} g (ξ) | F_{ξ}^{0} | d Δ

Where,

$| F |$: The Jacobian determinant of the transformation between the current and the initial configuration
$| F_{ξ} |$: The Jacobian determinant of the transformation between the current configuration and the domain in the intrinsic coordinate system
$| F_{ξ}^{0} |$: The Jacobian determinant of the transformation between the reference configuration and the intrinsic coordinate system

On the other hand, it comes from Vicinity Transformation, Equation 2 and Element Coordinates, Equation 35:(42)

F_{ξ k j} = \frac{\partial x_{k}}{\partial ξ_{j}} = \frac{\partial Φ_{I} (ξ)}{\partial ξ_{j}} x_{k l}

So, using Equation 40, internal forces computed by integration over the current domain will be obtained by the following quadrature:(43)

f_{i I}^{int} = \int_{Ω} σ_{i j} \frac{\partial Φ_{I}}{\partial x_{j}} d Ω = \int_{Δ} σ_{i j} \frac{\partial Φ_{I}}{\partial x_{j}} | F_{ξ} | d Δ

and $| F_{ξ} |$ obtained from Equation 42.

External forces and the mass matrix can similarly be integrated over the domain in the intrinsic coordinate system.

Function Derivatives

The definition of internal forces also shows that derivatives of the form:(44)

\frac{\partial}{\partial x_{j}}

need to be computed. These spatial derivatives are obtained by implicit differentiation. Considering the velocity gradient such as:(45)

L_{i j} = \frac{\partial v_{i}}{\partial x_{j}}

one has:(46)

L_{i j} = \frac{\partial v_{i}}{\partial ξ_{k}} \frac{\partial ξ_{k}}{\partial x_{j}} = \frac{\partial v_{i}}{\partial ξ_{k}} F_{ξ k j}^{- 1} = v_{i I} \frac{\partial Φ_{I}}{\partial ξ_{k}} F_{ξ k j}^{- 1}

Where, the Jacobian matrix of the map between the current coordinates and the intrinsic coordinates is:(47)

F_{ξ k j} = \frac{\partial x_{k}}{\partial ξ_{j}} = \frac{\partial Φ_{I} (ξ)}{\partial ξ_{j}} x_{k I}

Usually, it is not possible to have closed form expression of the Jacobian matrix. As a result the inversion will be performed numerically and numerical quadrature will be necessary for the evaluation of nodal forces.

Numerical Quadrature: Reduced Integration

All elements in Radioss are integrated numerically. Hence, the integrals for nodal forces are replaced by a summation:(48)

\int f (ξ) d ξ = \sum_{j = 1}^{n} w_{j} f (ξ_{j})

Where,

$n$: Number of integration points in the element
$w_{j}$: Eeight associated to the integration point $j$

Values of

w_{j}

and locations of

ξ_{j}

are given in tables according to the numerical quadrature approach. Radioss uses either full or reduced integration schemes.

For full integration, the number of integration points is sufficient for the exact integration of the virtual work expression. The full integration scheme is often used in programs for static or dynamic problems with implicit time integration. It presents no problem for stability, but sometimes involves locking and the computation is often expensive.

Reduced integration can also be used. In this case, the number of integration points is sufficient for the exact integration of the contributions of the strain field that are one order less than the order of the shape functions. The incomplete higher order contributions to the strain field present in these elements are not integrated.

The reduced integration scheme, especially with one-point quadrature is widely used in programs with explicit time integration to compute the force vectors. It drastically decreases the computation time, and is very competitive if the spurious singular modes (often called hourglass modes which result from the reduced integration scheme) are properly stabilized. In two dimensions, a one point integration scheme will be almost four times less expensive than a four point integration scheme. The savings are even greater in three dimensions. The use of one integration point is recommended to save CPU time, but also to avoid "locking" problems, for example, shear locking or volume locking.

Shear locking is related to bending behavior. In the stress analysis of relatively thin members subjected to bending, the strain variation through the thickness must be at least linear, so constant strain first order elements are not well suited to represent this variation, leading to shear locking. Fully integrated first-order isoparametric elements (tetrahedron) also suffer from shear locking in the geometries where they cannot provide the pure bending solution because they must shear at the numerical integration points to represent the bending kinematic behavior. This shearing then locks the element, that is, the response is far too stiff.

On the other hand, most fully integrated solid elements are unsuitable for the analysis of approximately incompressible material behavior (volume locking). The reason for this is that the material behavior forces the material to deform approximately without volume changes. Fully-integrated solid elements, and in particular low-order elements do not allow such deformations. This is another reason for using selectively reduced integration. Reduced integration is used for volume strain and full integration is used for the deviatoric strains.

However, as mentioned above, the disadvantage of reduced integration is that the element can admit deformation modes that are not causing stresses at the integration points. These zero-energy modes make the element rank-deficient which causes a phenomenon called hour-glassing; the zero-energy modes start propagating through the mesh, leading to inaccurate solutions. This problem is particularly severe in first-order quadrilaterals and hexahedra.

To prevent these excessive deformations, a small artificial stiffness or viscosity associated with the zero-energy deformation modes is added, leading in Equation of Motion for Translational Velocities, Equation 22 and Equation of Motion for Angular Velocities, Equation 30 to anti-hourglass force and moment vectors:(49)

M \frac{d v}{d t} = f^{e x t} - f^{int} + f^{h g r}

(50)

I \frac{d ω}{d t} = M^{e x t} - M^{int} + M^{h g r}

Zero-energy or hourglass modes are controlled using a perturbation stabilization as described by Flanagan-Belytschko ¹, or physical stabilization as described in ² (Element Library).

So, for isoparametric elements, reduced integration allows simple and cost effective computation of the volume integrals, in particular on vectorized supercomputers, and furnishes a simple way to cope with locking aspects, but at the cost of allowing hour-glassing.

Numerical Procedures

The Radioss numerical solver can be summarized by the flow chart in Figure 1. For each time step in a particular analysis, the algorithm used to compute results is:

For the displacement, velocity and acceleration at a particular time step, the external force vector is constructed and applied.
A loop over element is performed, in which the internal and hourglass forces are computed, along with the size of the next time step. The procedure for this loop is:
1. The Jacobian matrix is used to relate displacements in the intrinsic coordinates system to the physical space:
  (51) ${\frac{\partial Φ}{\partial x_{j}} |}_{t} = {F_{ξ}^{- 1} \frac{\partial Φ}{\partial ξ} |}_{t}$
2. The strain rate is calculated:
  (52) ${\dot{ε}}_{i j} = (\frac{\partial Φ_{I}}{\partial x_{j}}) \dot{x} = \frac{1}{2} (\frac{\partial v_{i}}{\partial x_{j}} + \frac{\partial v_{j}}{\partial x_{i}})$
3. The stress rate is calculated:
  (53) ${\dot{σ}}_{i j} = f (\dot{ε}, m a t e r i a l - l a w)$
4. Cauchy stresses are computed using explicit time integration: (54) $σ (t + Δ t) = σ (t) + \dot{σ} Δ t$
5. The internal and hourglass force vectors are computed.
6. The next time step size is computed, using element or nodal time step methods (Dynamic Analysis)
After the internal and hourglass forces are calculated for each element, the algorithm proceeds by computing the contact forces between any interfaces.
With all forces known, the new accelerations are calculated using the mass matrix and the external and internal force vectors: (55) ${\dot{v}}_{i} = M^{- 1} (f_{e x t_{i}} - f_{{int}_{i}})$
Finally, time integration of velocity and displacement is performed using the new value.

Figure 1. Numerical Procedure