7.11 Nonlinear Constraints

Nonlinear constraints are equations of the form \(f\left(x_{i}\right)=0\), where \(x_{i}\) represent nodal degrees of freedom. These constraints may be enforced using the penalty method, with or without augmented Lagrangian, or the method of Lagrange multipliers. Each method has its own advantages and disadvantages.

Penalty Method with Augmented Lagrangian

The penalty method employs the penalty parameter \(\varepsilon\) to add \(W_{p}=\frac{1}{2}\varepsilon f^{2}\) to the total work, where \(\varepsilon\) serves as the 'stiffness' or spring constant in the 'spring' whose 'elongation' is \(f\), such that \(W_{p}\) represents the stored energy in the spring. The contribution to the virtual work for this type of analysis is \(\delta W_{p}=\varepsilon f\,\delta f\), where

\[ \begin{equation} \delta f=\sum_{i}\frac{\partial f}{\partial x_{i}}\delta x_{i}\label{eq:nlc-virtual-force} \end{equation} \]

If we are concerned that the value of \(\varepsilon\) may become too large to produce good numerical convergence, we may augment this expression by using \(\delta W_{p}=\left(\lambda+\varepsilon f\right)\delta f\), where \(\lambda\) is the Lagrangian augmentation parameter, which gets updated at each iteration to help us avoid using a value of \(\varepsilon\) that's too large. We may let \(c=\lambda+\varepsilon f\) represent the constraint that needs to be satisfied. Substituing eq.\eqref{eq:nlc-virtual-force} into this relation allows us to express \(\delta W_{p}\) in matrix form,

\[ \begin{equation} \delta W_{p}=\left[\begin{array}{ccc} \delta x_{1} & \cdots & \delta x_{n}\end{array}\right]\left[\begin{array}{c} c\frac{\partial f}{\partial x_{1}}\\ \vdots\\ c\frac{\partial f}{\partial x_{n}} \end{array}\right]=\delta\boldsymbol{\mathbf{r}}\cdot\mathbf{\boldsymbol{f_{\mathit{p}}}}\,.\label{eq:nlc-pm-virtual-work} \end{equation} \]

The linearization of \(\delta W_{p}\) along increments \(\Delta x_{i}\) of the degrees of freedom employs

\[ \begin{equation} Df=\sum_{j}\frac{\partial f}{\partial x_{j}}\Delta x_{j}\label{eq:nlc-pm-f-linearization} \end{equation} \]

and takes the form

\[ \begin{equation} D\delta W_{p}=\left[\begin{array}{ccc} \delta x_{1} & \cdots & \delta x_{n}\end{array}\right]\left[\begin{array}{ccc} \varepsilon\left(\frac{\partial f}{\partial x_{1}}\right)^{2}+c\frac{\partial^{2}f}{\partial x_{1}^{2}} & \cdots & \varepsilon\frac{\partial f}{\partial x_{1}}\frac{\partial f}{\partial x_{n}}+c\frac{\partial^{2}f}{\partial x_{1}\partial x_{n}}\\ \vdots & \ddots & \vdots\\ \varepsilon\frac{\partial f}{\partial x_{n}}\frac{\partial f}{\partial x_{1}}+c\frac{\partial^{2}f}{\partial x_{1}\partial x_{n}} & \cdots & \varepsilon\left(\frac{\partial f}{\partial x_{n}}\right)^{2}+c\frac{\partial^{2}f}{\partial x_{n}^{2}} \end{array}\right]\left[\begin{array}{c} \Delta x_{1}\\ \vdots\\ \Delta x_{n} \end{array}\right]=\delta\boldsymbol{\mathbf{r}}\cdot\mathbf{\boldsymbol{K_{\mathit{p}}}\boldsymbol{\Delta r}}\,,\label{eq:nlc-pm-stiffness} \end{equation} \]

The stiffness matrix for the penalty constraint is symmetric.

Lagrangian augmentation is performed after the (quasi-)Newton method has converged. At that time, the augmented Lagrangians are updated as follows

\[ \begin{equation} \lambda^{\left(k+1\right)}\leftarrow\lambda^{\left(k\right)}+\varepsilon f\label{eq:nlc-pm-augment} \end{equation} \]

The penalty method with augmented Lagrangian does not increase the number of equations that need to be solved globally. However, it does require users to specify the parameter \(\varepsilon\), and to provide a convergence tolerance for Lagrange augmentation.

Example 1. For linear constraints of the form

\[ \begin{equation} f=a_{0}+\sum_{j}a_{j}x_{j}=0\label{eq:nlc-linear-constraint} \end{equation} \]

we find that \(\frac{\partial f}{\partial x_{i}}=a_{i}\) and the second derivatives are all zero, so that the constributions to the global system of equations take the form

\[ \begin{equation} \mathbf{f}_{\mathit{p}}=\left[\begin{array}{c} ca_{1}\\ \vdots\\ ca_{n} \end{array}\right],\;\mathbf{K}_{\mathit{p}}=\left[\begin{array}{ccc} \varepsilon a_{1}^{2} & \cdots & \varepsilon a_{1}a_{n}\\ \vdots & \ddots & \vdots\\ \varepsilon a_{n}a_{1} & \cdots & \varepsilon a_{n}^{2} \end{array}\right].\label{eq:nlc-pm-linear} \end{equation} \]

