Skip to content

Multi Point Constraints (MPC)

Multi Point Constraints (MPCs) enforce relations among the degrees of freedom (DOF) at two or more distinct nodes in a Finite Element model. In numgeo MPCs are imposed using the Penalty function method.

The penalty functions imposing the constraints can be physically interpreted as forces acting in fictitious high stiffness elastic elements connecting the targeted DOF and enforcing the constraints approximately. The equations \(\boldsymbol{z}(\boldsymbol{d})\) constraining the DOF in \(\boldsymbol{d}\) can be written as follows 1:

\[ \begin{equation} \boldsymbol{z} = \boldsymbol{A} \boldsymbol{d} - \boldsymbol{z}_0 = \boldsymbol{0} \end{equation} \]

where \(\boldsymbol{A}\) and \(\boldsymbol{z}_0\) contain constants. Note that there are more DOF than constraint equations (\(\boldsymbol{A}\) has more columns than rows). If the constraints are satisfied \(\boldsymbol{z}\) vanishes (\(\boldsymbol{z}=0\)). To enforce the constrained, the overall system of equations:

\[ \begin{equation}\label{eq:soe} \boldsymbol{K} \boldsymbol{d} = \boldsymbol{f} \end{equation} \]

is rewritten as total energy potential and augmented by a penalty function:

\[ \begin{align} \nonumber \Pi &= \frac{1}{2} \boldsymbol{d}^T \boldsymbol{K} \boldsymbol{d} - \boldsymbol{d}^T \boldsymbol{f} + \frac{1}{2} \boldsymbol{z}^T(\boldsymbol{d}) \alpha \boldsymbol{z}(\boldsymbol{d}) \\ &= \frac{1}{2} \boldsymbol{d}^T \boldsymbol{K} \boldsymbol{d} - \boldsymbol{d}^T \boldsymbol{f} + \frac{1}{2} \alpha \left(\boldsymbol{A} \boldsymbol{d} - \boldsymbol{z}_0\right)^T \left(\boldsymbol{A} \boldsymbol{d} - \boldsymbol{z}_0\right) \end{align} \]

with the penalty factor \(\alpha\). As \(\alpha\) increases, the penalty constraint violation becomes greater. The minimum of the potential \(\Pi\) (and thus satisfaction of the constraint equations) is found by minimizing the variation of \(\Pi\) with respect to \(\boldsymbol{d}\):

\[ \begin{align} \nonumber \frac{\partial \Pi}{\partial \boldsymbol{d}}& = \boldsymbol{K} \boldsymbol{d} - \boldsymbol{f} + \alpha \boldsymbol{A}^T\boldsymbol{d}\boldsymbol{A} - \alpha \boldsymbol{A}^T\boldsymbol{z}_0 = 0 \\ \frac{\partial \Pi}{\partial \boldsymbol{d}}& = \left(\boldsymbol{K} + \alpha \boldsymbol{A}^T\boldsymbol{A}\right)\boldsymbol{d} - \boldsymbol{f} - \alpha \boldsymbol{A}^T\boldsymbol{z}_0 = 0 \end{align} \]

Accounting for the constraint equations, Equation (\(\ref{eq:soe}\)) now reads:

\[ \begin{equation}\label{eq:soe_mpc} \left(\boldsymbol{K} + \boldsymbol{K}^p \right) \boldsymbol{d} = \boldsymbol{f} - \boldsymbol{f}^p \end{equation} \]

with the penalty matrix \(\boldsymbol{K}^p = \alpha \boldsymbol{A}^T\boldsymbol{A}\) and the corresponding contribution to the out of balance forces \(\boldsymbol{f}^p = \alpha \boldsymbol{A}^T\boldsymbol{z}_0\).

In comparison with Lagrange multipliers, penalty functions have the advantage of introducing no new variables to the system of equations. However, the following points should be kept in mind:

  • For penalty weights \(\alpha\) going to infinity, equation (\(\ref{eq:soe}\)) and (\(\ref{eq:soe_mpc}\)) coincide.
  • Particular care is needed in the choice of the penalty weight \(\alpha\) to balance the trade-off between reducing the constraint violation, which requires to increase the penalty weights, and limiting the solution error due to ill-conditioning with respect to inversion of the modified tangent stiffness matrix. This ill-conditioning increases for increasing values of the penalty weights.
  • The number of equations solved in remains the same as for unconstrained systems but the penalty matrix may significantly increase the bandwidth of the structural equations, depending on how DOFs are numbered and what DOFs are coupled by the constraint equations.

Penalty functions have the disadvantage that the penalty weight \(\alpha\) must be chosen in an allowable range: large enough to be effective but not so large as to provoke numerical difficulties. In numgeo the penalty weight is defined as follows:

\[ \begin{equation}\label{eq:penalty_alpha} \alpha = \max(K_{ii}) \alpha_0 \end{equation} \]

The parameter \(\alpha_0\) can be prescribed by the user, the default value is \(\alpha_0=10^8\). Note that in numgeo we rigorously enforce the the heuristic square root rule. The heuristic square root rule states, that the penalty factor \(\alpha\) should not exceed \(\sqrt{10^p}\), where \(p\) is the working machine precision (see 1). Thus, instead of equation (\(\ref{eq:penalty_alpha}\)) we use:

\[ \begin{equation} \alpha = \max \left\{ \max(K_{ii}) \alpha_0, 10^8 \right\} \end{equation} \]

Remarks

Implicit Dynamic Analysis: In implicit dynamic analysis numgeo enforces MPCs rigorously for the displacements. The velocities and accelerations are derived from the displacements with the relations defined by the dynamic integration operator (see Implicit dynamic analysis). Therefore, the velocities and accelerations satisfy the constraint only approximately. In the special case of linear MPCs (such as equalDof) and geometrically linear analysis the velocities obtained in this way satisfy this constraint exactly.

  1. Robert D. Cook and others. Concepts and applications of finite element analysis. John Wiley & Sons, 2007. ↩↩