Solid Mechanics Model¶
The solid mechanics model is a specific implementation of the 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 FEEngine
object to compute the
interpolation, gradient, integration and assembly operations. A
SolidMechanicsModel
object can simply
be created like this:
SolidMechanicsModel model(mesh);
where mesh
is the mesh for which the equations are to be
solved. A second parameter called spatial_dimension
can be
added after mesh
if the spatial dimension of the problem is
different than that of the mesh.
This model contains at least the following six Arrays
:
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 blocked_dofs value of a degree of freedom to
true
. A Neumann boundary condition can be applied by setting the blocked_dofs value of a degree of freedom tofalse
. The displacement, velocity and acceleration are computed for all degrees of freedom for which the blocked_dofs value is set tofalse
. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept.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).
velocity
:contains the velocities of all degrees of freedom. As displacement, it contains computed or imposed velocities depending on the nature of the degrees of freedom (\(\dot{\vec{u}}\) in the following).
acceleration
:contains the accelerations of all degrees of freedom. As displacement, it contains computed or imposed accelerations depending on the nature of the degrees of freedom (\(\ddot{\vec{u}}\) in the following).
external_force
:contains the external forces applied on the nodes (\(\vec{f}_{\st{ext}}\) in the following).
internal_force
:contains the internal forces on the nodes (\(\vec{f}_{\mathrm{int}}\) in the following).
Some examples to help to understand how to use this model will be presented in the next sections.
Model Setup¶
Setting Initial Conditions¶
For a unique solution of the equations of motion, initial displacements and velocities for all degrees of freedom must be specified:
The solid mechanics model can be initialized as follows:
model.initFull()
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:
auto & disp = model.getDisplacement();
auto & velo = model.getVelocity();
for (UInt node = 0; node < mesh.getNbNodes(); ++node) {
disp(node, 0) = 0.1;
velo(node, 1) = 1.;
}
Setting Boundary Conditions¶
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:smm:boundaries
):
A Neumann boundary condition imposes the value of the gradient of the
solution at the boundary \(\Gamma_t\) of the problem domain
(fig:smm:boundaries
):
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. Fig. 8 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 using the function blocked
. While,
in order to impose a load on the right side, the nodes are not marked.
The detail codes are shown as follows
auto & blocked = model.getBlockedDOFs();
const auto & pos = mesh.getNodes();
UInt nb_nodes = mesh.getNbNodes();
for (UInt node = 0; node < nb_nodes; ++node) {
if(Math::are_float_equal(pos(node, _x), 0)) {
blocked(node, _x) = true; // block dof in xdirection
blocked(node, _y) = true; // block dof in ydirection
disp(node, _x) = 0.; // fixed displacement in xdirection
disp(node, _y) = 0.; // fixed displacement in ydirection
} else if (Math::are_float_equal(pos(node, _y), 0)) {
blocked(node, _x) = false; // unblock dof in xdirection
forces(node, _x) = 10.; // force in xdirection
}
}
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 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 tag_0
, tag_1
or
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
physical_names
to the function
createGroupsFromMeshData
mesh.createGroupsFromMeshData<std::string>("physical_names").
Boundary conditions support can also be created from the geometry by calling
createBoundaryGroupFromGeometry
. This function gathers all the elements on
the boundary of the geometry.
To apply the required boundary conditions, the function applyBC
needs to be called on a
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 builtin functors for Dirichlet exist:
FlagOnly
, FixedValue
and IncrementValue
. The functor 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 FixedValue
functor is used when a
displacement value is applied in a fixed direction. The 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
model.applyBC(BC::Dirichlet::FixedValue(13.0, _y), "Top");
model.applyBC(BC::Dirichlet::FlagOnly(_x), "Bottom");
model.applyBC(BC::Dirichlet::IncrementValue(13.0, _x), "Side");
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
FromTraction
of Neumann
boundary conditions is called. Otherwise, the functor FromStress
should be called which gets the stress tensor
as an input parameter
Vector<Real> surface_traction{0., 0., 1.};
auto surface_stress(3, 3) = Matrix<Real>::eye(3);
model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom");
model.applyBC(BC::Neumann::FromStress(surface_stress), "Top");
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
model.applyBC(BC::Neumann::FreeBoundary(), "Side");
User specified functors can also be implemented. A full example for
setting both initial and boundary conditions can be found in
examples/boundary_conditions.cc
. The problem solved
in this example is shown in figsmmbc_and_id
. 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.
As it is mentioned in Section ref{sect:common:groups}, node and element groups can be used to assign the boundary conditions. A generic example is given below with a Dirichlet boundary condition:
// create a node group
NodeGroup & node_group = mesh.createNodeGroup("nodes_fix");
/* fill the node group with the nodes you want */
// create an element group using the existing node group
mesh.createElementGroupFromNodeGroup("el_fix",
"nodes_fix",
spatial_dimension1);
// boundary condition can be applied using the element group name
model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "el_fix");
Material Selector¶
If the user wants to assign different materials to different
finite elements groups 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 physical names can be used.
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:
auto mat_selector =
std::make_shared<MeshDataMaterialSelector<std::string>>("physical_names",
model);
model.setMaterialSelector(mat_selector);
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 to use the first (tag_0
) or the second
(tag_1
) tag associated to each elements in the mesh:
auto mat_selector = std::make_shared<MeshDataMaterialSelector<UInt>>(
"tag_1", model, first_index);
model.setMaterialSelector(*mat_selector);
where first_index
(default is 1) is the value of tag_1
that will
be associated to the first material in the material input file. The following
values of the tag will be associated with the following materials.
There are four different material selectors predefined in Akantu
.
MaterialSelector
and
DefaultMaterialSelector
is used
to assign a material to regular elements by default. For the regular elements,
as in the example above, MeshDataMaterialSelector
can be used to assign different materials to
different elements.
Apart from the Akantu
’s default material selectors, users can always
develop their own classes in the main code to tackle various
multimaterial assignment situations.
Static Analysis¶
The SolidMechanicsModel
class can
handle different analysis methods, the first one being presented is the static
case. In this case, the equation to solve is
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
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
SolidMechanicsModel
is obtained, the
easiest way to initialize it is to use the initFull
method by giving the
SolidMechanicsModelOptions
.
These options specify the type of analysis to be performed and whether the
materials should be initialized with initMaterials
or not
SolidMechanicsModel model(mesh);
model.initFull(_analysis_method = _static);
Here, a static analysis is chosen by passing the argument
_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 initFull
. The method 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 Setting Boundary Conditions.
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
solveStep
can be called
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 1);
solver.set("threshold", 1e4);
solver.set("convergence_type", SolveConvergenceCriteria::_residual);
model.solveStep();
This function is templated by the solving method and the convergence criterion and takes two arguments: the tolerance and the maximum number of iterations (100 by default), which are \(10^{4}\) and \(1\) for this example. The modified NewtonRaphson method is chosen to solve the system. In this method, the equilibrium equation (Eq. 3) is modified in order to apply a NewtonRaphson convergence algorithm:
where \(\delta\vec{u}\) is the increment of displacement to be added from
one iteration to the other, and \(i\) is the NewtonRaphson iteration
counter. By invoking the solveStep
method in the first step, the global
stiffness matrix \(\mat{K}\) from (Eq. 3) is automatically
assembled. A NewtonRaphson 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
(SolveConvergenceCriteria::_residual
), i.e. \(\vec{r} <\)
threshold
. One can also
iterate until the increment of displacement is zero
(SolveConvergenceCriteria::_solution
) 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 nonlinear 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\)
solver.set("max_iterations", 100);
model.solveStep();
At the end of the analysis, the final solution is stored in the
displacement vector. A full example of how to solve a static
problem is presented in the code examples/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 Fig. 10.
The nodes from the right side of the sample are displaced by \(0.01\%\)
of the length of the plate.
The results of this analysis is depicted in Fig. 11.
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:
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, (Eq. 4) becomes a system of three equations (see [curnier92a][hughes83a] for more details):
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:equationmotiondiscret}) is the equation of motion discretized in space (finiteelement discretization), and the equations above 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\).
\(alpha\) 
Method (\(beta= 1/2\)) 
Type 

