Tutorial 2: a toy problem
We consider a regular domain Ω ⊂ Rd (d = 2,3), filled with a linear isotropic elastic material. This body is clamped on a part ΓD ⊂ ∂Ω of its boundary, and is submitted to body forces f, as well as to traction loads g, applied on a part ΓN ⊂ ∂Ω disjoint from ΓD. The remaining, traction-free region Γ:=∂Ω\(ΓD ∪ ΓN ) is called the free boundary of Ω.
The deformation of the body Ω is the unique solution u to the linear elasticity system:
−div(Ae(u)) = f on Ω, f are body forces applied on the shape
u = 0 on ΓD, ΓD⊂∂Ω is a clamping region
Ae(u)n = g on ΓN, (1) where g are surface loads applied on a subset ΓN⊂∂Ω
Ae(u)n = 0 on Γ, Γ ⊂ ∂Ω is traction-free
where e(u) := (∇u^t+∇u)/2 is the linearized strain tensor, and n is the unit normal vector to ∂Ω, pointing outward Ω. The mechanical model uses also the stress tensor A which is related to e(u) by Hooke’s law
Ae(u) = 2μe(u) + λ tr(e(u)) I
where λ and μ are the Lamé coefficients of the material.
Equation (1) can be also written in terms of the modulus of elasticity (Young modulus) E and the Poisson ratio (contraction ratio) ν. These constants are related to λ and μ via
E = μ (3λ + 2μ) / (λ+μ) ν = λ / (2(λ+μ))
λ = Eν / ((1+ν)(1-2ν)) μ = E / (2(1+ν))
If ΓD = ∅ (resp. ΓN = ∅), the boundary value problem (2) is called a pure traction (resp. pure displacement) problem.
Suppose we like to solve problem (1) in a 2d domain discretized by a mesh Th shown hereafter.
A unit square plate [0,1]x[0,1] fixed at its bottom side labelled '1' and a constant traction load g is applied perpendicularly (along th normal direction) on its top side labelled '2'. The body is subjected to gravity force. Lamé coefficients λ and μ are respectively set to 18600 and 3400.
(Note: lateral sides of the square are not subjected to any surface force (g=0 in the equation) and are free to move. Such natural Neumann condition will be automatically satisfied by the program and do not need to be specified).
The parameter file carre.elas describing this test conditions must then contain the following lines:
Dirichlet
1 # one boundary condition specified
1 vertex v 0 0 # vertices labelled '1' set to value (0,0)
Load
1 # one traction force applied
2 edges n 1000 # edges labelled '2' enjoy a force of modulus 1000 perpendicularly
gravity
0 -9.81e2
Lame
1
0 186000. 3400. # material properties assigned to triangles labelled '0'
Running the program elastic with the following arguments
elastic carre
or equivalently
elastic carre.mesh -p carre.elas -o carre.sol
produces the following displacement field:
showing the isovalues of the modulus of the displacement (left) and the deformed body (right).
Similar results can be obtained on a 3d body. Go to the directory demos 3d:
cd demos/3d
and run elastic with the billed example
elastic biel3d.mesh
This will produces the output file shown below:
