Nonlinear finite elements/Updated Lagrangian formulation

Updated Lagrangian formulationEdit

We will now derive the finite element equations for the updated Lagrangian formulation for three-dimensional problems in solid mechanics.

Updated Lagrangian governing equationsEdit

The updated Lagrangian equations are

 

The boundary conditions are

 

The initial conditions are

 

Dependent variablesEdit

Note that the dependent variables in this case are all expressed as functions of the Lagrangian (material) coordinates, i.e.,

  1. Velocity ( )
  2. Cauchy stress ( )
  3. Rate of deformation ( )

Also note that instead of   and  , we may use   and   in our calculations and then transform to  .

Weak form of the balance of linear momentumEdit

Another name for the weak form is the Principle of Virtual Power (or Principle of Virtual Work if rates are not involved).

The strong form consists of the balance of linear momentum and the traction boundary conditions. That is,

 

where   is the acceleration.

The weak form of this equation can be written as

 

where   is a vector valued weighting function (also called a test function) that is zero at the points on the boundary where the essential boundary conditions on   are applied. Therefore we can think of   as a variation of  .

Derivation of the weak formEdit

We proceed in the usual manner to derive the weak form. We multiply the differential equation with the weighting function and integrate over the volume to get

 

Using the identity

 

and noting that   (conservation of angular momentum), we have

 

From the divergence theorem we have

 

where   is the outward unit normal to the boundary. So,

 

Now, on the boundary

 

On the part of the boundary where no tractions are applied, we have   and hence  . Therefore we can write

 

where   is the part of the boundary where the tractions are applied.

Plugging this expression into the weighted average form of the balance of linear momentum we get

 

Rearranging terms we get the required weak form

 

Physical interpretation of termsEdit

To provide a physical interpretation of the four terms in the weak form, we write the weighting function as a variation of the velocity, i.e.,  . Then we have the familiar principle of virtual power:

 

The first term on the left can then be identified as the virtual kinetic power, the second term on the left as the virtual internal power, and the sum of the two terms on the right as the virtual external power. In other words, we have

 

and

 

We can see why the second terms is identified as the virtual internal power as follows. Note that

 

where   is the velocity gradient,   is symmetric rate of deformation tensor, and   is the skew symmetric spin tensor. Therefore, since the product of a symmetric and a skew symmetric tensor is zero, we have

 

Hence,

 

This has the same form as the expression for internal power that we discussed earlier.

Finite element discretizationEdit

Let us go back to the original weak form that we derived, i.e.,

 

Trial functionEdit

To start the discretization process we have to first choose an approximate solution   from the space   of admissible finite element trial functions. These functions must satisfy at least   continuity and the essential boundary condition   on  .

More formally, we choose   where

 

If we discretize the body into a number of elements with   nodes, then we can write the approximate solution as

 

The quantities   are the shape functions which provide the required level of continuity and are chosen so that  .

Test functionEdit

We also have to choose appropriate weighting functions ( ) which are also known as test functions.

In Galerkin finite elements we choose these functions from the same space as the trial functions with the additional restriction that these functions go to zero at points on the boundary where essential boundary conditions are applied.

Formally, choose   where

 

For finite element analysis, these weighting functions have the same form as the trial function, i.e.,

 

MotionEdit

Instead of a trial function for the velocity, we can start with one for the motion and then take time derivatives. In this case, we choose the motion to be approximated by

 

where   are the nodal shape functions (note that these are functions of  ) and   are the positions of nodes   in the current configuration. Therefore, nodes remain coincident with material point labels at all times.

Therefore, in the reference configuration, we have

 

The shape functions are also chosen such that

 

Then the displacement   can be approximated by

 

The velocity   is given by the time derivative of   keeping   fixed. So we have

 

Note that this approximation has the same form as the trial function for the velocity that we started with.

The acceleration   is given by the time derivative of  . Hence

 

Velocity gradient and rate of deformationEdit

The spatial velocity gradient   is given by

 

Note that the derivatives are with respect to   and not  .

Therefore, the approximation that we seek has the form

 

Using the identity

 

we get

 

In index notation,

 

