file by <a href="step-18.data/intro.pdf">clicking here</a>]
</p>
+
<P>
This tutorial program is another one in the series on the elasticity problem
that we have already started with step-8 and step-17. It extends it into two
WIDTH="47" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
SRC="step-18.data/intro/img27.png"
ALT="$ \vec f(\vec x,t)$">
-.
+. The changes in
+configuration can therefore be considered as being stationary
+instantaneously. An alternative view of this is that <IMG
+ WIDTH="10" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
+ SRC="step-18.data/intro/img28.png"
+ ALT="$ t$">
+ is not really a time
+variable, but only a time-like parameter that governs the evolution of the
+problem.
<P>
While these equations are sufficient to describe small deformations, computing
large deformations is a little more complicated. To do so, let us first
-introduce a stress variable <IMG
+introduce a tensorial stress variable <IMG
WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img28.png"
+ SRC="step-18.data/intro/img29.png"
ALT="$ \sigma$">
-, and write the differential equations in
-terms of the stress:
+, and write the differential
+equations in terms of the stress:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
WIDTH="54" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img29.png"
+ SRC="step-18.data/intro/img30.png"
ALT="$\displaystyle - \div\sigma$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
WIDTH="29" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">in <IMG
WIDTH="34" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img30.png"
+ SRC="step-18.data/intro/img31.png"
ALT="$ \Omega(t)$">
<IMG
WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
-->
<IMG
WIDTH="86" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img31.png"
+ SRC="step-18.data/intro/img32.png"
ALT="$ \Gamma_D\subset\partial\Omega(t)$">
<IMG
WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
-->
<IMG
WIDTH="116" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img32.png"
+ SRC="step-18.data/intro/img33.png"
ALT="$ \Gamma_N=\partial\Omega(t)\backslash\Gamma_D$">
<IMG
WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
<BR CLEAR="ALL"><P></P>
Note that these equations are posed on a domain <IMG
WIDTH="34" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img30.png"
+ SRC="step-18.data/intro/img31.png"
ALT="$ \Omega(t)$">
that
changes with time, with the boundary moving according to the
-->
<IMG
WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img33.png"
+ SRC="step-18.data/intro/img34.png"
ALT="$ \vec u(\vec x,t)$">
of the points on the boundary. To
-complete this system, we have to specify the relationship between the
-stress and the strain, as follows:
+complete this system, we have to specify the incremental relationship between
+the stress and the strain, as follows:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="82" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img34.png"
+ SRC="step-18.data/intro/img35.png"
ALT="$\displaystyle \dot\sigma = C \varepsilon (\dot{\vec u}),$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<A NAME="eq:stress-strain">(11)</A></TD></TR>
<BR CLEAR="ALL"><P></P>
where a dot indicates a time derivative. Both the stress <IMG
WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img28.png"
+ SRC="step-18.data/intro/img29.png"
ALT="$ \sigma$">
and the
strain <!-- MATH
-->
<IMG
WIDTH="34" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img35.png"
+ SRC="step-18.data/intro/img36.png"
ALT="$ \varepsilon(\vec u)$">
are symmetric tensors of rank 2.
the time component using a backward Euler scheme. This leads to a
discrete equilibrium of force at time step <IMG
WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img36.png"
+ SRC="step-18.data/intro/img37.png"
ALT="$ n$">
:
<P></P>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="106" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img37.png"
+ SRC="step-18.data/intro/img38.png"
ALT="$\displaystyle -\div\sigma^n = f^n,$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(12)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="167" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img38.png"
+ SRC="step-18.data/intro/img39.png"
ALT="$\displaystyle \sigma^n = \sigma^{n-1} + C \varepsilon (\Delta \vec u^n),$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(13)</TD></TR>
-->
<IMG
WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img39.png"
+ SRC="step-18.data/intro/img40.png"
ALT="$ \Delta \vec u^n$">
the incremental displacement for time step
<IMG
WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img36.png"
+ SRC="step-18.data/intro/img37.png"
ALT="$ n$">
-. This way, if we want to solve for the displacement increment, we
+. In addition, we have to specify initial data <!-- MATH
+ $\vec u(\cdot,0)=\vec u_0$
+ -->
+<IMG
+ WIDTH="85" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-18.data/intro/img41.png"
+ ALT="$ \vec u(\cdot,0)=\vec u_0$">
+.
+This way, if we want to solve for the displacement increment, we
have to solve the following system:
<P></P>
<DIV ALIGN="CENTER"><TABLE CELLPADDING="0" WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
WIDTH="108" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img40.png"
+ SRC="step-18.data/intro/img42.png"
ALT="$\displaystyle - \div C \varepsilon(\Delta\vec u^n)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
WIDTH="110" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img41.png"
+ SRC="step-18.data/intro/img43.png"
ALT="$\displaystyle = \vec f + \div\sigma^{n-1}$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">in <!-- MATH
-->
<IMG
WIDTH="59" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img42.png"
+ SRC="step-18.data/intro/img44.png"
ALT="$ \Omega(t_{n-1})$">
<IMG
WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
WIDTH="71" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img43.png"
+ SRC="step-18.data/intro/img45.png"
ALT="$\displaystyle \Delta \vec u^n(\vec x,t)$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
WIDTH="196" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img44.png"
+ SRC="step-18.data/intro/img46.png"
ALT="$\displaystyle = \vec d(\vec x,t_n) - \vec d(\vec x,t_{n-1}) \qquad$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">on <!-- MATH
-->
<IMG
WIDTH="111" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img45.png"
+ SRC="step-18.data/intro/img47.png"
ALT="$ \Gamma_D\subset\partial\Omega(t_{n-1})$">
<IMG
WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="RIGHT"><IMG
WIDTH="119" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img46.png"
+ SRC="step-18.data/intro/img48.png"
ALT="$\displaystyle \vec n \ C \varepsilon(\Delta \vec u^n(\vec x,t))$"></TD>
<TD NOWRAP ALIGN="LEFT"><IMG
WIDTH="196" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img47.png"
+ SRC="step-18.data/intro/img49.png"
ALT="$\displaystyle = \vec b(\vec x,t_n)-\vec b(\vec x,t_{n-1}) \qquad$"></TD>
<TD> </TD>
<TD NOWRAP ALIGN="LEFT">on <!-- MATH
-->
<IMG
WIDTH="141" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img48.png"
+ SRC="step-18.data/intro/img50.png"
ALT="$ \Gamma_N=\partial\Omega(t_{n-1})\backslash\Gamma_D$">
<IMG
WIDTH="9" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
-->
<IMG
WIDTH="394" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img49.png"
+ SRC="step-18.data/intro/img51.png"
ALT="$ \Delta \vec u^n \in
\{v\in H^1(\Omega(t_{n-1}))^d: v\vert _{\Gamma_D}=\vec d(\cdot,t_n) - \vec d(\cdot,t_{n-1})\}$">
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="412" HEIGHT="78" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img50.png"
+ SRC="step-18.data/intro/img52.png"
ALT="\begin{gather*}\begin{split}(C \varepsilon(\Delta\vec u^n), \varepsilon(\varphi)...
...in H^1(\Omega(t_{n-1}))^d: \vec v\vert _{\Gamma_D}=0\}. \end{split}\end{gather*}"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
-->
<IMG
WIDTH="44" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img51.png"
+ SRC="step-18.data/intro/img53.png"
ALT="$ \vec b = 0$">
, and that the deformation of the
body is driven by body forces <IMG
and prescribed boundary displacements
<IMG
WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img52.png"
+ SRC="step-18.data/intro/img54.png"
ALT="$ \vec d$">
alone. It is also worth noting that when integrating by parts, we
would get terms of the form <!-- MATH
-->
<IMG
WIDTH="158" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img53.png"
+ SRC="step-18.data/intro/img55.png"
ALT="$ (C \varepsilon(\Delta\vec u^n), \nabla \varphi
)_{\Omega(t_{n-1})}$">
, but that we replace it with the term involving the
-->
<IMG
WIDTH="35" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img54.png"
+ SRC="step-18.data/intro/img56.png"
ALT="$ \varepsilon(\varphi)$">
instead of <!-- MATH
$\nabla\varphi$
-->
<IMG
WIDTH="28" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img55.png"
+ SRC="step-18.data/intro/img57.png"
ALT="$ \nabla\varphi$">
. Due to
the symmetry of <IMG
<P>
The system at time step <IMG
WIDTH="14" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img36.png"
+ SRC="step-18.data/intro/img37.png"
ALT="$ n$">
, to be solved on the old domain
<!-- MATH
-->
<IMG
WIDTH="59" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img42.png"
+ SRC="step-18.data/intro/img44.png"
ALT="$ \Omega(t_{n-1})$">
, has exactly the form of a stationary elastic
problem, and is therefore similar to what we have already implemented
There are differences, however:
<OL>
-<LI>We have to move the mesh after each time step, in order to be
+<LI>We have to move (update) the mesh after each time step, in order to be
able to solve the next time step on a new domain;
<P>
-->
<IMG
WIDTH="39" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img56.png"
+ SRC="step-18.data/intro/img58.png"
ALT="$ \sigma^{n-1}$">
to compute the next incremental
displacement, i.e. we need to compute it at the end of the time step
<P>
As indicated above, we need to have the stress variable <IMG
WIDTH="23" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img57.png"
+ SRC="step-18.data/intro/img59.png"
ALT="$ \sigma^n$">
available
when computing time step <IMG
WIDTH="41" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img58.png"
+ SRC="step-18.data/intro/img60.png"
ALT="$ n+1$">
, and we can compute it using
<P></P>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="167" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img59.png"
+ SRC="step-18.data/intro/img61.png"
ALT="$\displaystyle \sigma^n = \sigma^{n-1} + C \varepsilon (\Delta \vec u^n).$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<A NAME="eq:stress-update">(18)</A></TD></TR>
There are, despite the apparent simplicity of this equation, two questions
that we need to discuss. The first concerns the way we store <IMG
WIDTH="23" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img57.png"
+ SRC="step-18.data/intro/img59.png"
ALT="$ \sigma^n$">
: even
if we compute the incremental updates <!-- MATH
-->
<IMG
WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img39.png"
+ SRC="step-18.data/intro/img40.png"
ALT="$ \Delta \vec u^n$">
using lowest-order
finite elements, then its symmetric gradient <!-- MATH
-->
<IMG
WIDTH="56" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img60.png"
+ SRC="step-18.data/intro/img62.png"
ALT="$ \varepsilon(\Delta\vec u^n)$">
is
in general still a function that is not easy to describe. In particular, it is
it is not even a bi- or trilinear function. Thus, it
is a priori not clear how to store <IMG
WIDTH="23" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img57.png"
+ SRC="step-18.data/intro/img59.png"
ALT="$ \sigma^n$">
in a computer program.
-->
<IMG
WIDTH="135" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img61.png"
+ SRC="step-18.data/intro/img63.png"
ALT="$ (\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})}$">
. In practice, we of
course replace this term by numerical quadrature:
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="518" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img62.png"
+ SRC="step-18.data/intro/img64.png"
ALT="$\displaystyle (\sigma^{n-1},\varepsilon(\varphi))_{\Omega(t_{n-1})} = \sum_{K\s...
...athbb{T}}} \sum_q w_q \ \sigma^{n-1}(\vec x_q) : \varepsilon(\varphi(\vec x_q),$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<BR CLEAR="ALL"><P></P>
where <IMG
WIDTH="23" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img63.png"
+ SRC="step-18.data/intro/img65.png"
ALT="$ w_q$">
are the quadrature weights and <IMG
WIDTH="21" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img64.png"
+ SRC="step-18.data/intro/img66.png"
ALT="$ \vec x_q$">
the quadrature points on
cell <IMG
WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img65.png"
+ SRC="step-18.data/intro/img67.png"
ALT="$ K$">
. This should make clear that what we really need is not the stress
<!-- MATH
-->
<IMG
WIDTH="39" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img56.png"
+ SRC="step-18.data/intro/img58.png"
ALT="$ \sigma^{n-1}$">
in itself, but only the values of the stress in the quadrature
points on all cells. This, however, is a simpler task: we only have to provide
-->
<IMG
WIDTH="85" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img66.png"
+ SRC="step-18.data/intro/img68.png"
ALT="$ \varepsilon(\Delta \vec u^n(\vec x_q))$">
, multiply it by the stress-strain tensor <IMG
WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
-->
<IMG
WIDTH="52" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img67.png"
+ SRC="step-18.data/intro/img69.png"
ALT="$ \sigma^n(\vec x_q)$">
at quadrature point <IMG
WIDTH="12" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img68.png"
+ SRC="step-18.data/intro/img70.png"
ALT="$ q$">
.
-->
<IMG
WIDTH="36" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img69.png"
+ SRC="step-18.data/intro/img71.png"
ALT="$ \Delta u^n$">
on the domain <!-- MATH
$\Omega(t_{n-1})$
-->
<IMG
WIDTH="59" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img42.png"
+ SRC="step-18.data/intro/img44.png"
ALT="$ \Omega(t_{n-1})$">
,
and then use this displacement increment to both update the stress as well as
-->
<IMG
WIDTH="43" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img70.png"
+ SRC="step-18.data/intro/img72.png"
ALT="$ \Omega(t_n)$">
on which the next increment
is computed. What we have to make sure, in this context, is that moving the
-->
<IMG
WIDTH="28" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img71.png"
+ SRC="step-18.data/intro/img73.png"
ALT="$ \Delta\vec u$">
can be decomposed into three parts, a linear translation (the constant part
-of the displacement field in the neighborhood of a point), a dilational
+of the displacement increment field in the neighborhood of a point), a
+dilational
component (that part of the gradient of the displacement field that has a
nonzero divergence), and a rotation. A linear translation of the material does
not affect the stresses that are frozen into it - the stress values are
-->
<IMG
WIDTH="108" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img72.png"
+ SRC="step-18.data/intro/img74.png"
ALT="$ \Delta\vec u=(y, -x)^T$">
, with which <!-- MATH
$\varepsilon(\Delta \vec u)=0$
-->
<IMG
WIDTH="77" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img73.png"
+ SRC="step-18.data/intro/img75.png"
ALT="$ \varepsilon(\Delta \vec u)=0$">
). Nevertheless, if the the material was pre-stressed in a certain
direction, then this direction will be rotated along with the material. To
-->
<IMG
WIDTH="61" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img74.png"
+ SRC="step-18.data/intro/img76.png"
ALT="$ R(\Delta \vec u^n)$">
that
describes, in each point the rotation due to the displacement increments. It
is not hard to see that the actual dependence of <IMG
WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img75.png"
+ SRC="step-18.data/intro/img77.png"
ALT="$ R$">
on <!-- MATH
$\Delta \vec u^n$
-->
<IMG
WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img39.png"
+ SRC="step-18.data/intro/img40.png"
ALT="$ \Delta \vec u^n$">
can
only be through the curl of the displacement, rather than the displacement
the increment describe translations, its divergence the dilational modes, and
the curl the rotational modes). Since the exact form of <IMG
WIDTH="17" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img75.png"
+ SRC="step-18.data/intro/img77.png"
ALT="$ R$">
is cumbersome, we
only state it in the program code, and note that the correct updating formula
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="299" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img76.png"
+ SRC="step-18.data/intro/img78.png"
ALT="$\displaystyle \sigma^n = R(\Delta \vec u^n)^T [\sigma^{n-1} + C \varepsilon (\Delta \vec u^n)] R(\Delta \vec u^n).$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
<A NAME="eq:stress-update+rot">(20)</A></TD></TR>
<LI><TT>assemble_system ()</TT> [via <TT>solve_timestep ()</TT>]:
This first function is also the most interesting one. It assembles the
linear system corresponding to the discretized version of equation
- (<A HREF="#eq:linear-system">17</A>). This leads to a system matrix <!-- MATH
+ (<A HREF="#eq:linear-system"><IMG ALIGN="BOTTOM" BORDER="1" ALT="[*]"
+ SRC="step-18.data/intro/file:/home/bangerth/bin/share/lib/latex2html/icons/crossref.png"></A>). This leads to a system matrix <!-- MATH
$A_{ij} = \sum_K
A^K_{ij}$
-->
<IMG
WIDTH="105" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img77.png"
+ SRC="step-18.data/intro/img79.png"
ALT="$ A_{ij} = \sum_K
A^K_{ij}$">
built up of local contributions on each cell <IMG
WIDTH="19" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img65.png"
+ SRC="step-18.data/intro/img67.png"
ALT="$ K$">
with entries
<P></P>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="171" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img78.png"
+ SRC="step-18.data/intro/img80.png"
ALT="$\displaystyle A^K_{ij} = (C \varepsilon(\varphi_j), \varepsilon(\varphi_i))_K;$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(21)</TD></TR>
<BR CLEAR="ALL"><P></P>
In practice, <IMG
WIDTH="29" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img79.png"
+ SRC="step-18.data/intro/img81.png"
ALT="$ A^K$">
is computed using numerical quadrature according to the
formula
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="277" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img80.png"
+ SRC="step-18.data/intro/img82.png"
ALT="$\displaystyle A^K_{ij} = \sum_q w_q [\varepsilon(\varphi_i(\vec x_q)) : C : \varepsilon(\varphi_j(\vec x_q))],$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(22)</TD></TR>
<BR CLEAR="ALL"><P></P>
with quadrature points <IMG
WIDTH="21" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img64.png"
+ SRC="step-18.data/intro/img66.png"
ALT="$ \vec x_q$">
and weights <IMG
WIDTH="23" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img63.png"
+ SRC="step-18.data/intro/img65.png"
ALT="$ w_q$">
. We have built these
contributions before, in step-8 and step-17, but in both of these cases we
-->
<IMG
WIDTH="84" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img81.png"
+ SRC="step-18.data/intro/img83.png"
ALT="$ \varepsilon(\varphi_i),\varepsilon(\varphi_j)$">
. This is not really
convenient, in particular if we want to consider more complicated elasticity
-->
<IMG
WIDTH="241" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img82.png"
+ SRC="step-18.data/intro/img84.png"
ALT="$ C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl}
+ \delta_{il} \delta_{jk})$">
. While we in fact do not use a more complicated
as well as to create a symmetric tensor of rank 2 (the strain
tensor) from the gradients of a shape function <IMG
WIDTH="20" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img83.png"
+ SRC="step-18.data/intro/img85.png"
ALT="$ \varphi_i$">
at a quadrature
point <IMG
WIDTH="21" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img64.png"
+ SRC="step-18.data/intro/img66.png"
ALT="$ \vec x_q$">
on a given cell. At the top of the implementation of this
example program, you will find such functions. The first one,
<TT>get_stress_strain_tensor</TT>, takes two arguments corresponding to
the Lamé constants <IMG
WIDTH="14" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img84.png"
+ SRC="step-18.data/intro/img86.png"
ALT="$ \lambda$">
and <IMG
WIDTH="14" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img85.png"
+ SRC="step-18.data/intro/img87.png"
ALT="$ \mu$">
and returns the stress-strain tensor
for the isotropic case corresponding to these constants (in the program, we
<TT>get_strain</TT> takes an object of type <TT>FEValues</TT> and indices
<IMG
WIDTH="10" HEIGHT="17" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img86.png"
+ SRC="step-18.data/intro/img88.png"
ALT="$ i$">
and <IMG
WIDTH="12" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img68.png"
+ SRC="step-18.data/intro/img70.png"
ALT="$ q$">
and returns the symmetric gradient, i.e. the strain,
corresponding to shape function <!-- MATH
-->
<IMG
WIDTH="49" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img87.png"
+ SRC="step-18.data/intro/img89.png"
ALT="$ \varphi_i(\vec x_q)$">
, evaluated on the cell
on which the <TT>FEValues</TT> object was last reinitialized.
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="331" HEIGHT="71" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img88.png"
+ SRC="step-18.data/intro/img90.png"
ALT="\begin{gather*}\begin{split}f^K_i &= (\vec f, \varphi_i)_K -(\sigma^{n-1},\varep...
...gma^{n-1}_q : \varepsilon(\varphi_i(\vec x_q)) \right\} \end{split}\end{gather*}"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(23)</TD></TR>
</TABLE></DIV>
<BR CLEAR="ALL"><P></P>
-to the right hand side of (<A HREF="#eq:linear-system">17</A>) is equally
+to the right hand side of (<A HREF="#eq:linear-system"><IMG ALIGN="BOTTOM" BORDER="1" ALT="[*]"
+ SRC="step-18.data/intro/file:/home/bangerth/bin/share/lib/latex2html/icons/crossref.png"></A>) is equally
straightforward (note that we do not consider any boundary tractions <IMG
WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img89.png"
+ SRC="step-18.data/intro/img91.png"
ALT="$ \vec b$">
here). Remember that we only had to store the old stress in the
quadrature points of cells. In the program, we will provide a variable
-->
<IMG
WIDTH="39" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img90.png"
+ SRC="step-18.data/intro/img92.png"
ALT="$ \sigma^{n-1}_q$">
in each quadrature point. With this the code for the right
hand side looks as this, again rather elegant:
-->
<IMG
WIDTH="97" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img91.png"
+ SRC="step-18.data/intro/img93.png"
ALT="$ \vec f(\vec x_q) \cdot \varphi_i(\vec x_q)$">
, we have made use of the fact that for the chosen finite element, only
one vector component (namely <TT>component_i</TT>) of <IMG
WIDTH="20" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img83.png"
+ SRC="step-18.data/intro/img85.png"
ALT="$ \varphi_i$">
is
nonzero, and that we therefore also have to consider only one component of
-->
<IMG
WIDTH="41" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img92.png"
+ SRC="step-18.data/intro/img94.png"
ALT="$ \vec f(\vec x_q)$">
.
-->
<IMG
WIDTH="37" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img39.png"
+ SRC="step-18.data/intro/img40.png"
ALT="$ \Delta \vec u^n$">
computed before, we update the stress values in all quadrature points
- according to (<A HREF="#eq:stress-update">18</A>) and (<A HREF="#eq:stress-update+rot">20</A>),
+ according to (<A HREF="#eq:stress-update"><IMG ALIGN="BOTTOM" BORDER="1" ALT="[*]"
+ SRC="step-18.data/intro/file:/home/bangerth/bin/share/lib/latex2html/icons/crossref.png"></A>) and (<A HREF="#eq:stress-update+rot"><IMG ALIGN="BOTTOM" BORDER="1" ALT="[*]"
+ SRC="step-18.data/intro/file:/home/bangerth/bin/share/lib/latex2html/icons/crossref.png"></A>),
including the rotation of the coordinate system.
<P>
</H4> The most obvious extension is to use a more
realistic material model for large-scale quasistatic deformation. The natural
choice for this would be plasticity, in which a nonlinear relationship between
-stress and strain replaces equation (<A HREF="#eq:stress-strain">11</A>). Plasticity
+stress and strain replaces equation (<A HREF="#eq:stress-strain"><IMG ALIGN="BOTTOM" BORDER="1" ALT="[*]"
+ SRC="step-18.data/intro/file:/home/bangerth/bin/share/lib/latex2html/icons/crossref.png"></A>). Plasticity
models are usually rather complicated to program since the stress-strain
dependence is generally non-smooth. The material can be thought of being able
to withstand only a maximal stress (the yield stress) after which it starts to
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="311" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img93.png"
+ SRC="step-18.data/intro/img95.png"
ALT="$\displaystyle E(\vec u) = (\varepsilon(\vec u), C\varepsilon(\vec u))_{\Omega} - (\vec f, \vec u)_{\Omega} - (\vec b, \vec u)_{\Gamma_N},$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(24)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="88" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img94.png"
+ SRC="step-18.data/intro/img96.png"
ALT="$\displaystyle f(\sigma(\vec u)) \le 0$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
(25)</TD></TR>
<TR VALIGN="MIDDLE">
<TD NOWRAP ALIGN="CENTER"><IMG
WIDTH="101" HEIGHT="52" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img95.png"
+ SRC="step-18.data/intro/img97.png"
ALT="$\displaystyle \nu = \frac{\lambda}{2(\lambda+\mu)},$"></TD>
<TD NOWRAP WIDTH="10" ALIGN="RIGHT">
</TD></TR>
-->
<IMG
WIDTH="30" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img96.png"
+ SRC="step-18.data/intro/img98.png"
ALT="$ \lambda,\mu$">
are the Lamé constants of the material.
Physical constraints indicate that <!-- MATH
-->
<IMG
WIDTH="86" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img97.png"
+ SRC="step-18.data/intro/img99.png"
ALT="$ -1\le \nu\le \tfrac 12$">
-. If <IMG
+ (the condition
+also follows from mathematical stability considerations). If <IMG
WIDTH="13" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-18.data/intro/img98.png"
+ SRC="step-18.data/intro/img100.png"
ALT="$ \nu$">
approaches <IMG
WIDTH="15" HEIGHT="35" ALIGN="MIDDLE" BORDER="0"
- SRC="step-18.data/intro/img99.png"
+ SRC="step-18.data/intro/img101.png"
ALT="$ \tfrac 12$">
, then the material becomes incompressible. In that
case, pure displacement-based formulations are no longer appropriate for the