FreeFEM Documentation on GitHub

stars - forks

The System of elasticity

Elasticity

Solid objects deform under the action of applied forces:

a point in the solid, originally at \((x,y,z)\) will come to \((X,Y,Z)\) after some time; the vector \(\mathbf{u}=(u_1,u_2,u_3) = (X-x, Y-y, Z-z)\) is called the displacement. When the displacement is small and the solid is elastic, Hooke’s law gives a relationship between the stress tensor \(\sigma(u)=(\sigma_{ij}(u) )\) and the strain tensor \(\epsilon(u)=\epsilon_{ij}(u)\)

\[\sigma_{ij}(u) = \lambda \delta_{ij} \nabla.\mathbf{u}+ 2\mu\epsilon_{ij}(u),\]

where the Kronecker symbol \(\delta_{ij} = 1\) if \(i=j\), \(0\) otherwise, with

\[\epsilon_{ij}(u) = {1\over 2}({\partial u_i\over\partial x_j} + {\partial u_j\over\partial x_i} ),\]

and where \(\lambda, \mu\) are two constants that describe the mechanical properties of the solid, and are themselves related to the better known constants \(E\), Young’s modulus, and \(\nu\), Poisson’s ratio:

\[\mu = {E\over 2( 1+\nu)}, \quad \lambda = {E\nu\over (1+\nu)(1-2\nu)}.\]

Lamé’s system

Let us consider a beam with axis \(Oz\) and with perpendicular section \(\Omega\). The components along \(x\) and \(y\) of the strain \({\bf u}(x)\) in a section \(\Omega\) subject to forces \({\bf f}\) perpendicular to the axis are governed by:

\[-\mu \Delta {\bf u} - (\mu+\lambda) \nabla (\nabla .{\bf u})={\bf f}~~\hbox{in}~~\Omega,\]

where \(\lambda\) ,\(\mu\) are the Lamé coefficients introduced above.

Remark, we do not use this equation because the associated variational form does not give the right boundary condition, we simply use:

\[- div( \sigma ) = \mathbf{f} \quad \mbox{in}~~\Omega\]

where the corresponding variational form is:

\[\int_{\Omega} \sigma(u) : \epsilon(\mathbf{v})\;dx - \int_{\Omega} \mathbf{v} f \;dx =0;\]

where \(:\) denotes the tensor scalar product, i.e. \(a: b = \sum_{i,j} a_{ij}b_{ij}\).

So the variational form can be written as :

\[\int_{\Omega} \lambda \nabla.u \nabla.v + 2 \mu \epsilon(\mathbf{u}):\epsilon(\mathbf{v}) \; dx - \int_{\Omega} \mathbf{v} f \;dx =0;\]

Tip

Consider an elastic plate with the undeformed rectangle shape \([0,20]\times [-1,1]\).

The body force is the gravity force \(\mathbf{f}\) and the boundary force \(\mathbf{g}\) is zero on lower, upper and right sides. The left vertical side of the beam is fixed. The boundary conditions are:

\[\begin{split}\begin{array}{rcll} \sigma . {\bf n} &= \mathbf{g} &= 0 & \hbox{ on }\Gamma_1, \Gamma_4, \Gamma_3, \\ {\bf u} &= \mathbf{0} && \hbox{ on }\Gamma_2 \end{array}\end{split}\]

Here \({\bf u}=(u,v)\) has two components.

The above two equations are strongly coupled by their mixed derivatives, and thus any iterative solution on each of the components is risky. One should rather use FreeFEM’s system approach and write:

 1 // Parameters
 2 real E = 21e5;
 3 real nu = 0.28;
 4 
 5 real f = -1;
 6 
 7 // Mesh
 8 mesh Th = square(10, 10, [20*x,2*y-1]);
 9 
10 // Fespace
11 fespace Vh(Th, P2);
12 Vh u, v;
13 Vh uu, vv;
14 
15 // Macro
16 real sqrt2=sqrt(2.);
17 macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] //
18 // The sqrt2 is because we want: epsilon(u1,u2)'* epsilon(v1,v2) = epsilon(u): epsilon(v)
19 macro div(u,v) ( dx(u)+dy(v) ) //
20 
21 // Problem
22 real mu= E/(2*(1+nu));
23 real lambda = E*nu/((1+nu)*(1-2*nu));
24 
25 solve lame([u, v], [uu, vv])
26    = int2d(Th)(
27         lambda * div(u, v) * div(uu, vv)
28       + 2.*mu * ( epsilon(u,v)' * epsilon(uu,vv) )
29    )
30    - int2d(Th)(
31         f*vv
32    )
33    + on(4, u=0, v=0)
34    ;
35 
36 // Plot
37 real coef=100;
38 plot([u, v], wait=1, ps="lamevect.eps", coef=coef);
39 
40 // Move mesh
41 mesh th1 = movemesh(Th, [x+u*coef, y+v*coef]);
42 plot(th1,wait=1,ps="lamedeform.eps");
43 
44 // Output
45 real dxmin = u[].min;
46 real dymin = v[].min;
47 
48 cout << " - dep. max x = "<< dxmin << " y=" << dymin << endl;
49 cout << "   dep. (20, 0) = " << u(20, 0) << " " << v(20, 0) << endl;

The output is:

1 -- square mesh : nb vertices  =121 ,  nb triangles = 200 ,  nb boundary edges 40
2 -- Solve :           min -0.00174137  max 0.00174105
3          min -0.0263154  max 1.47016e-29
4 - dep.  max   x = -0.00174137 y=-0.0263154
5    dep.  (20,0)  = -1.8096e-07 -0.0263154
6 times: compile 0.010219s, execution 1.5827s

Solution of Lamé’s equations for elasticity for a 2D beam deflected by its own weight and clamped by its left vertical side is shown Fig. 19 and Fig. 20. Result are shown with a amplification factor equal to 100. The size of the arrow is automatically bound, but the color gives the real length.

LameVector

Fig. 19 Vector

LameDeformation

Fig. 20 Deformation

Elasticity

Table of content