\(0\) 
central difference 
explicit 
\(\frac{1}{6}\) 
FoxGoodwin(royal road) 
implicit 
\(\frac{1}{3}\) 
Linear acceleration 
implicit 
\(\frac{1}{2}\) 
Average acceleration (trapeziodal rule) 
implicit 
The solution of this system of equations, (Eq. 5)) is split into a predictor and a corrector system of equations. Moreover, in the case of a nonlinear equations, an iterative algorithm such as the NewtonRaphson method is applied. The system of equations can be written as:
Predictor:
\[\begin{split}\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{split}\]Solve:
\[\begin{split}\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{split}\]Corrector:
\[\begin{split}\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{split}\]
where \(i\) is the NewtonRaphson iteration counter and \(c\), \(d\) and \(e\) are parameters depending on the method used to solve the equations
\(vec{w}\) 
\(e\) 
\(d\) 
\(c\) 


in acceleration 
\(\delta\ddot{\vec{u}}\) 
\(\alpha\beta\Delta t^2\) 
\(\beta\Delta t\) 
\(1\) 
in velocity 
\(\delta\dot{\vec{u}}\) 
\(\alpha\Delta t\) 
\(1\) 
\(\frac{1}{\beta\Delta t}\) 
in displacement 
\(\delta\vec{u}\) 
\(1\) 
\(\frac{1}{\alpha\Delta t}\) 
\(\frac{1}{\alpha\beta \Delta t^2}\) 
Note
If you want to use the implicit solver Akantu
should be compiled at
least with one sparse matrix solver such as Mumps [mumps].
Implicit Time Integration¶
To solve a problem with an implicit time integration scheme, first a
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:
SolidMechanicsModel model(mesh);
model.initFull(_analysis_method = _implicit_dynamic);
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:
model.setTimeStep(time_step);
Since the system has to be solved for a given amount of time steps, the
method solveStep()
, (which has already been used in the static
example in Section~ref{sect:smm:static}), is called inside a time
loop:
/// time loop
Real time = 0.;
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 100);
solver.set("threshold", 1e12);
solver.set("convergence_type", SolveConvergenceCriteria::_solution);
for (UInt s = 1; time <max_time; ++s, time += time_step) {
model.solveStep();
}
An example of solid mechanics with an implicit time integration scheme is
presented in examples/implicit/implicit_dynamic.cc
. This example consists of
a 3D beam of
\(10\mathrm{m}\times1\mathrm{m}\times1\mathrm{m}\) blocked
on one side and is on a roller on the other side. A constant force of
\(5\mathrm{kN}\) is applied in its middle.
Fig. 12 presents the geometry of this case. The
material used is a fictitious linear elastic material with a density of
\(1000 \mathrm{kg/m}^3\), a Young’s Modulus of
\(120 \mathrm{MPa}\) 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:
Figure Fig. 13 presents the deformed beam at 3 different times during the simulation: time steps 0, 1000 and 2000.
Explicit Time Integration¶
The explicit dynamic time integration scheme is based on the
Newmark\(\beta\) scheme with \(\alpha=0\) (see equations
ref{eqn:equationmotiondiscret}ref{eqn:finitedifference2}). 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
SolidMechanicsModel
class. In the initialization, the explicit
scheme is selected using the _explicit_lumped_mass
constant:
SolidMechanicsModel model(mesh);
model.initFull(_analysis_method = _explicit_lumped_mass);
Note
Writing model.initFull()
or model.initFull();
is
equivalent to use the _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 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:
critical_time_step = model.getStableTimeStep();
The stable time step corresponds to the time the fastest wave (the compressive wave) needs to travel the characteristic length of the mesh:
where \(\Delta x\) is a characteristic length (eg the inradius in the case of linear triangle element) and \(c\) is the celerity of the fastest wave in the material. It is generally the compressive wave of celerity \(c = \sqrt{\frac{2 \mu + \lambda}{\rho}}\), \(\mu\) and \(\lambda\) are the first and second Lame’s coefficients and \(\rho\) is the density. 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:
const Real safety_time_factor = 0.8;
Real applied_time_step = critical_time_step * safety_time_factor;
model.setTimeStep(applied_time_step);
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:math:beta equations with \(\beta=1/2\) and \(\alpha=0\). In Akantu
these computations at each time step are invoked by calling the
function solveStep
:
for (UInt s = 1; (s1)*applied_time_step < total_time; ++s) {
model.solveStep();
}
The method solveStep
wraps the four following functions:
model.explicitPred()
allows to compute the displacementfield 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
model.explicitCorr()
. The equations are:\[\begin{split}\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}\end{split}\]
model.updateResidual()
andmodel.updateAcceleration()
compute the acceleration increment\(\delta \vec{\ddot{u}}\):
\[\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}^{\st{p}}}_{n+1}  \mat{M} \vec{\ddot{u}}_{n}\] The internal force \(\vec{f}_{\st{int}\, n+1}\) is computed from the
displacement \(\vec{u}_{n+1}\) based on the constitutive law.
model.explicitCorr()
computes the velocity andacceleration fields at \(t+1\):
\[\begin{split}\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{split}\]
The use of an explicit time integration scheme is illustrated by the example:
examples/explicit/explicit_dynamic.cc
. 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 Fig. 14.
The length and height of the beam are \(L={10}\textrm{m}\) and \(h = {1}\textrm{m}\), respectively. The material is linear elastic, homogeneous and isotropic (density: SI{7800}{kilogrampercubicmetre}, Young’s modulus: SI{210}{gigapascal} and Poisson’s ratio: \(0.3\)). The imposed displacement follow a Gaussian function with a maximum amplitude of \(A = {0.01}\textrm{m}\). The potential, kinetic and total energies are computed. The safety factor is equal to \(0.8\).
Constitutive Laws¶
In order to compute an element’s response to deformation, one needs to use an appropriate constitutive relationship. The constitutive law is used to compute the element’s stresses from the element’s strains.
In the finiteelement discretization, the constitutive formulation is applied to every quadrature point of each element. When the implicit formulation is used, the tangent matrix has to be computed.
element_material
vector. For
every material assigned to the problem one has to specify the material
characteristics (constitutive behavior and material properties) using
the text input file (see [sect:io:material]).Akantu
provides a special data structure, the at InternalField
. The internal fields are inheriting from the at
ElementTypeMapArray
. Furthermore,
it provides several functions for initialization, autoresizing and auto
removal of quadrature points.Sometimes it is also desired to generate random distributions of internal parameters. An example might be the critical stress at which the material fails. To generate such a field, in the text input file, a random quantity needs be added to the base value:
All parameters are real numbers. For the uniform distribution, minimum and maximum values have to be specified. Random parameters are defined as a \(base\) value to which we add a random number that follows the chosen distribution.
The Uniform distribution is gives a random values between in \([min, max)\). The Weibull distribution is characterized by the following cumulative distribution function:
which depends on \(m\) and \(\lambda\), which are the shape parameter and the scale parameter. These random distributions are different each time the code is executed. In order to obtain always the same one, it possible to manually set the seed that is the number from which these pseudorandom distributions are created. This can be done by adding the following line to the input file outside the material parameters environments:
seed = 1.0
where the value 1.0 can be substituted with any number. Currently
Akantu
can reproduce always the same distribution when the seed is
specified only in serial. The value of the seed can be also
specified directly in the code (for instance in the main file) with the
command:
RandGenerator<Real>::seed(1.0)
The same command, with empty brackets, can be used to check the value of the seed used in the simulation.
The following sections describe the constitutive models implemented in
Akantu
. In Appendix 7 a summary of
the parameters for all materials of Akantu
is provided.
Elasticity¶
The elastic law is a commonly used constitutive relationship that can be used for a wide range of engineering materials (e.g., metals, concrete, rock, wood, glass, rubber, etc.) provided that the strains remain small (i.e., small deformation and stress lower than yield strength).
The elastic laws are often expressed as \(\boldsymbol{\sigma} = \boldsymbol{C}:\boldsymbol{\varepsilon}\) with where \(\boldsymbol{\sigma}\) is the Cauchy stress tensor, \(\boldsymbol{\varepsilon}\) represents the infinitesimal strain tensor and \(\boldsymbol{C}\) is the elastic modulus tensor.
Linear isotropic¶
The linear isotropic elastic behavior is described by Hooke’s law, which states that the stress is linearly proportional to the applied strain (material behaves like an ideal spring), as illustrated in Figure [fig:smm:cl:elastic].
The equation that relates the strains to the displacements is: point) from the displacements as follows:
where \(\boldsymbol{\varepsilon}\) represents the infinitesimal strain tensor, \(\nabla_{0}\boldsymbol{u}\) the displacement gradient tensor according to the initial configuration. The constitutive equation for isotropic homogeneous media can be expressed as:
where \(\boldsymbol{\sigma}\) is the Cauchy stress tensor (\(\lambda\) and \(\mu\) are the the first and second Lame’s coefficients).
In Voigt notation this correspond to
Linear anisotropic¶
This formulation is not sufficient to represent all elastic material behavior. Some materials have characteristic orientation that have to be taken into account. To represent this anisotropy a more general stressstrain law has to be used. For this we define the elastic modulus tensor as follow:
To simplify the writing of input files the \(\boldsymbol{C}\) tensor is expressed in the material basis. And this basis as to be given too. This basis \(\Omega_{{\mathrm{mat}}} = \{\boldsymbol{n_1}, \boldsymbol{n_2}, \boldsymbol{n_3}\}\) is used to define the rotation \(R_{ij} = \boldsymbol{n_j} . \boldsymbol{e_i}\). And \(\boldsymbol{C}\) can be rotated in the global basis \(\Omega = \{\boldsymbol{e_1}, \boldsymbol{e_2}, \boldsymbol{e_3}\}\) as follow:
Linear orthotropic¶
A particular case of anisotropy is when the material basis is orthogonal in which case the elastic modulus tensor can be simplified and rewritten in terms of 9 independents material parameters.
The Poisson ratios follow the rule \(\nu_{ij} = \nu_{ji} E_i / E_j\).
NeoHookean¶
The hyperelastic NeoHookean constitutive law results from an extension of the linear elastic relationship (Hooke’s Law) for large deformation. Thus, the model predicts nonlinear stressstrain behavior for bodies undergoing large deformations.
As illustrated in Figure 4.6, the behavior is initially linear and the mechanical behavior is very close to the corresponding linear elastic material. This constitutive relationship, which accounts for compressibility, is a modified version of the one proposed by Ronald Rivlin [Belytschko:2000].
The strain energy stored in the material is given by:
where \(\lambda_0\) and \(\mu_0\) are, respectively, Lamé’s first parameter and the shear modulus at the initial configuration. \(J\) is the jacobian of the deformation gradient (\(\boldsymbol{F}=\nabla_{\!\!\boldsymbol{X}}\boldsymbol{x}\)): \(J=\text{det}(\boldsymbol{F})\). Finally \(\boldsymbol{C}\) is the right CauchyGreen deformation tensor.
Since this kind of material is used for large deformation problems, a finite deformation framework should be used. Therefore, the Cauchy stress (\(\boldsymbol{\sigma}\)) should be computed through the second PiolaKirchhoff stress tensor \(\boldsymbol{S}\):
Finally the second PiolaKirchhoff stress tensor is given by:
The parameters to indicate in the material file are the same as those
for the elastic case: E
(Young’s modulus), nu
(Poisson’s ratio).
ViscoElasticity¶
Viscoelasticity is characterized by strain rate dependent behavior. Moreover, when such a material undergoes a deformation it dissipates energy. This dissipation results in a hysteresis loop in the stressstrain curve at every loading cycle (see Figure [fig:smm:cl:viscoelastic:hyst]). In principle, it can be applied to many materials, since all materials exhibit a viscoelastic behavior if subjected to particular conditions (such as high temperatures).
The standard rheological linear solid model (see Sections 10.2 and 10.3
of [simo92]) has been implemented in Akantu
. This
model results from the combination of a spring mounted in parallel with
a spring and a dashpot connected in series, as illustrated in
Figure [fig:smm:cl:viscoelastic:model].
The advantage of this model is that it allows to account for creep or
stress relaxation. The equation that relates the stress to the strain is
(in 1D):
where \(\eta\) is the viscosity. The equilibrium condition is unique and is
attained in the limit, as \(t \to \infty\). At this stage, the response is
elastic and depends on the Young’s modulus \(E\). The mandatory parameters
for the material file are the following: rho
(density), E
(Young’s
modulus), nu
(Poisson’s ratio), Plane_Stress
(if set to zero plane
strain, otherwise plane stress), eta
(dashpot viscosity) and Ev
(stiffness of the viscous element).
Note that the current standard linear solid model is applied only on the deviatoric part of the strain tensor. The spheric part of the strain tensor affects the stress tensor like an linear elastic material.
SmallDeformation Plasticity¶
The smalldeformation plasticity is a simple plasticity material formulation which accounts for the additive decomposition of strain into elastic and plastic strain components. This formulation is applicable to infinitesimal deformation where the additive decomposition of the strain is a valid approximation. In this formulation, plastic strain is a shearing process where hydrostatic stress has no contribution to plasticity and consequently plasticity does not lead to volume change. Figure 4.7 shows the linear strain hardening elastoplastic behavior according to the additive decomposition of strain into the elastic and plastic parts in infinitesimal deformation as
In this class, the von Mises yield criterion is used. In the von Mises yield criterion, the yield is independent of the hydrostatic stress. Other yielding criteria such as Tresca and Gurson can be easily implemented in this class as well.
In the von Mises yield criterion, the hydrostatic stresses have no effect on the plasticity and consequently the yielding occurs when a critical elastic shear energy is achieved.
where \(\sigma_y\) is the yield strength of the material which can be function of plastic strain in case of hardening type of materials and \({\boldsymbol{\sigma}}^{{\mathrm{tr}}}\) is the deviatoric part of stress given by
After yielding \((f = 0)\), the normality hypothesis of plasticity determines the direction of plastic flow which is normal to the tangent to the yielding surface at the load point. Then, the tensorial form of the plastic constitutive equation using the von Mises yielding criterion (see equation 4.34) may be written as
In these expressions, the direction of the plastic strain increment (or equivalently, plastic strain rate) is given by \(\frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}}\) while the magnitude is defined by the plastic multiplier \(\Delta p\). This can be obtained using the consistency condition which impose the requirement for the load point to remain on the yielding surface in the plastic regime.
Here, we summarize the implementation procedures for the smalldeformation plasticity with linear isotropic hardening:
Compute the trial stress:
\[{\boldsymbol{\sigma}}^{{\mathrm{tr}}} = {\boldsymbol{\sigma}}_t + 2G\Delta \boldsymbol{\varepsilon} + \lambda \mathrm{tr}(\Delta \boldsymbol{\varepsilon})\boldsymbol{I}\]Check the Yielding criteria:
\[f = (\frac{3}{2} {\boldsymbol{\sigma}}^{{\mathrm{tr}}} : {\boldsymbol{\sigma}}^{{\mathrm{tr}}})^{1/2}\sigma_y (\boldsymbol{\varepsilon}^p)\]Compute the Plastic multiplier:
\[\begin{split}\begin{aligned} d \Delta p &= \frac{\sigma^{tr}_{eff}  3G \Delta P^{(k)} \sigma_y^{(k)}}{3G + h}\\ \Delta p^{(k+1)} &= \Delta p^{(k)}+ d\Delta p\\ \sigma_y^{(k+1)} &= (\sigma_y)_t+ h\Delta p \end{aligned}\end{split}\]Compute the plastic strain increment:
\[\Delta {\boldsymbol{\varepsilon}}^p = \frac{3}{2} \Delta p \frac{{\boldsymbol{\sigma}}^{{\mathrm{tr}}}}{\sigma_{{\mathrm{eff}}}}\]Compute the stress increment:
\[{\Delta \boldsymbol{\sigma}} = 2G(\Delta \boldsymbol{\varepsilon}\Delta \boldsymbol{\varepsilon}^p) + \lambda \mathrm{tr}(\Delta \boldsymbol{\varepsilon}\Delta \boldsymbol{\varepsilon}^p)\boldsymbol{I}\]Update the variables:
\[\begin{split}\begin{aligned} {\boldsymbol{\varepsilon^p}} &= {\boldsymbol{\varepsilon}}^p_t+{\Delta {\boldsymbol{\varepsilon}}^p}\\ {\boldsymbol{\sigma}} &= {\boldsymbol{\sigma}}_t+{\Delta \boldsymbol{\sigma}} \end{aligned}\end{split}\]
We use an implicit integration technique called the radial return method to
obtain the plastic multiplier. This method has the advantage of being
unconditionally stable, however, the accuracy remains dependent on the step
size. The plastic parameters to indicate in the material file are:
\(\sigma_y\) (Yield stress) and h
(Hardening modulus). In addition, the
elastic parameters need to be defined as previously mentioned: E
(Young’s
modulus), nu
(Poisson’s ratio).
Damage¶
In the simplified case of a linear elastic and brittle material, isotropic damage can be represented by a scalar variable \(d\), which varies from \(0\) to \(1\) for no damage to fully broken material respectively. The stressstrain relationship then becomes:
where \(\boldsymbol{\sigma}\),
\(\boldsymbol{\varepsilon}\) are the Cauchy stress and
strain tensors, and \(\boldsymbol{C}\) is the elastic
stiffness tensor. This formulation relies on the definition of an
evolution law for the damage variable. In Akantu
, many possibilities
exist and they are listed below.
Marigo¶
This damage evolution law is energy based as defined by Marigo [marigo81a][ lemaitre96a]. It is an isotropic damage law.
In this formulation, \(Y\) is the strain energy release rate, \(Y_d\) the rupture criterion and \(S\) the damage energy. The nonlocal version of this damage evolution law is constructed by averaging the energy \(Y\).
Mazars¶
This law introduced by Mazars [mazars84a] is a behavioral model to represent damage evolution in concrete. This model does not rely on the computation of the tangent stiffness, the damage is directly evaluated from the strain.
The governing variable in this damage law is the equivalent strain \(\varepsilon_{{\mathrm{eq}}} = \sqrt{<\boldsymbol{\varepsilon}>_+:<\boldsymbol{\varepsilon}>_+}\), with \(<.>_+\) the positive part of the tensor. This part is defined in the principal coordinates (I, II, III) as \(\varepsilon_{{\mathrm{eq}}} = \sqrt{<\boldsymbol{\varepsilon_I}>_+^2 + <\boldsymbol{\varepsilon_{II}}>_+^2 + <\boldsymbol{\varepsilon_{III}}>_+^2}\). The damage is defined as:
With \(\kappa_0\) the damage threshold, \(A_t\) and \(B_t\) the damage parameter in traction, \(A_c\) and \(B_c\) the damage parameter in compression, \(\beta\) is the shear parameter. \(\alpha_t\) is the coupling parameter between traction and compression, the \(\varepsilon_i\) are the eigenstrain and the \(\varepsilon_{{\mathrm{nd}}\;i}\) are the eigenvalues of the strain if the material were undamaged.
The coefficients \(A\) and \(B\) are the postpeak asymptotic value and the decay shape parameters.
Adding a New Constitutive Law¶
There are several constitutive laws in Akantu
as described in the previous
Section Constitutive Laws. It is also possible to use a userdefined 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 (material_XXX.hh
and
material_XXX.cc
) have to be provided where XXX
is the name of the new
material. The header file material_XXX.hh
defines the interface of your
custom material. Its implementation is provided in the material_XXX.cc
. The
new law must inherit from the 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:
auto & solver = model.getNonLinearSolver();
solver.set("max_iterations", 1);
solver.set("threshold", 1e4);
solver.set("convergence_type", SolveConvergenceCriteria::_residual);
model.solveStep();
/*  */
#include "material.hh"
/*  */
#ifndef __AKANTU_MATERIAL_XXX_HH__
#define __AKANTU_MATERIAL_XXX_HH__
namespace akantu {
class MaterialXXX : public Material {
/// declare here the interface of your material
};
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é
coefficients and its bulk modulus. Furthermore, the user has to add a new
internal variable 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:
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;
};
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:
for (UInt m = 0; m < model.getNbMaterials(); ++m)
std::cout << model.getMaterial(m) << std::endl;
the standard output stream operator has to be redefined. This should be done at the end of the header file:
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;
}
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 material_XXX.cc
file. However, the declaration has to be provided in the
material_XXX.hh
file:
class LocalMaterialDamage : public Material {
/*  */
/* Constructors/Destructors */
/*  */
public:
LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = "");
};
The user can now define the implementation of the constructor in the
material_XXX.cc
file:
/*  */
#include "local_material_damage.hh"
#include "solid_mechanics_model.hh"
namespace 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();
}
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:
registerParam
:
void registerParam(name of the parameter (key in the material file),
member variable,
default value (optional parameter),
access permissions,
description);
The available access permissions are as follows:
 _pat_internal
