manual-solidmechanicsmodel.tex
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Wed, Dec 11, 20:18

manual-solidmechanicsmodel.tex
View Options

	\chapter{Solid Mechanics Model\index{SolidMechanicsModel}\label{sect:smm}}

	The solid mechanics model is a specific implementation of the
	\code{Model} interface dedicated to handle the equations of motion or
	equations of equilibrium. The model is created for a given mesh. It
	will create its own \code{FEEngine} object to compute the interpolation,
	gradient, integration and assembly operations. A
	\code{SolidMechanicsModel} object can simply be created like this:
	\begin{cpp}
	SolidMechanicsModel model(mesh);
	\end{cpp}
	where \code{mesh} is the mesh for which the equations are to be
	solved. A second parameter called \code{spatial\_dimension} can be
	added after \code{mesh} if the spatial dimension of the problem is
	different than that of the mesh.

	This model contains at least the the following six \code{Arrays}:
	\begin{description}
	\item[blocked\_dofs] contains a Boolean value for each degree of
	freedom specifying whether that degree is blocked or not. A
	Dirichlet boundary condition can be prescribed by setting the
	\textbf{blocked\_dofs} value of a degree of freedom to \code{true}.
	A Neumann boundary condition can be applied by setting the
	\textbf{blocked\_dofs} value of a degree of freedom to \code{false}.
	The \textbf{displacement}, \textbf{velocity} and
	\textbf{acceleration} are computed for all degrees of freedom for
	which the \textbf{blocked\_dofs} value is set to \code{false}. For
	the remaining degrees of freedom, the imposed values (zero by
	default after initialization) are kept.
	\item[displacement] contains the displacements of all degrees of
	freedom. It can be either a computed displacement for free degrees
	of freedom or an imposed displacement in case of blocked ones
	($\vec{u}$ in the following).
	\item[velocity] contains the velocities of all degrees of freedom. As
	\textbf{displacement}, it contains computed or imposed velocities
	depending on the nature of the degrees of freedom ($\dot{\vec{u}}$
	in the following).
	\item[acceleration] contains the accelerations of all degrees of
	freedom. As \textbf{displacement}, it contains computed or imposed
	accelerations depending on the nature of the degrees of freedom
	($\ddot{\vec{u}}$ in the following).
	\item[force] contains the external forces applied on the nodes
	($\vec{f}_{\st{ext}}$ in the following).
	\item[residual] contains the difference between external and internal
	forces. On blocked degrees of freedom, \textbf{residual} contains
	the support reactions. ($\vec{r}$ in the following). It should be
	mentioned that at equilibrium \textbf{residual} should be zero on
	free degrees of freedom.
	\end{description}

	Some examples to help to understand how to use this model will be
	presented in the next sections.

	\section{Model Setup}


	\subsection{Setting Initial Conditions \label{sect:smm:initial_condition}}

	For a unique solution of the equations of motion, initial
	displacements and velocities for all degrees of freedom must be
	specified:
	\begin{eqnarray}
	\vec{u}(t=0) = \vec{u}_0\\
	\dot{\vec u}(t=0) =\vec{v}_0
	\end{eqnarray} The solid mechanics model can be initialized as
	follows:
	\begin{cpp}
	model.initFull()
	\end{cpp}
	This function initializes the internal arrays and sets them to
	zero. Initial displacements and velocities that are not equal to zero
	can be prescribed by running a loop over the total number of
	nodes. Here, the initial displacement in $x$-direction and the
	initial velocity in $y$-direction for all nodes is set to $0.1$ and $1$,
	respectively.
	\begin{cpp}
	Array<Real> & disp = model.getDisplacement();
	Array<Real> & velo = model.getVelocity();
	for (UInt i = 0; i < mesh.getNbNodes(); ++i) {
	disp(i, 0) = 0.1;
	velo(i, 1) = 1.;
	}
	\end{cpp}

	\subsection{Setting Boundary Conditions\label{sect:smm:boundary}}
	This section explains how to impose Dirichlet or Neumann boundary
	conditions. A Dirichlet boundary condition specifies the values that
	the displacement needs to take for every point $x$ at the boundary
	($\Gamma_u$) of the problem domain (Fig.~\ref{fig:smm:boundaries}):
	\begin{equation}
	\vec{u} = \bar{\vec u} \quad \forall \vec{x}\in
	\Gamma_{u}
	\end{equation}
	A Neumann boundary condition imposes the value of the gradient of the
	solution at the boundary $\Gamma_t$ of the problem domain
	(Fig.~\ref{fig:smm:boundaries}):
	\begin{equation}
	\vec{t} = \mat{\sigma} \vec{n} = \bar{\vec t} \quad
	\forall \vec{x}\in \Gamma_{t}
	\end{equation}

	\begin{figure} \centering
	\def\svgwidth{0.5\columnwidth}
	\input{figures/problemDomain.pdf_tex}
	\caption{Problem domain $\Omega$ with boundary in three
	dimensions. The Dirchelet and the Neumann regions of the boundary
	are denoted with $\Gamma_u$ and $\Gamma_t$,
	respecitvely.\label{fig:smm:boundaries}}
	\label{fig:problemDomain}
	\end{figure}

	Different ways of imposing these boundary conditions exist. A basic
	way is to loop over nodes or elements at the boundary and apply local
	values. A more advanced method consists of using the notion of the
	boundary of the mesh. In the following both ways are presented.

	Starting with the basic approach, as mentioned, the Dirichlet boundary
	conditions can be applied by looping over the nodes and assigning the
	required values. Figure~\ref{fig:smm:dirichlet_bc} shows a beam with a
	fixed support on the left side. On the right end of the beam, a load
	is applied. At the fixed support, the displacement has a given
	value. For this example, the displacements in both the $x$ and the
	$y$-direction are set to zero. Implementing this displacement boundary
	condition is similar to the implementation of initial displacement
	conditions described above. However, in order to impose a displacement
	boundary condition for all time steps, the corresponding nodes need to
	be marked as boundary nodes as shown in the following code:

	\begin{cpp}
	Array<bool> & blocked = model.getBlockedDOFs();
	const Array<Real> & pos = mesh.getNodes();

	UInt nb_nodes = mesh.getNbNodes();

	for (UInt i = 0; i < nb_nodes; ++i) {
	if(Math::are_float_equal(pos(i, 0), 0)) {
	blocked(i, 0) = true; //block displacement in x-direction
	blocked(i, 1) = true; //block displacement in y-direction
	disp(i, 0) = 0.; //fixed displacement in x-direction
	disp(i, 1)= 0.; //fixed displacement in y-direction
	}
	}
	\end{cpp}

	\begin{figure}[!htb]
	\centering
	\includegraphics[scale=0.4]{figures/dirichlet}
	\caption{Beam with fixed support.\label{fig:smm:dirichlet_bc}}
	\end{figure}

	For the more advanced approach, one needs the notion of a boundary in
	the mesh. Therefore, the boundary should be created before boundary
	condition functors can be applied. Generally the boundary can be
	specified from the mesh file or the geometry. For the first case, the
	function \code{createGroupsFromMeshData} is called. This function
	can read any types of mesh data which are provided in the mesh
	file. If the mesh file is created with Gmsh, the function takes one
	input strings which is either \code{tag\_0}, \code{tag\_1} or
	\code{physical\_names}. The first two tags are assigned by Gmsh to
	each element which shows the physical group that they belong to. In
	Gmsh, it is also possible to consider strings for different groups of
	elements. These elements can be separated by giving a string
	\code{physical\_names} to the function
	\code{createGroupsFromMeshData}. Boundary conditions can also be
	created from the geometry by calling
	\code{createBoundaryGroupFromGeometry}. This function gathers all the
	elements on the boundary of the geometry.

	To apply the required boundary conditions, the function \code{applyBC}
	needs to be called on a \code{SolidMechanicsModel}. This function
	gets a Dirichlet or Neumann functor and a string which specifies the
	desired boundary on which the boundary conditions is to be
	applied. The functors specify the type of conditions to apply. Three
	built-in functors for Dirichlet exist: \code{FlagOnly, FixedValue,}
	and \code{IncrementValue}. The functor \code{FlagOnly} is used if a
	point is fixed in a given direction. Therefore, the input parameter to
	this functor is only the fixed direction. The \code{FixedValue}
	functor is used when a displacement value is applied in a fixed
	direction. The \code{IncrementValue} applies an increment to the
	displacement in a given direction. The following code shows the
	utilization of three functors for the top, bottom and side surface of
	the mesh which were already defined in the Gmsh file:

	\begin{cpp}
	mesh.createGroupsFromMeshData<std::string>("physical_names");

	model.applyBC(BC::Dirichlet::FixedValue(13.0, BC::_y), "Top");

	model.applyBC(BC::Dirichlet::FlagOnly(BC::_x), "Bottom");

	model.applyBC(BC::Dirichlet::IncrementValue(13.0, BC::_x), "Side");
	\end{cpp}

	To apply a Neumann boundary condition, the applied traction or stress
	should be specified before. In case of specifying the traction on the
	surface, the functor \code{FromTraction} of Neumann boundary
	conditions is called. Otherwise, the functor \code{FromStress} should
	be called which gets the stress tensor as an input parameter.

	\begin{cpp}
	Array<Real> surface_traction(3);
	surface_traction(0)=0.0;
	surface_traction(1)=0.0;
	surface_traction(2)=-1.0;

	Array<Real> surface_stress(3, 3, 0.0);
	surface_stress(0,0)=0.0;
	surface_stress(1,1)=0.0;
	surface_stress(2,2)=-1.0;

	model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom");

	model.applyBC(BC::Neumann::FromStress(surface_stress), "Top");
	\end{cpp}

	If the boundary conditions need to be removed during the simulation, a
	functor is called from the Neumann boundary condition to free those
	boundary conditions from the desired boundary.

	\begin{cpp}
	model.applyBC(BC::Neumann::FreeBoundary(), "Side");
	\end{cpp}

	User specified functors can also be implemented. A full example for
	setting both initial and boundary conditions can be found in
	\shellcode{\examplesdir/boundary\_conditions.cc}. The problem solved
	in this example is shown in Fig.~\ref{fig:smm:bc_and_ic}. It consists
	of a plate that is fixed with movable supports on the left and bottom
	side. On the right side, a traction, which increases linearly with the
	number of time steps, is applied. The initial displacement and
	velocity in $x$-direction at all free nodes is zero and two
	respectively.
	\begin{figure}[!htb]
	\centering
	\includegraphics[scale=0.8]{figures/bc_and_ic_example}
	\caption{Plate on movable supports.\label{fig:smm:bc_and_ic}}
	\end{figure}

	\subsection{Material Selector\label{sect:smm:materialselector}}

	If the user wants to assign different materials to the different
	finite elements of choice in \akantu, a material selector has to be
	used. By default, \akantu assigns the first valid material in the
	material file to all elements present in the model (regular continuum
	materials are assigned to the regular elements and cohesive materials
	are assigned to cohesive elements or element facets).

	To assign different materials to specific elements, mesh data
	information such as tag information or specified pyhsical names can be
	used. \code{MeshDataMaterialSelector} class uses this information to
	assign different materials. With the proper physical name or tag name
	and index, different materials can be assigned as demonstrated in the
	examples below.

	\begin{cpp}
	MeshDataMaterialSelector<std::string> * mat_selector;
	mat_selector = new MeshDataMaterialSelector<std::string>("physical_names", model);
	model.setMaterialSelector(*mat_selector);
	\end{cpp}

	In this example the physical names specified in a GMSH geometry file will by
	used to match the material names in the input file.

	Another example would be:

	\begin{cpp}
	MeshDataMaterialSelector<UInt> * mat_selector;
	mat_selector = new MeshDataMaterialSelector<UInt>("tag_1", model);
	model.setMaterialSelector(*mat_selector);
	\end{cpp}

	where \code{tag\_1} of the mesh file is used as the classifier index
	and the elements that have index value equal to one will be assigned
	to the second material in the material file.

	There are four different material selectors pre-defined in
	\akantu. \code{MaterialSelector} and \code{DefaultMaterialSelector} is
	used to assign a material to regular elements by default. For the
	regular elements, as in the example above,
	\code{MeshDataMaterialSelector} can be used to assign different
	materials to different elements. However, for cohesive elements, it is
	not possible to assign multiple cohesive materials automatically and
	user has to write a custom class. \akantu has a pre-defined material
	selector to assign the first cohesive material by default to the
	cohesive elements which is called
	\code{DefaultMaterialCohesiveSelector} and it inherits its properties
	from \code{DefaultMaterialSelector}. This material selector first
	checks if an element is a cohesive element or a facet. If there is,
	then the first cohesive material in the material file is assigned to
	that element. Otherwise, this material is assigned to the element
	facet, to which, a cohesive element might be inserted later. Hence,
	this material selector works for the extrinsic cohesive elements which
	are dynamically inserted throughout the analysis. If the element of
	interest is not a cohesive element or an element facet, then
	\code{DefaultMaterialSelector} is used by default.

	Apart from the \akantu's default material selectors, users can always
	develop their own classes in the main code to tackle various
	multi-material assignment situations.

	% An application of \code{DefaultMaterialCohesiveSelector} and usage in
	% a customly generated material selector class can be seen in
	% \shellcode{\examplesdir/cohesive\_element/cohesive\_extrinsic\_IG\_TG/cohesive\_extrinsic\_IG\_TG.cc}.

	\IfFileExists{manual-cohesive_elements_insertion.tex}{\input{manual-cohesive_elements_insertion}}{}

	\section{Static Analysis\label{sect:smm:static}}

	The \code{SolidMechanicsModel} class can handle different analysis
	methods, the first one being presented is the static case. In this
	case, the equation to solve is
	\begin{equation}
	\label{eqn:smm:static} \mat{K} \vec{u} =
	\vec{f}_{\st{ext}}
	\end{equation}
	where $\mat{K}$ is the global stiffness matrix, $\vec{u}$ the
	displacement vector and $\vec{f}_{\st{ext}}$ the vector of external
	forces applied to the system.

	To solve such a problem, the static solver of the
	\code{SolidMechanicsModel}\index{SolidMechanicsModel} object is used.
	First, a model has to be created and initialized. To create the
	model, a mesh (which can be read from a file) is needed, as explained
	in Section~\ref{sect:common:mesh}. Once an instance of a
	\code{SolidMechanicsModel} is obtained, the easiest way to initialize
	it is to use the \code{initFull}\index{SolidMechanicsModel!initFull}
	method by giving the \code{SolidMechanicsModelOptions}. These options
	specify the type of analysis to be performed and whether the materials
	should be initialized or not.
	\begin{cpp}
	SolidMechanicsModel model(mesh);
	model.initFull(SolidMechanicsModelOptions(_static));
	\end{cpp}
	Here, a static analysis is chosen by passing the argument
	\code{\_static} to the method. By default, the Boolean for no
	initialization of the materials is set to false, so that they are
	initialized during the \code{initFull}. The method \code{initFull}
	also initializes all appropriate vectors to zero. Once the model is
	created and initialized, the boundary conditions can be set as
	explained in Section~\ref{sect:smm:boundary}. Boundary conditions
	will prescribe the external forces for some free degrees of freedom
	$\vec{f}_{\st{ext}}$ and displacements for some others. At this point
	of the analysis, the function
	\code{solveStep}\index{SolidMechanicsModel!solveStep} can be called:
	\begin{cpp}
	model.solveStep<_scm_newton_raphson_tangent_modified, _scc_residual>(1e-4, 1);
	\end{cpp}
	This function is templated by the solving method and the convergence
	criterion and takes two arguments: the tolerance and the maximum
	number of iterations, which are $\num{1e-4}$ and $1$ for this example. The
	modified Newton-Raphson method is chosen to solve the system. In this
	method, the equilibrium equation (\ref{eqn:smm:static}) is modified in
	order to apply a Newton-Raphson convergence algorithm:
	\begin{align}\label{eqn:smm:static-newton-raphson}
	\mat{K}^{i+1}\delta\vec{u}^{i+1} &= \vec{r} \\
	&= \vec{f}_{\st{ext}} -\vec{f}_{\st{int}}\\
	&= \vec{f}_{\st{ext}} - \mat{K}^{i} \vec{u}^{i}\\
	\vec{u}^{i+1} &= \vec{u}^{i} + \delta\vec{u}^{i+1}~,\nonumber
	\end{align}
	where $\delta\vec{u}$ is the increment of displacement to be added
	from one iteration to the other, and $i$ is the Newton-Raphson
	iteration counter. By invoking the \code{solveStep} method in the
	first step, the global stiffness matrix $\mat{K}$ from
	Equation~(\ref{eqn:smm:static}) is automatically assembled. A
	Newton-Raphson iteration is subsequently started, $\mat{K}$ is updated
	according to the displacement computed at the previous iteration and
	one loops until the forces are balanced (\code{\_scc\_residual}), \ie
	$\|\|\vec{r}\|\| < \mbox{\code{\_scc\_residual}}$. One can also iterate
	until the increment of displacement is zero (\code{\_scc\_increment})
	which also means that the equilibrium is found. For a linear elastic
	problem, the solution is obtained in one iteration and therefore the
	maximum number of iterations can be set to one. But for a non-linear
	case, one needs to iterate as long as the norm of the residual exceeds
	the tolerance threshold and therefore the maximum number of iterations
	has to be higher, e.g. $100$:
	\begin{cpp}
	model.solveStep<_scm_newton_raphson_tangent_modified,_scc_residual>(1e-4, 100)
	\end{cpp}
	At the end of the analysis, the final solution is stored in the
	\textbf{displacement} vector. A full example of how to solve a static
	problem is presented in the code \code{\examplesdir/static/static.cc}.
	This example is composed of a 2D plate of steel, blocked with rollers
	on the left and bottom sides as shown in Figure \ref{fig:smm:static}.
	The nodes from the right side of the sample are displaced by $0.01\%$
	of the length of the plate.

	\begin{figure}[!htb]
	\centering
	\includegraphics[scale=1.05]{figures/static}
	\caption{Numerical setup\label{fig:smm:static}}
	\end{figure}

	The results of this analysis is depicted in
	Figure~\ref{fig:smm:implicit:static_solution}.

	\begin{figure}[!htb]
	\centering
	\includegraphics[width=.7\linewidth]{figures/static_analysis}
	\caption{Solution of the static analysis. Left: the initial
	condition, right: the solution (deformation magnified 50 times)}
	\label{fig:smm:implicit:static_solution}
	\end{figure}

	\section{Dynamic Methods} \label{sect:smm:Dynamic_methods}

	Different ways to solve the equations of motion are implemented in the
	solid mechanics model. The complete equations that should be solved
	are:
	\begin{equation}
	\label{eqn:equation-motion}
	\mat{M}\ddot{\vec{u}} +
	\mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = \vec{f}_{\st{ext}}~,
	\end{equation}
	where $\mat{M}$, $\mat{C}$ and $\mat{K}$ are the mass,
	damping and stiffness matrices, respectively.

	In the previous section, it has already been discussed how to solve this
	equation in the static case, where $\ddot{\vec{u}} = \dot{\vec{u}} = 0$. Here
	the method to solve this equation in the general case will be presented. For
	this purpose, a time discretization has to be specified. The most common
	discretization method in solid mechanics is the Newmark-$\beta$ method, which is
	also the default in \akantu.

	For the Newmark-$\beta$ method, (\ref{eqn:equation-motion}) becomes a
	system of three equations (see \cite{curnier92a} \cite{hughes-83a} for
	more details):
	\begin{align}
	\mat{M} \ddot{\vec{u}}_{n+1} + \mat{C}\dot{\vec{u}}_{n+1} + \mat{K} \vec{u}_{n+1} &={\vec{f}_{\st{ext}}}_{\, n+1}
	\label{eqn:equation-motion-discret} \\
	\vec{u}_{n+1} &=\vec{u}_{n} + \left(1 - \alpha\right) \Delta t \dot{\vec{u}}_{n} +
	\alpha \Delta t \dot{\vec{u}}_{n+1} + \left(\frac{1}{2} -
	\alpha\right) \Delta t^2
	\ddot{\vec{u}}_{n} \label{eqn:finite-difference-1}\\
	\dot{\vec{u}}_{n+1} &= \dot{\vec{u}}_{n} + \left(1 - \beta\right)
	\Delta t \ddot{\vec{u}}_{n} + \beta \Delta t
	\ddot{\vec{u}}_{n+1} \label{eqn:finite-difference-2}
	\end{align}

	In these new equations, $\ddot{\vec{u}}_{n}$, $\dot{\vec{u}}_{n}$ and
	$\vec{u}_{n}$ are the approximations of $\ddot{\vec{u}}(t_n)$,
	$\dot{\vec{u}}(t_n)$ and $\vec{u}(t_n)$.
	Equation~(\ref{eqn:equation-motion-discret}) is the equation of motion
	discretized in space (finite-element discretization), and equations
	(\ref{eqn:finite-difference-1}) and (\ref{eqn:finite-difference-2})
	are discretized in both space and time (Newmark discretization). The
	$\alpha$ and $\beta$ parameters determine the stability and the
	accuracy of the algorithm. Classical values for $\alpha$ and $\beta$
	are usually $\beta = 1/2$ for no numerical damping and $0 < \alpha <
	1/2$.

	\begin{center}
	\begin{tabular}{cll}
	\toprule
	$\alpha$ & Method ($\beta = 1/2$) & Type\\
	\midrule
	$0$ & central difference & explicit\\
	$1/6$ & Fox-Goodwin (royal road) &implicit\\
	$1/3$ & Linear acceleration &implicit\\
	$1/2$ & Average acceleration (trapezoidal rule)& implicit\\
	\bottomrule
	\end{tabular}
	\end{center}

	The solution of this system of equations,
	(\ref{eqn:equation-motion-discret})-(\ref{eqn:finite-difference-2}) is
	split into a predictor and a corrector system of equations. Moreover,
	in the case of a non-linear equations, an iterative algorithm such as
	the Newton-Raphson method is applied. The system of equations can be
	written as:

	\begin{enumerate}
	\item \textit{Predictor:}
	\begin{align}
	\vec{u}_{n+1}^{0} &= \vec{u}_{n} + \Delta t
	\dot{\vec{u}}_{n} + \frac{\Delta t^2}{2} \ddot{\vec{u}}_{n} \\
	\dot{\vec{u}}_{n+1}^{0} &= \dot{\vec{u}}_{n} + \Delta t
	\ddot{\vec{u}}_{n} \\
	\ddot{\vec{u}}_{n+1}^{0} &= \ddot{\vec{u}}_{n}
	\end{align}

	\item \textit{Solve:}
	\begin{align}
	\left(c \mat{M} + d \mat{C} + e \mat{K}_{n+1}^i\right)
	\vec{w} = {\vec{f}_{\st{ext}}}_{\,n+1} - {\vec{f}_{\st{int}}}_{\,n+1}^i -
	\mat{C} \dot{\vec{u}}_{n+1}^i - \mat{M} \ddot{\vec{u}}_{n+1}^i = \vec{r}_{n+1}^i
	\end{align}

	\item \textit{Corrector:}
	\begin{align}
	\ddot{\vec{u}}_{n+1}^{i+1} &= \ddot{\vec{u}}_{n+1}^{i} +c \vec{w} \\
	\dot{\vec{u}}_{n+1}^{i+1} &= \dot{\vec{u}}_{n+1}^{i} + d\vec{w} \\
	\vec{u}_{n+1}^{i+1} &= \vec{u}_{n+1}^{i} + e \vec{w}
	\end{align}
	\end{enumerate}

	where $i$ is the Newton-Raphson iteration counter and $c$, $d$ and $e$
	are parameters depending on the method used to solve the equations

	\begin{center}
	\begin{tabular}{lcccc}
	\toprule
	& $\vec{w}$ & $e$ & $d$ & $c$\\
	\midrule
	in acceleration &$ \delta\ddot{\vec{u}}$ & $\alpha \beta\Delta t^2$ &$\beta \Delta t$ &$1$\\
	in velocity & $ \delta\dot{\vec{u}}$& $\frac{1}{\beta} \Delta t$ & $1$ & $\alpha\Delta t$\\
	in displacement &$\delta\vec{u}$ & $ 1$ & $\frac{1}{\alpha} \Delta t$ & $\frac{1}{\alpha \beta} \Delta t^2$\\
	\bottomrule
	\end{tabular}
	\end{center}

	%\note{If you want to use the implicit solver \akantu should be compiled at least with one sparse matrix solver such as Mumps\cite{mumps}.}


	\subsection{Implicit Time Integration}
	To solve a problem with an implicit time integration scheme, first a
	\code{SolidMechanicsModel} object has to be created and initialized.
	Then the initial and boundary conditions have to be set. Everything
	is similar to the example in the static case
	(Section~\ref{sect:smm:static}), however, in this case the implicit
	dynamic scheme is selected at the initialization of the model.

	\begin{cpp}
	SolidMechanicsModel model(mesh);
	model.initFull(SolidMechanicsModelOptions(_implicit_dynamic));
	/Boundary conditions see Section~%\ref{sect:smm:boundary}% /
	\end{cpp}
	Because a dynamic simulation is conducted, an integration time step
	$\Delta t$ has to be specified. In the case of implicit simulations,
	\akantu implements a trapezoidal rule by default. That is to say
	$\alpha = 1/2$ and $\beta = 1/2$ which is unconditionally
	stable. Therefore the value of the time step can be chosen arbitrarily
	within reason. \index{SolidMechanicsModel!setTimeStep}
	\begin{cpp}
	model.setTimeStep(time_step);
	\end{cpp}
	Since the system has to be solved for a given amount of time, the
	method \code{solveStep()}, (which has already been used in the static
	example in Section~\ref{sect:smm:static}), is called inside a time
	loop:
	\begin{cpp}
	/// time loop
	Real time = 0.;
	for (UInt s = 1; time <max_time; ++s, time += time_step) {
	model.solveStep<_scm_newton_raphson_tangent_modified,_scc_increment>(1e-12, 100);
	}
	\end{cpp}
	An example of solid mechanics with an implicit time integration scheme
	is presented in
	\shellcode{\examplesdir/implicit/implicit\_dynamic.cc}. This example
	consists of a 3D beam of
	$\SI{10}{\metre}\,\times\,\SI{1}{\metre}\,\times\,\SI{1}{\metre}$
	blocked on one side and is on a roller on the other side. A constant
	force of \SI{5}{\kilo\newton} is applied in its middle.
	Figure~\ref{fig:smm:implicit:dynamic} presents the geometry of this
	case. The material used is a fictitious linear elastic material with a
	density of \SI{1000}{\kilo\gram\per\cubic\metre}, a Young's Modulus of
	\SI{120}{\mega\pascal} and Poisson's ratio of $0.3$. These values
	were chosen to simplify the analytical solution.

	An approximation of the dynamic response of the middle point of the
	beam is given by:
	\begin{equation}
	\label{eqn:smm:implicit}
	u\left(\frac{L}{2}, t\right)
	= \frac{1}{\pi^4} \left(1 - cos\left(\pi^2 t\right) +
	\frac{1}{81}\left(1 - cos\left(3^2 \pi^2 t\right)\right) +
	\frac{1}{625}\left(1 - cos\left(5^2 \pi^2 t\right)\right)\right)
	\end{equation}

	\begin{figure}[!htb]
	\centering
	\includegraphics[scale=.6]{figures/implicit_dynamic}
	\caption{Numerical setup}
	\label{fig:smm:implicit:dynamic}
	\end{figure}

	Figure \ref{fig:smm:implicit:dynamic_solution} presents the deformed
	beam at 3 different times during the simulation: time steps 0, 1000 and
	2000.

	\begin{figure}[!htb]
	\centering
	\setlength{\unitlength}{0.1\textwidth}
	\begin{tikzpicture}
	\node[above right] (img) at (0,0)
	{\includegraphics[width=.6\linewidth]{figures/dynamic_analysis}};
	\node[left] at (0pt,20pt) {$0$}; \node[left] at (0pt,60pt) {$1000$};
	\node[left] at (0pt,100pt) {$2000$};
	\end{tikzpicture}

	\caption{Deformed beam at 3 different times (displacement are
	magnified by a factor 10).}
	\label{fig:smm:implicit:dynamic_solution}
	\end{figure}

	\subsection{Explicit Time Integration}
	\label{ssect:smm:expl-time-integr}

	The explicit dynamic time integration scheme is based on the
	Newmark-$\beta$ scheme with $\alpha=0$ (see equations
	\ref{eqn:equation-motion-discret}-\ref{eqn:finite-difference-2}). In
	\akantu, $\beta$ is defaults to $\beta=1/2$, see section
	\ref{sect:smm:Dynamic_methods}.

	The initialization of the simulation is similar to the static and
	implicit dynamic version. The model is created from the
	\code{SolidMechanicsModel} class. In the initialization, the explicit
	scheme is selected using the \code{\_explicit\_lumped\_mass} constant.

	\begin{cpp}
	SolidMechanicsModel model(mesh);
	model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass));
	\end{cpp}
	\index{SolidMechanicsModel!initFull}
	\note{Writing \code{model.initFull()} or \code{model.initFull(SolidMechanicsModelOptions());} is
	equivalent to use the \code{\_explicit\_lumped\_mass} keyword, as this
	is the default case.}

	The explicit time integration scheme implemented in \akantu uses a
	lumped mass matrix $\mat{M}$ (reducing the computational cost). This
	matrix is assembled by distributing the mass of each element onto its
	nodes. The resulting $\mat{M}$ is therefore a diagonal matrix stored
	in the \textbf{mass} vector of the model.

	The explicit integration scheme is conditionally stable. The time step
	has to be smaller than the stable time step which is obtained in
	\akantu as follows:

	\begin{cpp}
	time_step = model.getStableTimeStep();
	\end{cpp} \index{SolidMechanicsModel!StableTimeStep}

	The stable time step is defined as:
	\begin{equation}
	\label{eqn:smm:explicit:stabletime}
	\Delta t_{\st{crit}} = \Delta x \sqrt{\frac{\rho}{2 \mu + \lambda}}
	\end{equation}
	where $\Delta x$ is a characteristic length (\eg the inradius in the
	case of linear triangle element), $\mu$ and $\lambda$ are the first
	and second Lame's coefficients and $\rho$ is the density. The stable
	time step corresponds to the time the fastest wave (the compressive
	wave) needs to travel the characteristic length of the mesh.
	However, it is recommended to impose a time step that is smaller than the stable time step, for
	instance, by multiplying the stable time step by a safety factor
	smaller than one.

	\begin{cpp}
	const Real safety_time_factor = 0.8;
	Real applied_time_step = time_step * safety_time_factor;
	model.setTimeStep(applied_time_step);
	\end{cpp}
	\index{SolidMechanicsModel!setTimeStep} The initial displacement and
	velocity fields are, by default, equal to zero if not given
	specifically by the user (see \ref{sect:smm:initial_condition}).

	Like in implicit dynamics, a time loop is used in which the
	displacement, velocity and acceleration fields are updated at each
	time step. The values of these fields are obtained from the
	Newmark$-\beta$ equations with $\beta=1/2$ and $\alpha=0$. In \akantu
	these computations at each time step are invoked by calling the
	function \code{solveStep}:
	\begin{cpp}
	for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) {
	model.solveStep();
	}
	\end{cpp} \index{SolidMechanicsModel!solveStep}
	The method
	\code{solveStep} wraps the four following functions:
	\begin{itemize}
	\item \code{model.explicitPred()} allows to compute the displacement
	field at $t+1$ and a part of the velocity field at $t+1$, denoted by
	$\vec{\dot{u}^{\st{p}}}_{n+1}$, which will be used later in the method
	\code{model.explicitCorr()}. The equations are:
	\begin{align}
	\vec{u}_{n+1} &= \vec{u}_{n} + \Delta t
	\vec{\dot{u}}_{n} + \frac{\Delta t^2}{2} \vec{\ddot{u}}_{n}\\
	\vec{\dot{u}^{\st{p}}}_{n+1} &= \vec{\dot{u}}_{n} + \Delta t
	\vec{\ddot{u}}_{n}
	\label{eqn:smm:explicit:onehalfvelocity}
	\end{align}

	\item \code{model.updateResidual()} and
	\code{model.updateAcceleration()} compute the acceleration increment
	$\delta \vec{\ddot{u}}$:
	\begin{equation}
	\left(\mat{M} + \frac{1}{2} \Delta t \mat{C}\right)
	\delta \vec{\ddot{u}} = \vec{f_{\st{ext}}} - \vec{f}_{\st{int}\, n+1}
	- \mat{C} \vec{\dot{u}}_{n} - \mat{M} \vec{\ddot{u}}_{n}
	\end{equation}

	\note{The internal force $\vec{f}_{\st{int}\, n+1}$ is computed from
	the displacement $\vec{u}_{n+1}$ based on the constitutive law.}

	\item \code{model.explicitCorr()} computes the velocity and
	acceleration fields at $t+1$:
	\begin{align}
	\vec{\dot{u}}_{n+1} &= \vec{\dot{u}^{\st{p}}}_{n+1} + \frac{\Delta t}{2}
	\delta \vec{\ddot{u}} \\ \vec{\ddot{u}}_{n+1} &=
	\vec{\ddot{u}}_{n} + \delta \vec{\ddot{u}}
	\end{align}
	\end{itemize}

	The use of an explicit time integration scheme is illustrated by the
	example:\par
	\noindent \shellcode{\examplesdir/explicit/explicit\_dynamic.cc}\par
	\noindent This example models the propagation of a wave in a steel beam. The
	beam and the applied displacement in the $x$ direction are shown in
	Figure~\ref{fig:smm:explicit}.

	\begin{figure}[!htb] \centering
	\begin{tikzpicture}
	\coordinate (c) at (0,2);
	\draw[shift={(c)},thick, color=blue] plot [id=x, domain=-5:5, samples=50] ({\x, {(40 * sin(0.1pi3\x) exp(- (0.1pi3\x)(0.1pi3*\x) / 4))}});
	\draw[shift={(c)},-latex] (-6,0) -- (6,0) node[right, below] {$x$};
	\draw[shift={(c)},-latex] (0,-0.7) -- (0,1) node[right] {$u$};
	\draw[shift={(c)}] (-0.1,0.6) node[left] {$A$}-- (1.5,0.6);

	\coordinate (l) at (0,0.6);
	\draw[shift={(0,-0.7)}] (-5, 0) -- (5,0) -- (5, 1) -- (-5, 1) -- cycle;
	\draw[shift={(l)}, latex-latex] (-5,0)-- (5,0) node [midway, above] {$L$};
	\draw[shift={(l)}] (5,0.2)-- (5,-0.2);
	\draw[shift={(l)}] (-5,0.2)-- (-5,-0.2);

	\coordinate (h) at (5.3,-0.7);
	\draw[shift={(h)}, latex-latex] (0,0)-- (0,1) node [midway, right] {$h$};
	\draw[shift={(h)}] (-0.2,1)-- (0.2,1);
	\draw[shift={(h)}] (-0.2,0)-- (0.2,0);
	\end{tikzpicture}

	\caption{Numerical setup \label{fig:smm:explicit}}
	\end{figure}

	The length and height of the beam are $L=\SI{10}{\metre}$ and
	$h = \SI{1}{\metre}$, respectively. The material is linear elastic,
	homogeneous and isotropic (density:
	\SI{7800}{\kilo\gram\per\cubic\metre}, Young's modulus:
	\SI{210}{\giga\pascal} and Poisson's ratio: $0.3$). The imposed
	displacement follow a Gaussian function with a maximum amplitude of $A = \SI{0.01}{\meter}$. The
	potential, kinetic and total energies are computed. The safety factor
	is equal to $0.8$.

	\input{manual-constitutive-laws}

	\section{Adding a New Constitutive Law}\index{Material!create a new
	material}

	There are several constitutive laws in \akantu as described in the
	previous Section~\ref{sect:smm:CL}. It is also possible to use a
	user-defined material for the simulation. These materials are referred
	to as local materials since they are local to the example of the user
	and not part of the \akantu library. To define a new local material,
	two files (\code {material\_XXX.hh} and \code{material\_XXX.cc}) have
	to be provided where \code{XXX} is the name of the new material. The
	header file \code {material\_XXX.hh} defines the interface of your
	custom material. Its implementation is provided in the
	\code{material\_XXX.cc}. The new law must inherit from the
	\code{Material} class or any other existing material class. It is
	therefore necessary to include the interface of the parent material
	in the header file of your local material and indicate the inheritance
	in the declaration of the class:
	\begin{cpp}
	/* ---------------------------------------------------------------------- */
	#include "material.hh"
	/* ---------------------------------------------------------------------- */

	#ifndef __AKANTU_MATERIAL_XXX_HH__
	#define __AKANTU_MATERIAL_XXX_HH__

	__BEGIN_AKANTU__

	class MaterialXXX : public Material {

	/// declare here the interface of your material

	};
	\end{cpp}
	In the header file the user also needs to declare all the members of the new
	material. These include the parameters that a read from the
	material input file, as well as any other material parameters that will be computed during the simulation and internal variables.\\
	In the following the example of adding a new damage material will be presented. In this
	case the parameters in the material will consist of the Young's modulus, the
	Poisson coefficient, the resistance to damage and the damage threshold. The
	material will then from these values compute its Lam\'{e} coefficients and its
	bulk modulus. Furthermore, the user has to add a new internal variable
	\code{damage} in order to store the amount of damage at each quadrature point in
	each step of the simulation. For this specific material the member declaration
	inside the class will look as follows:
	\begin{cpp}
	class LocalMaterialDamage : public Material {

	/// declare constructors/destructors here

	/// declare methods and accessors here

	/* -------------------------------------------------------------------- */
	/* Class Members */
	/* -------------------------------------------------------------------- */

	AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real);
	private:

	/// the young modulus
	Real E;

	/// Poisson coefficient
	Real nu;

	/// First Lame coefficient
	Real lambda;

	/// Second Lame coefficient (shear modulus)
	Real mu;

	/// resistance to damage
	Real Yd;

	/// damage threshold
	Real Sd;

	/// Bulk modulus
	Real kpa;

	/// damage internal variable
	InternalField<Real> damage;

	};
	\end{cpp}
	In order to enable to print the material parameters at any point in
	the user's example file using the standard output stream by typing:
	\begin{cpp}
	for (UInt m = 0; m < model.getNbMaterials(); ++m)
	std::cout << model.getMaterial(m) << std::endl;
	\end{cpp}
	the standard output stream operator has to be redefined. This should be done at the end of the header file:
	\begin{cpp}
	class LocalMaterialDamage : public Material {

	/// declare here the interace of your material

	}:
	/* ---------------------------------------------------------------------- */
	/* inline functions */
	/* ---------------------------------------------------------------------- */
	/// standard output stream operator
	inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this)
	{
	_this.printself(stream);
	return stream;
	}
	\end{cpp}
	However, the user still needs to register the material parameters that
	should be printed out. The registration is done during the call of the
	constructor. Like all definitions the implementation of the
	constructor has to be written in the \code{material\_XXX.cc}
	file. However, the declaration has to be provided in the
	\code{material\_XXX.hh} file:
	\begin{cpp}
	class LocalMaterialDamage : public Material {
	/* -------------------------------------------------------------------- */
	/* Constructors/Destructors */
	/* -------------------------------------------------------------------- */
	public:

	LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = "");
	};
	\end{cpp}
	The user can now define the implementation of the constructor in the
	\code{material\_XXX.cc} file:
	\begin{cpp}
	/* ---------------------------------------------------------------------- */
	#include "local_material_damage.hh"
	#include "solid_mechanics_model.hh"

	__BEGIN_AKANTU__

	/* ---------------------------------------------------------------------- */
	LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model,
	const ID & id) :
	Material(model, id),
	damage("damage", *this) {
	AKANTU_DEBUG_IN();

	this->registerParam("E", E, 0., _pat_parsable, "Young's modulus");
	this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio");
	this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient");
	this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient");
	this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient");
	this->registerParam("Yd", Yd, 50., _pat_parsmod);
	this->registerParam("Sd", Sd, 5000., _pat_parsmod);

	damage.initialize(1);

	AKANTU_DEBUG_OUT();
	}
	\end{cpp}
	During the intializer list the reference to the model and the material id are
	assigned and the constructor of the internal field is called. Inside the scope
	of the constructor the internal values have to be initialized and the
	parameters, that should be printed out, are registered with the function:
	\code{registerParam}\index{Material!registerParam}:
	\begin{cpp}
	void registerParam(name of the parameter (key in the material file),
	member variable,
	default value (optional parameter),
	access permissions,
	description);
	\end{cpp}
	The available access permissions are as follows:
	\begin{itemize}
	\item \code{\_pat\_internal}: Parameter can only be output when the material is printed.
	\item \code{\_pat\_writable}: User can write into the parameter. The parameter is output when the material is printed.
	\item \code{\_pat\_readable}: User can read the parameter. The parameter is output when the material is printed.
	\item \code{\_pat\_modifiable}: Parameter is writable and readable.
	\item \code{\_pat\_parsable}: Parameter can be parsed, \textit{i.e.} read from the input file.
	\item \code{\_pat\_parsmod}: Parameter is modifiable and parsable.
	\end{itemize}

	In order to implement the new constitutive law the user needs to
	specify how the additional material parameters, that are not
	defined in the input material file, should be calculated. Furthermore,
	it has to be defined how stresses and the stable time step should be
	computed for the new local material. In the case of implicit
	simulations, in addition, the computation of the tangent stiffness needs
	to be defined. Therefore, the user needs to redefine the following
	functions of the parent material:
	\begin{cpp}
	void initMaterial();

	// for explicit and implicit simulations void
	computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);

	// for implicit simulations
	void computeTangentStiffness(const ElementType & el_type,
	Array<Real> & tangent_matrix,
	GhostType ghost_type = _not_ghost);

	// for explicit and implicit simulations
	Real getStableTimeStep(Real h, const Element & element);
	\end{cpp}
	In the following a detailed description of these functions is provided:
	\begin{itemize}

	\item \code{initMaterial}:~ This method is called after the material
	file is fully read and the elements corresponding to each material
	are assigned. Some of the frequently used constant parameters are
	calculated in this method. For example, the Lam\'{e} constants of
	elastic materials can be considered as such parameters.

	\item \code{computeStress}:~ In this method, the stresses are
	computed based on the constitutive law as a function of the strains of the
	quadrature points. For example, the stresses for the elastic
	material are calculated based on the following formula:
	\begin{equation}
	\label{eqn:smm:constitutive_elastic}
	\mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon}
	\end{equation}

	Therefore, this method contains a loop on all quadrature points
	assigned to the material using the two macros:\par
	\code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_BEGIN}\par
	\code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_END}

	\begin{cpp}
	MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type);

	// sigma <- f(grad_u)

	MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
	\end{cpp}

	\note{The strain vector in \akantu contains the values of $\nabla
	\vec{u}$, i.e. it is really the \emph{displacement gradient},}

	\item \code{computeTangentStiffness}:~ This method is called when
	the tangent to the stress-strain curve is desired (see Fig \ref
	{fig:smm:AL:K}). For example, it is called in the implicit solver
	when the stiffness matrix for the regular elements is assembled
	based on the following formula:
	\begin{equation}
	\label{eqn:smm:constitutive_elasc} \mat{K }
	=\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}}
	\end{equation}

	Therefore, in this method, the \code{tangent} matrix (\mat{D}) is
	computed for a given strain.

	\note{ The \code{tangent} matrix is a $4^{th}$ order tensor which is
	stored as a matrix in Voigt notation.}

	\begin{figure}[!htb]
	\begin{center}
	\includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf}
	\caption{Tangent to the stress-strain curve.}
	\label{fig:smm:AL:K}
	\end{center}
	\end{figure}


	\item \code{getStableTimeStep}:~ The stable time step should be
	calculated based on \eqref{eqn:smm:explicit:stabletime} for the
	conditionally stable methods. This function depends on the longitudinal wave
	speed which changes for different materials and strains. Therefore,
	the value of this velocity should be defined in this method for
	each new material.

	\note {In case of need, the calculation of the longitudinal and
	shear wave speeds can be done in separate functions
	(\code{getPushWaveSpeed} and \code{getShearWaveSpeed}).
	Therefore, the first function can be called for calculation of the
	stable time step.}
	\end{itemize}
	Once the declaration and implementation of the new material has been
	completed, this material can be used in the user's example by including the header file:
	\begin{cpp}
	#include "material_XXX.hh"
	\end{cpp}
	For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by
	default, the materials are initialized inside the method
	\code{initFull}. If a local material should be used instead, the
	initialization of the material has to be postponed until the local
	material is registered in the model. Therefore, the model is
	initialized with the boolean for skipping the material initialization
	equal to true:
	\begin{cpp}
	/// model initialization
	model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass,true));
	\end{cpp}
	Once the model has been initialized, the local material needs
	to be registered in the model:
	\begin{cpp}
	model.registerNewCustomMaterials<XXX>("name_of_local_material");
	\end{cpp}
	Only at this point the material can be initialized:
	\begin{cpp}
	model.initMaterials();
	\end{cpp}
	A full example for adding a new damage law can be found in
	\shellcode{\examplesdir/new\_material}.

	%%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End:

manual-solidmechanicsmodel.texNo OneTemporaryActions

File Metadata

manual-solidmechanicsmodel.texView Options

Event Timeline

manual-solidmechanicsmodel.tex
No OneTemporary
Actions

manual-solidmechanicsmodel.tex
View Options