1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
\subsection{The incompressible elasticity problem}
\label{sec-incompressible-elasticity}
\cindex{problem!elasticity!incompressible}%
\cindex{boundary condition!Dirichlet}%
\cindex{boundary condition!Neumann}%
\cindex{boundary condition!mixed}%
\cindex{stabilization}%
\subsubsection*{Formulation}
Let us go back to the linear elasticity problem.
When $\lambda$ becomes large, this problem is related
to the \emph{incompressible elasticity} and
cannot be solved as it was previously done.
To overcome this difficulty, the pressure is introduced~:
\[
p = - \lambda {\rm div}\, {\bf u}
\]
and the problem becomes:
\ \ \ $(E)$ {\it find ${\bf u}$ and $p$ defined in $\Omega$ such that:}
\[
\begin{array}{ccccl}
-\ {\bf div}(2D({\bf u})) &+& \bnabla p &=& {\bf f} \ {\rm in}\ \Omega, \\
-\ {\rm div}\,{\bf u} &-& \Frac{1}{\lambda} p &=& 0 \ {\rm in}\ \Omega, \\
+ B. C.
\end{array}
\]
\cindex{form!{$2D({\bf u}):D({\bf v})$}}%
\cindex{form!{${\rm div}({\bf u})\,q$}}%
The variational formulation of this problem expresses:
\ \ \ $(VFE)$ {\it find ${\bf u}\in V(1)$ and $p \in L^2(\Omega)$ such that:}
\[
\begin{array}{lcccl}
a({\bf u},{\bf v}) &+& b({\bf v}, p) &=& m({\bf f},{\bf v}), \ \forall {\bf v}\in V(0), \\
b({\bf u},q) &-& c(p,q)&=& 0, \ \forall q \in L^2_0(\Omega),
\end{array}
\]
where
\begin{eqnarray*}
m({\bf u},{\bf v}) &=& \int_\Omega {\bf u} . {\bf v} \, dx, \\
a({\bf u},{\bf v}) &=& \int_\Omega D({\bf u}) : D({\bf v}) \, dx, \\
b({\bf v}, q) &=& - \int_\Omega {\rm div}({\bf v}) \, q \, dx. \\
c(p, q) &=& \Frac{1}{\lambda} \int_\Omega p \, q \, dx. \\
V &=& \{ {\bf v} \in (H^1(\Omega))^2; \
{\bf v} = 0 \ {\rm on}\ \Gamma_{left} \cup \Gamma_{bottom} \}
\end{eqnarray*}
When $\lambda$ becomes large, we obtain the incompressible
elasticity problem, that coincides with the Stokes problem.
\subsubsection*{Approximation}
\apindex{P1}%
\apindex{P2}%
As for the Stokes problem,
the \citet*{hood-taylor-73} finite element approximation
is considered.
We introduce a mesh ${\cal T}_h$ of $\Omega$
and the following finite dimensional spaces:
\begin{eqnarray*}
X_h &=& \{ {\bf v} \in (H^1(\Omega)); \
{\bf v}_{/K} \in (P_2)^2, \
\forall K \in {\cal T}_h \}, \\
V_h(\alpha) &=& X_h \cap V, \\
Q_h &=& \{ q \in L^2(\Omega))\cap C^0(\bar{\Omega}); \
q_{/K} \in P_1, \
\forall K \in {\cal T}_h \},
\end{eqnarray*}
The approximate problem writes:
\ \ \ $(VFE)_h$ {\it find ${\bf u}_h \in V_h(1)$ and $p \in Q_h$ such that:}
\[
\begin{array}{lcccl}
a({\bf u}_h,{\bf v}) &+& b({\bf v}, p_h) &=& 0, \ \forall {\bf v}\in V_h(0), \\
b({\bf u}_h,q) &-& c(p,q) &=& 0, \ \forall q \in Q_h.
\end{array}
\]
% --------------------------------------
\myexamplelicense{incompressible-elasticity.cc}
% --------------------------------------
% ======================
\subsubsection*{Comments}
% ======================
\cindex{method!conjugate gradient algorithm}%
\cindex{preconditioner!for nearly incompressible elasticity}%
The problem admits the following matrix form:
\[
\left( \begin{array}{cc}
{\tt a} & {\tt trans(b)} \\
{\tt b} & -{\tt c}
\end{array} \right)
\left( \begin{array}{c}
{\tt uh} \\
{\tt ph}
\end{array} \right)
=
\left( \begin{array}{c}
{\tt lh} \\
{\tt 0}
\end{array} \right)
\]
The problem is similar to the Stokes one (see page~\pageref{ref-pcg-abtb}).
This system is solved by:
\clindex{problem_mixed}%
\begin{lstlisting}[numbers=none,frame=none]
problem_mixed elasticity (a, b, c);
elasticity.solve (lh, field(Qh,0), uh, ph);
\end{lstlisting}
For two-dimensional problems, a direct solver is used by default.
In the three-dimensional case, an iterative algorithm is the default:
the preconditioned conjugate gradient.
The preconditioner is here the mass matrix \code{mp} for the
pressure. As showed by~\citet{Kla-1998}, the number of iterations need by
the conjugate gradient algorithm to reach a given precision
is then independent of the mesh size
and is uniformly bounded when $\lambda$ becomes small,
i.e. in the incompressible case.
% ======================
\subsubsection*{How to run the program}
% ======================
\begin{figure}[htb]
%\begin{center}
\mbox{}\hspace{-0.5cm}
\begin{tabular}{ccc}
\includegraphics[height=7cm]{incompressible-elasticity-square-fig.pdf} &
\includegraphics[height=7.5cm]{incompressible-elasticity-cube-fig.png} &
\stereoglasses
\end{tabular}
%\end{center}
\caption{The incompressible linear elasticity ($\lambda=+\infty$) for $N=2$ and $N=3$.}
\label{fig-incompressible-elasticity-deformation}
\end{figure}
We assume that the previous code is contained in
the file \reffile{incompressible-elasticity.cc}.
Compile the program as usual (see page~\pageref{makefile}):
\begin{verbatim}
make incompressible-elasticity
\end{verbatim}
and enter the commands:
\pindexopt{field}{-scale}%
\begin{verbatim}
mkgeo_grid -t 10 > square.geo
./incompressible-elasticity square.geo 0 > square.field
field square.field -nofill
mkgeo_grid -T 10 > cube.geo
./incompressible-elasticity cube.geo 0 > cube.field
field cube.field -fill -scale 2
\end{verbatim}
The visualization is performed as usual: see section~\ref{sec-howtorun-elasticity},
page~\pageref{sec-howtorun-elasticity}.
Compare the results on {\sc Fig}.~\ref{fig-incompressible-elasticity-deformation},
obtained for $\lambda=+\infty$ with those of
{\sc Fig}.~\ref{fig-embankment-deformation}, page~\pageref{fig-embankment-deformation},
obtained for $\lambda=1$.
Finally, the stress computation and the mesh adaptation loop is left as an exercise to the reader.
|