: Parameter can only be output when the material is printed.
 _pat_writable
: User can write into the parameter. The parameter is output when the material is printed.
 _pat_readable
: User can read the parameter. The parameter is output when the material is printed.
 _pat_modifiable
: Parameter is writable and readable.
 _pat_parsable
: Parameter can be parsed, i.e. read from the input file.
 _pat_parsmod
: Parameter is modifiable and parsable.
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:
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);
In the following a detailed description of these functions is provided:
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.computeStress
: In this method, the stresses are computed based on theconstitutive 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:
\[\mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon}\]Therefore, this method contains a loop on all quadrature points assigned to the material using the two macros:
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN
andMATERIAL_STRESS_QUADRATURE_POINT_LOOP_END
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type); // sigma < f(grad_u) MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
The strain vector in Akantu contains the values of \(\nabla \vec{u}\), i.e. it is really the displacement gradient,
computeTangentStiffness
: This method is called when the tangent to thestressstrain 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:
\[\label{eqn:smm:constitutive_elasc} \mat{K } =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}}\]Therefore, in this method, the
tangent
matrix (mat{D}) is computed for a given strain.The
tangent
matrix is a \(4^{th}\) order tensor which is stored as a matrix in Voigt notation.
getCelerity
: The stability criterion of the explicit integration scheme depend on the fastest wave celerity~eqref{eqn:smm:explicit:stabletime}. This celerity depend on the material, and therefore the value of this velocity should be defined in this method for each new material. By default, the fastest wave speed is the compressive wave whose celerity can be defined ingetPushWaveSpeed
.
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:
#include "material_XXX.hh"
For existing materials, as mentioned in Section~ref{sect:smm:CL}, by
default, the materials are initialized inside the method
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:
/// model initialization
model.initFull(_analysis_method = _explicit_lumped_mass);
Once the model has been initialized, the local material needs to be registered in the model:
model.registerNewCustomMaterials<XXX>("name_of_local_material");
Only at this point the material can be initialized:
model.initMaterials();
A full example for adding a new damage law can be found in
examples/new_material
.
Adding a New NonLocal Constitutive Law¶
In order to add a new nonlocal material we first have to add the local
constitutive law in Akantu (see above). We can then add the nonlocal version
of the constitutive law by adding the two files (material_XXX_non_local.hh
and material_XXX_non_local.cc
) where XXX
is the name of the
corresponding local material. The new law must inherit from the two classes,
nonlocal parent class, such as the MaterialNonLocal
class, and from the
local version of the constitutive law, i.e. MaterialXXX
. It is therefore
necessary to include the interface of those classes in the header file of your
custom material and indicate the inheritance in the declaration of the class:
/*  */
#include "material_non_local.hh" // the nonlocal parent
#include "material_XXX.hh"
/*  */
#ifndef __AKANTU_MATERIAL_XXX_HH__
#define __AKANTU_MATERIAL_XXX_HH__
namespace akantu {
class MaterialXXXNonLocal : public MaterialXXX,
public MaterialNonLocal {
/// declare here the interface of your material
};
As members of the class we only need to add the internal fields to store the nonlocal quantities, which are obtained from the averaging process:
/*  */
/* Class members */
/*  */
protected:
InternalField<Real> grad_u_nl;
The following four functions need to be implemented in the nonlocal material:
/// initialization of the material
void initMaterial();
/// loop over all element and invoke stress computation
virtual void computeNonLocalStresses(GhostType ghost_type);
/// compute stresses after local quantities have been averaged
virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type)
/// compute all local quantities
void computeStress(ElementType el_type, GhostType ghost_type);
In the intialization of the nonlocal material we need to register the local quantity for the averaging process. In our example the internal field grad_u_nl is the nonlocal counterpart of the gradient of the displacement field (grad_u_nl):
void MaterialXXXNonLocal::initMaterial() {
MaterialXXX::initMaterial();
MaterialNonLocal::initMaterial();
/// register the nonlocal variable in the manager
this>model>getNonLocalManager().registerNonLocalVariable(
this>grad_u.getName(),
this>grad_u_nl.getName(),
spatial_dimension * spatial_dimension);
}
The function to register the nonlocal variable takes as parameters the name of the local internal field, the name of the nonlocal counterpart and the number of components of the field we want to average. In the computeStress we now need to compute all the quantities we want to average. We can then write a loop for the stress computation in the function computeNonLocalStresses and then provide the constitutive law on each integration point in the function computeNonLocalStress.