The contact formulation corresponds to the augmented-Lagrangian method \cite{Laursen:2002}, which seeks to minimize the following energy functional between two bodies with contact interface $\Gamma_c^{\left(1\right)} $:
\begin{equation}\label{eq:functinal}
\Pi^{al}\left( \boldsymbol u, \lambda_N \right) := \sum_{i=1}^2 \Pi^{\left( i \right) } ( \boldsymbol u ^{\left( i \right) } ) + \int_{\Gamma_c^{\left( 1 \right) }}\left[ \frac{1}{2 \epsilon_N } \left< \lambda_N + \epsilon_N g \right>^ 2 - \frac{1}{2 \epsilon_N} \lambda^2_N \right]\, d\Gamma,
\end{equation}
where $\Pi^{\left( i \right) } $ corresponds to the energy functional of the $ i $-th body in contact, which is a function of the displacement field $\boldsymbol u$ and the Lagrangian multiplier $\lambda_N$. In the equation above, $\epsilon_N$ is a penalty parameter, $g$ the contact gap function, and $\left< \cdot\right>$ the Macaulay bracket.
The solution procedure seeks to make the equation above stationary with respect to both $\boldsymbol u$ and $\lambda_N$:
\begin{equation}
\begin{aligned}
0 &= D_{\boldsymbol u}\Pi^{al}\cdot\boldsymbol w = G^{int,ext}\left( \boldsymbol u, \boldsymbol w \right) + \int_{\Gamma_c^{\left( 1 \right) }}\left< \lambda_N + \epsilon_N g \right> \delta g \, d\Gamma\qquad\forall\boldsymbol w \in\mathcal{V}\\
A solution is then found by using Uzawa's method, for which we solve for $\boldsymbol u ^ {\left( k \right) }$, with $\lambda_N^{\left( k \right)}$ fixed:
\begin{equation}
G^{int,ext} ( \boldsymbol u ^{\left( k \right) }, \boldsymbol w ) + \int_{\Gamma_c^{\left( 1 \right) }}\left< \lambda_N^{\left( k \right) } + \epsilon_N g( \boldsymbol u ^{\left( k \right) } ) \right> \delta g \, d\Gamma = 0 \qquad\forall\boldsymbol w \in\mathcal{V},
\end{equation}
followed by an update of the multipliers on $\Gamma_c ^ {\left(1\right) }$:
\begin{equation}
\lambda_N ^{\left( k+1 \right) } = \left< \lambda_N ^{\left( k \right) } + \epsilon_N g ( \boldsymbol u ^{\left( k \right) } ) \right>.
\end{equation}
It is worth noting that the results of the implicit contact resolution depend largely on the choice of the penalty parameter $\epsilon_N $. Depending on this parameter, the computational time needed to obtain a converged solution can be increased dramatically, or a convergence solution could not even be obtained at all.
The code provides a flag that allows the user to rely on an automatic value of $\epsilon_N $ for each slave node. Yet, this value should be used as a reference only, since for some problems it is actually overestimated and convergence cannot be obtained.
\subsection{Implementation}
In \akantu, the object that handles the implicit contact can be found in \code{implicit\_contact\_manager.hh}.
The object that handles the implicit contact resolution stage is the class template
\begin{cpp}
template <int dim, class model_type> struct ContactData;
\end{cpp}
This object takes the command line parameters during construction, which can be used to set up the behavior during contact resolution. The object can take the following parameters (default values in brackets):
\begin{tabular}{lrl}
-e & [auto] & Penalty parameter for augmented-Lagrangian formulation \\
-alpha & [1] & Multiplier for values of the penalty parameter\\
-utol & [0.001] & Tolerance used for multipliers in the Uzawa method\\
-ntol &[0.001]& Tolerance used in the Newton-Raphson inner convergence loop\\
-usteps &[100]& Maximum number of steps allowed in the Uzawa loop\\
-nsteps & [100]& Maximum number of steps allowed in the Newton-Raphson loop\\
\end{tabular}\\
Also, the following flags can be given to the command line
\begin{tabular}{ll}
-dump& Dumping within Newton iterations \\
-v& Verbose output
\end{tabular}\\
The \code{ContactData} object stores the state of the contact mechanics simulation. The state is contained within the following variables:
\begin{tabular}{ll}
\code{sm\_}& slave-master map \\
\code{multipliers\_}& Lagrange multiplier map \\
\code{areas\_}& slave areas map \\
\code{penalty\_}& penalty parameter map \\
\code{gaps\_}& gap function map \\
\code{model\_}& reference to solid mechanics model \\
\code{multiplier\_dumper\_}& structures used to dump multipliers \\
\code{pressure\_dumper\_}& structures used to dump pressure \\
\code{options\_}& options map \\
\code{flags\_}& flags map \\
\code{uiter\_}, \code{niter\_}& Uzawa and Newton iteration counters
\end{tabular}\\
The interface of the \code{ContactData} object contains three methods to solve for each contact step, which is overloaded depending on the parameters passed. Their signatures are as follows
The second method allows the user to provide a pointer to an object that is used to search slave-master pairs. This can be done, for example, when due to the deformed configuration current slave-master pairs are no longer valid.
The last method in the snippet above allows the user to provide a functor that is called after the assembly of the contact contributions to the stiffness matrix and the force vector. The last takes place within the method \code{computeTangentAndResidual()}.
\subsubsection{Hertz Example}
Here we outline, step by step, the use of the implicit contact solver to obtain the solution of Hertzian contact. The complete implementation can be found in \code{examples/contact/hertz\_3D.cc}.
The following class is used as the object that will perform the search for new contact elements when a slave node is found to lie outside its master element. The class derives from a \code{SearchBase} class, and implements the virtual method \code{search}.
\begin{cpp}
struct Assignator : public ContactData<3,SolidMechanicsModel>::SearchBase {
We initialize the library, create a mesh, and set the solid mechanics model up:
\begin{cpp}
initialize("steel.dat", argc, argv);
// create mesh
Mesh mesh(dim);
// read mesh
mesh.read("hertz_3D.msh");
// create model
model_type model(mesh);
SolidMechanicsModelOptions opt(_static);
// initialize material
model.initFull(opt);
\end{cpp}
Then we create the contact data, which can be used to solve the contact problem. Note that some of the parameters required by the contact object can be coded in the implementation file (these can also be passed as arguments).
\begin{cpp}
// create data structure that holds contact data
contact_type cd(argc, argv, model);
// optimal value of penalty multiplier
cd[Alpha] = 0.05;
cd[Uzawa_tol] = 1.e-2;
cd[Newton_tol] = 1.e-2;
\end{cpp}
Next we find the area that corresponds to each slave node. For this we use the fact that if we apply a unit distributed load over the contact surface, the resulting force vector at each slave node has magnitude that is equal to the area.
Next we add the slave-master pairs for the analysis. We use a bounding box to consider only a fraction of the slave nodes in the model. Those slave nodes that are not within the bounding box are not considered in the analysis:
\begin{cpp}
// set-up bounding box to include slave nodes that lie inside it
Real l1 = 1.;
Real l2 = 0.2;
Real l3 = 1.;
point_type c1(-l1 / 2, -l2 / 2, -l3 / 2);
point_type c2(l1 / 2, l2 / 2, l3 / 2);
bbox_type bb(c1, c2);
// search policy for slave-master pairs
Assignator a(model);
// loop over nodes in contact surface to create contact elements
cout << "- Adding slave-master pairs" << endl;
for (auto nit = cs.node_begin(); nit != cs.node_end(); ++nit) {
point_type p(&coords(*nit));
// ignore slave node if it doesn't lie within the bounding box
We then start the loop over displacement increments, and at each step we call {{{solveContactStep}}} and save the displacement, resulting force, and maximum pressure, in an array that will be used to print the results at the end of the simulation:
\begin{cpp}
const size_t steps = 30;
Real data[3][steps]; // store results for printing