Handling of constraints in NuToPde

\(\newcommand\vec[1]{\boldsymbol{#1}}\)

Theory

From the principle of virtual displacements, it follows for a body in equilibrium:

\[ \delta W_{int}^{(t+1)} = \delta W_{ext}^{(t+1)}, \]

where the variation of the internal energy \(\delta W_{int}^{(t+1)}\) and the variation of the external energy \(\delta W_{ext}^{(t+1)}\) is given by

\begin{align*} \delta W_{int}^{(t+1)} &= \int_V\vec{\sigma}^T(\vec{d}^{(t+1)})\delta\vec{\varepsilon}^{(t+1)}\;dV\\ &= \int_V\left(\vec{\sigma}(\vec{d}^{(t)})+\Delta \vec\sigma\left(\vec{d}^{(t)}+\Delta \vec{d}^{(t+1)}\right)\right)^T\delta\vec\varepsilon^{(t+1)}\;dV \\ \delta W_{ext}^{(t+1)} &= \left(\vec F_{ext}^{(t)}+\Delta \vec{F}_{ext}^{(t+1)}\right)^T \delta \vec{d}^{(t+1)}. \end{align*}

Linearization of the system of equations gives

\[ \int_V \left( \dfrac{\partial \sigma(\vec{d}^{(t)})}{\partial \varepsilon}\Delta\vec\varepsilon^{(t+1)} \right)^T \delta\vec\varepsilon^{(t+1)}\;dV = \left( \Delta \vec F_{ext}^{(t+1)} \right)^T \delta \vec{d}^{(t+1)} - \int_V \left( \vec\sigma^{(t)} \right)^T \delta\vec\varepsilon^{(t+1)}\;dV + \left( \vec{F}_{ext}^{(t)} \right)^T \delta\vec d^{(t+1)}, \]

where the last two summands on the right hand side vanish, if the previous step is in equilibrium. In the general formulation, the vector \(\vec{d}\) includes all the nodal displacements. However, if kinematic boundary conditions have to be applied, the variations of the displacement state have to be kinematically compatible. In general, the linear constraint equations can be written as

\[ \vec{T}_1^{(t+1)} \vec{d}_1^{(t+1)}+\vec{T}_2^{(t+1)} \vec{d}_2^{(t+1)} = \vec{b}^{(t+1)}, \]

where the nodal displacements are separated into a vector of independent, unknown DOFs \(\vec{d}_1\) and a vector of dependent DOFs \(\vec{d}_2\).

We restrict the way, constraint equations can be build.

• Each constraint equation must only contain exactly one dependent dof. If you need to define a constraint equation with more than one dependent dof, you have to perform an elimination manually.
• The coefficient in front of the dependent dof term has to be 1.

With the correct DOF numbering, this ensures that \(\vec{T}_2\) is an identity matrix and we can trivially eliminate the dependent DOFs \(\vec{d}_2\).

\[ \vec{d}^{(t+1)} = \begin{bmatrix} \vec{d}_1^{(t+1)}\\ \vec{d}_2^{(t+1)} \end{bmatrix} = \begin{bmatrix} \vec{I}\\ -\vec{T}_{1} \end{bmatrix} \vec{d}_1^{(t+1)} + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix}, \]

The matrix \(\vec{T}_1\) is the so called constraint matrix and will now be called \(\vec{C}_{con} \).

\[ \vec{d}^{(t+1)} = \begin{bmatrix} \vec{d}_1^{(t+1)}\\ \vec{d}_2^{(t+1)} \end{bmatrix} = \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \vec{d}_1^{(t+1)} + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} \]

The matrix \(\vec{C}_{con} \) relates the dependent DOFs with the independent DOFs. In the further derivation it is assumed that the matrix \(\vec{C}_{con}\) does not depend on the load step, since the increase of the loading is only obtained by modification of the right hand side \(\vec{b}\).

Consequently, the variation of the nodal displacement vector is given by

\[ \delta\vec{d}^{(t+1)} = \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)}. \]

Using the element shape functions, the relation between nodal displacements \(\vec{d}\) and strains \(\vec{\varepsilon}\) can be expressed as

\begin{align*} \Delta\vec{\varepsilon}^{(t+1)} &= \vec{\varepsilon}^{(t+1)}- \vec{\varepsilon}^{(t)}\\ &= \vec{B}\left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \left( \vec{d}_1^{(t+1)}-\vec{d}_1^{(t)} \right) + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} - \begin{bmatrix} \vec{0}\\ \vec{b}^{(t)} \end{bmatrix} \right)\\ \delta\vec{\varepsilon}^{(t+1)} &= \vec{B}\left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)} \right). \end{align*}

Setting

\[ \Delta\vec{d}_1^{(t+1)}=\left(\vec{d}_1^{(t+1)}-\vec{d}_1^{(t)}\right) \]

and substituting of these equations into the linearized system of the virtual work equations gives