The rate of deformation is then given by

 

In index notation,

 

Approximate weak formEdit

We can now substitute the trial and test functions into the weak form of the momentum equation to get (using the same procedure as in one-dimension)

 

We can write this approximate weak form as

 

where the first term represents the inertial force, the second term the internal force, and the third term the external force. Thus,

 

The mass matrixEdit

The term

 

is the inertial force term.

If we substitute

 

into the inertial force terms, we get

 

We define the mass matrix as the quantity

 

Therefore, the inertial force can be expressed, in analogy with Newton's second law, as

 

Semi-discrete finite element equationsEdit

We can then write the semi-discrete version of the weak form in the shape of a matrix equation as

 

These equations are ordinary differential equations because we still have time derivatives on the left hand side - hence the system of equations is semi-discrete.

There are a number of ways in which the discretization in time can be effected. Some common approaches are the generalized midpoint rule, the Newmark   method, or Runge-Kutta explicit integration schemes. There is a large literature on the best way of doing the discretization in time with the goal of conserving both momentum and energy as accurately as possible.

Taking derivatives and calculating integralsEdit

At this stage we are faced with the problem of calculating derivatives and integrals with respect to   of quantities that depend on  . How can we proceed?

The standard way of resolving this issue is to do all our calculations with respect to a parent element and then map the results to the reference or current configurations as necessary.

To make things more concrete let us consider a single element and label the three domains as

  1. Parent element ( ).
  2. Reference element ( ).
  3. Current element ( ).

Let the maps that we need be

  1. Parent   Reference :  .
  2. Parent   Current :  .
  3. Reference   Current :  .

This situation is illustrated (for a two dimensional situation) in the following figure.

 
Maps from parent element to reference and current configurations

We define the map between the parent element and the element in the current configuration as

 

Then,

 

Let us now try to compute the spatial velocity gradient. We have

 

The quantity

 

is called the Jacobian of the motion and is a function of time.

Inverting, we get

 

But the spatial velocity gradient is

 

Therefore,

 

If we substitute the approximation

 

we get a formula for the velocity gradient

 

A similar procedure can be followed when other spatial gradients need to be calculated. The most frequently encountered situation is the calculation of the gradient of the shape functions. In that case, we have

 


For computing integrals, we observe that for a function   on   we can write

 

where

 

Similarly for a function   over  , we have

 

where

 

Recall that the internal force terms have the form

 

Expressed in terms of an integral over the parent element, we then have

 

where

 

Finite element implementationEdit

Recall that the finite element system of equations can be written as

 

where

 

To get a feel for the structure of the equation, let us consider a four-noded plane element. If the components of the velocity in the two coordinate directions are  , the acceleration vector has the form

 

The mass matrix has the form

 

In expanded form

 

Next, we consider the internal force term. Let the components of the internal force in the two coordinate directions be ( ). Then,

 

At this stage, recall that the velocity gradient is given by

 

In index notation,

 

We define

 

Then,

 

In matrix form

 

Expanded out

 

In more compact form, we can just write

 

Similarly, for the internal force term,

 

we have

 

In matrix form,

 

Expanded out,

 

In compact form, we have

 

We can express the external force term in a similar manner. Notice that we are not taking adavantage of the symmetry of the stress tensor in the above expressions.

A widely used alternative way of expressing the finite element system of equations is the Voigt notation. This notation can be found in discussions of introductory finite element analysis.

Algorithm for computing internal forcesEdit

The major complications and computational effort in the finite element implementation are encountered during the computation of internal forces.

The basic procedure for computing the internal force at the nodes of an element is given below.


  1. Set  .
  2. For all quadrature points in the parent element ( )
    1. Compute the gradient of the shape functions ( ) for all nodes  .
    2. Compute the velocity gradient  .
    3. Compute the rate of deformation tensor.
    4. Compute the deformation gradient/ the Lagrangian Green tensor.
    5. Compute the Cauchy stress or the 2nd P-K stress using the constitutive equation.
    6. Update the nodal internal force ( ) using   where   are the weights for Gaussian integration.


The mass matrix and integration in time will be discussed later.