Method of Lagrange Multipliers

When using the method of Lagrange multipliers we have to include equations that solve for the multipliers. We let the work done by the constraint be represented by \(W_{c}=\lambda f\), where \(\lambda\) is the multiplier (if \(f\) is a vector equation, we can use \(W_{c}=\boldsymbol{\lambda}\cdot\mathbf{f}\)). The total work is then \(\widetilde{W}=W+W_{c}\) where \(W\) is the work of the unconstrained problem (i.e. the contributions from the internal and external forces.) Then, we evaulate the first varition of \(\delta\widetilde{W}=\delta W+\delta W_{c}=0\), where \(\delta W_{c}=\delta f+\delta\lambda\,f\). Since \(\delta f\) is given as per eq.\eqref{eq:nlc-virtual-force}, \(\delta W_{c}\) may be written in matrix form as

\[ \begin{equation} \delta W_{c}=\left[\begin{array}{cccc} \delta x_{1} & \cdots & \delta x_{n} & \delta\lambda\end{array}\right]\left[\begin{array}{c} \lambda\frac{\partial f}{\partial x_{1}}\\ \vdots\\ \lambda\frac{\partial f}{\partial x_{n}}\\ f \end{array}\right]=\delta\boldsymbol{r\cdot}\mathbf{f}_{c}\,.\label{eq:nlc-lm-virtual-work} \end{equation} \]

The linearization of \(\delta W_{c}\) along increments \(\Delta x_{j}\) and \(\Delta\lambda\) takes the form

\[ \begin{equation} D\delta W_{c}=\sum_{i}\delta x_{i}\frac{\partial f}{\partial x_{i}}\Delta\lambda+\sum_{i}\delta x_{i}\sum_{j}\lambda\frac{\partial^{2}f}{\partial x_{i}\partial x_{j}}\Delta x_{j}+\delta\lambda\,\sum_{j}\frac{\partial f}{\partial x_{j}}\Delta x_{j}\,.\label{eq:nlc-lm-linearization} \end{equation} \]

The corresponding stiffness matrix is

\[ \begin{equation} D\delta W_{c}=\left[\begin{array}{cccc} \delta x_{1} & \cdots & \delta x_{n} & \delta\lambda\end{array}\right]\left[\begin{array}{cccc} \lambda\frac{\partial^{2}f}{\partial x_{1}^{2}} & \cdots & \lambda\frac{\partial^{2}f}{\partial x_{1}\partial x_{n}} & \frac{\partial f}{\partial x_{1}}\\ \vdots & \ddots & \vdots & \vdots\\ \lambda\frac{\partial^{2}f}{\partial x_{n}\partial x_{1}} & \cdots & \lambda\frac{\partial^{2}f}{\partial x_{n}^{2}} & \frac{\partial f}{\partial x_{n}}\\ \frac{\partial f}{\partial x_{1}} & \cdots & \frac{\partial f}{\partial x_{n}} & 0 \end{array}\right]\left[\begin{array}{c} \Delta x_{1}\\ \vdots\\ \Delta x_{n}\\ \Delta\lambda \end{array}\right]=\delta\boldsymbol{r\cdot}\mathbf{K}_{c}\boldsymbol{\Delta r}\,,\label{eq:nlc-lm-stiffness} \end{equation} \]

This resulting system of equations requires \(\lambda\) to part of the solution, and to be updated at each iteration \(k\) according to \(\lambda^{\left(k+1\right)}\leftarrow\lambda^{\left(k\right)}+\Delta\lambda\). The stiffness matrix is symmetric for this analysis.

The method of Lagrange multipliers increases the number of equations that need to be solved globally. However, it does not require users to specify additional parameters. Note that the contribution to the stiffness matrix introduces a zero coefficient on the diagonal of the stiffness matrix for the equation corresponding to \(\Delta\lambda\). This means that these equations can only be solved with linear solvers that support pivoting. (Most linear solvers in FEBio, including the default pardiso solver, can handle zero-diagonal elements.)

Example 2. For linear constraints of the form

\[ f=a_{0}+\sum_{j}a_{j}x_{j}=0 \]

we find that \(\frac{\partial f}{\partial x_{i}}=a_{i}\) and second derivatives are all zero, so that the contributions to the global residual and stiffness matrix take the form

\[ \begin{equation} \mathbf{f}_{c}=\left[\begin{array}{c} \lambda a_{1}\\ \vdots\\ \lambda a_{n}\\ f \end{array}\right],\:\boldsymbol{\mathbf{K}_{c}}=\left[\begin{array}{cccc} 0 & \cdots & 0 & a_{1}\\ \vdots & \ddots & \vdots & \vdots\\ 0 & \cdots & 0 & a_{n}\\ a_{1} & \cdots & a_{n} & 0 \end{array}\right]\label{eq:nlc-lm-linear} \end{equation} \]