\[ \left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \Delta\vec{d}_1^{(t+1)}+ \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} - \begin{bmatrix} \vec{0}\\ \vec{b}^{(t)} \end{bmatrix} \right)^{T} \vec{K}^T \left( \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)} \right) = \left( \Delta \vec F_{ext}^{(t+1)} + \vec{F}_{ext}^{(t)} - \vec F_{int}^{(t)} \right)^T \begin{bmatrix} \vec{I}\\ -\vec{C}_{con} \end{bmatrix} \delta\vec{d}_1^{(t+1)}. \]

The integrated stiffness matrix and the internal and external force vectors are split into submatrices similar to \(\vec d_1\) and \(\vec d_2\):

\begin{align*} \int_V \vec{B}^T\dfrac{\partial \sigma(\vec{d}^{(t)})}{\partial \varepsilon}\vec{B}\;dV &= \vec K = \begin{bmatrix} \vec K_{11} & \vec K_{12}\\ \vec K_{21} & \vec K_{22} \end{bmatrix}\\ \int_V\vec{B}^T\vec\sigma^{(t)}\;dV&=\vec F_{int}^{(t)} = \begin{bmatrix} \vec F_{1,int}^{(t)}\\ \vec F_{2,int}^{(t)} \end{bmatrix}. \end{align*}

Canceling out \(\delta\vec{d}_1^{(t+1)}\), transposing and splitting the external and internal forces into two parts yields:

\[ \begin{bmatrix} \vec{I}& -\vec{C}_{con}^{T} \end{bmatrix} \begin{bmatrix} \vec K_{11} & \vec K_{12}\\ \vec K_{21} & \vec K_{22} \end{bmatrix} \left( \begin{bmatrix} \vec{I} -\vec{C}_{con} \end{bmatrix} \Delta\vec{d}_1^{(t+1)} + \begin{bmatrix} \vec{0}\\ \vec{b}^{(t+1)} \end{bmatrix} - \begin{bmatrix} \vec{0}\\ \vec{b}^{(t)} \end{bmatrix} \right) = \begin{bmatrix} \vec{I}& -\vec{C}^T_{con} \end{bmatrix} \begin{bmatrix} \Delta \vec F_{1,ext}^{(t+1)}+\vec F_{1,ext}^{(t)}-\vec F_{1,int}^{(t)}\\ \Delta \vec F_{2,ext}^{(t+1)}+\vec F_{2,ext}^{(t)}-\vec F_{2,int}^{(t)} \end{bmatrix}. \]

Expanding this equation finally gives:

\begin{align*} \vec{K}^{mod}\Delta\vec{d}_1^{(t+1)} = \left( \vec{K}_{12}-\vec{C}^T_{con}\vec{K}_{22} \right) \left( \vec{b}^{(t)}-\vec{b}^{(t+1)}\right)+\Delta \vec F_{1,ext}^{(t+1)}+\vec F_{1,ext}^{(t)}\\ -\vec F_{1,int}^{(t)} -\vec{C}^T_{con}\left(\Delta \vec F_2^{(t+1)}+\vec F_{2,ext}^{(t)}-\vec F_{2,int}^{(t)} \right) \end{align*}

with the modified stiffness matrix

\[ \vec K^{mod} = \vec{K}_{11}-\vec{C}^T_{con}\vec{K}_{21}-\vec{K}_{12}\vec{C}_{con}+\vec{C}_{con}^T\vec{K}_{22}\vec{C}_{con}. \]

Implementation

Basics

• The class NuTo::ConstraintPde::Constraints has a container of NuTo::Constaint::Equation.
• Each NuTo::ConstraintPde::Equation has a time dependent right-hand-side (std::function<double(double time)>) and a container of NuTo::ConstraintPde::Term.
• Each term saves a reference to a NuTo::NodeSimple, a component and a coefficient.
• The first term of each equation is the dependent term. Its coefficient is always 1.
• The method NuTo::ConstraintPde::Constraints::BuildConstraintMatrix(...) returns a sparse matrix containing \(\vec{C}_{con}\)
• The method NuTo::ConstraintPde::Constraints::GetRhs() returns a the right hand side of the equations \(\vec{b}\)
• The constraints do not include a numbering of the dofs in NuTo::NodeSimple. This is done via NuTo::DofNumbering::Build() and has to be performed before you try to build \(\vec{C}_{con}\)

Usage

1) Create as many instances of NuTo::ConstraintPde::Equation as you need

• manually by instantiating an object of that class (not recommended)
• using mechanics/constraintsPde/ConstraintCompanion.h
  • Just have a look at this class, very intuitive, always returns an NuTo::ConstraintPde::Equation or a std::vector of those

2) Store the equation

• call NuTo::ConstraintPde::Constraints::Add(...)

3) Perform some kind of dof numbering

• manual
• using NuTo::DofNumbering::Build()