1   The equations of linear elasticity (100 pts)

Analysis of structures is one of the major activities of modern engineering, which likely makes the PDE modeling the deformation of elastic bodies the most popular PDE in the world. In this assignment, we will develop code to solve the equations of 2D or 3D elasticity in MATLAB. This will also be your first exposure to vector PDEs, i.e., the unknown u will not be a scalar at each node, but instead will be a d -dimensional vector.

1.1   PDE problem

The equations governing small elastic deformations of a body Ω can be written as,

σ   =   λ tr(e)I + 2µe                                                         (2)

e   =    ╱Vu + (Vu)T 、                                                     (3)

where σ is the stress tensor, f is the body force per unit volume, λ and µ are Lamé’s elasticity parameters for the material in Ω, I is the identity tensor, tr() is the trace operator on a tensor, e is the symmetric strain-rate tensor (symmetric gradient), and u is the displacement vector field. We have here assumed isotropicelastic conditions.

We combine (2) and (3) to obtain

σ = λ(V．u)I + µ(Vu + (Vu)T ).                                                   (4)

Note that (1)– (3) can easily be transformed to a single vector PDE for u, which is the governing PDE for the unknown u (Navier’s equation).  In the derivation of the variational formulation, however, it is convenient to keep the equations split as above.

1.2   Variational Formulation

The variational form of (1)– (3) consists of forming the inner product of (1) and a vector test func- tion v ← Vˆ , where Vˆ is a vector-valued test function space, and integrating over the domain Ω :

− lΩ (V．σ)．v dx = lΩ fv dx.

Since V．σ contains second-order derivatives of the primary unknown u, we integrate this term by parts:

− lΩ (V．σ)．v dx = lΩ σ : Vv dx − l合Ω (σ．n)．v ds.

where the colon operator is the inner product between tensors (summed pairwise product of all elements), and n is the outward unit normal at the boundary. The quantity σ．n is known as the traction or stress vector at the boundary, and is often prescribed as a boundary condition. We here assume that it is prescribed on a part ∂ΩT of the boundary as σ．n = T. On the remaining part of the boundary, we assume that the value of the displacement is given as a Dirichlet condition. We

thus obtain

lΩ σ : Vv dx = lΩ fv dx + l合ΩT (σ．n)．v ds.

Inserting the expression (4) for σ gives the variational form with u as unknown. Note that the boundary integral on the remaining part ∂Ω \ ∂ΩT vanishes due to the Dirichlet condition.

We can now summarize the variational formulation as: find u ← V such that

a(u, v) = L(v)       Av ← V,                                                       (5)

where

a(u, v)   =   lΩ σ(u) : Vv dx

σ(u)   =   λ(V．µ)I + µ(Vu + (Vu)T )

L(v)   =   lΩ fv dx + l合ΩT (σ．n)．v ds.

(6)

(7)

One can show that the inner product of a symmetric tensor A and an anti-symmetric tensor B vanishes. If we express Vv as a sum of its symmetric and anti-symmetric parts, only the symmetric part will survive in the product σ : Vv since σ is a symmetric tensor. Thus replacing Vu by the symmetric gradient e(u) gives rise to the slightly different variational form

a(u, v) = lΩ σ(u) : e(v) dx,

where e(v) is the symmetric part of Vv:

(8)

e(v) =  ╱Vv + (Vv)T、

The formulation (8) is what naturally arises from minimization of elastic potential energy and is a more popular formulation than (6).

1.3   Test Problem

As a test problem, you will model a clamped beam deformed under its own weight in 3D. This can be modeled by setting the right-hand side body force per unit volume to f = (0, 0, âL´ Šρg) with ρ the density of the beam and g the acceleration of gravity. The beam is box-shaped with length L and has a square cross section of width W. We set u = uD  = (0, 0, 0) at the clamped end, x = 0. The rest of the boundary is traction free; that is, we set T = 0.

Geometry   You need to create a simple geometry of the beam.  Start with L = 1, W = 0.2 and change these later.  It would be easiest to consider a mesh made of hexahedral elements, i.e., boxes. You should specify the number of elements in x, N! and the number of elements in y/z, N↓# .  Use the sample 1D code we did in class as a template for designing the appropriate data structures. You should use linear elements, so each element will have 8 nodes. You can start with N! = 20, N↓# = 4 and adjust accordingly based on the capabilities of your machine. Note that you should make the matrices sparse as it will be very expensive with full matrices.

Here are simple parameters to start with. Experiment with these.

Figure 1: The geometry for our test problem.

● L = 1, W = 0.2

●  µ = 1, ρ = 1, λ = 1.25

●  g = 10

Stress Computation   As soon as the displacement u is computed, we can compute various stress

measures. We will compute the von Mises stressdefined as σM = ′ s : s, where s is the deviatoric

stress tensor

s = σ —  tr(σ)I.

You should visualize the von Mises stress and the displacements. You can plot the slice through the beam along with the displacements using the matlab command quiver. An example, rendered using paraview in 3D is shown below.

● Implement the FEM to solve the linear elasticity problem for the clamped beam in Matlab.

● Compute the steady state displacements u for L = 1, W = 0.2, µ = 1, ρ = 1, λ = 1.25. Plot the magnitude of u along with the displacement vectors.

● Compute and plot the von Mises stress.

● what is the relationship between the displacement in the bar and changes in µ and λ . Change these, observe and report.