COMSOL 6.4 - Weak Inequality Constraint

To add a node on a domain, boundary, edge, or point in any physics interface, click the button (

) and select in the dialog. Then, depending on the geometric entity level, select > at the domain or boundary level, > , or > from the context menu. There is no global weak inequality constraint option.

The node allows specifying a complementary inequality constraint of the form

where mi is the constraint expression and fc is the Lagrange multiplier field enforcing the constraint. The Lagrange multiplier usually has a physical interpretation as a reaction force or reaction flux. The third condition is a complementarity condition. The meaning of this condition is that the Lagrange multiplier can only be nonzero at points where the constraint is active (mi = 0), and must be zero where the constraint is inactive (mi < 0). For example, in a structural contact problem, the contact force (the Lagrange multiplier) must only be nonzero at points where the contact gap is zero.

This type of inequality constraint is strongly nonlinear and cannot be satisfied exactly until one knows at which points the constraint is active in the converged solution. This is a combinatorial problem that cannot be handled by the standard nonlinear solver. The feature therefore adds a relaxed approximation to the above conditions to the overall system of equations, accepting a slight violation of the inequality constraint as well as the complementarity condition mifc = 0.

The default approximation to the inequality constraint conditions uses a penalty method, also know as a stiff spring method. The constraint force, which is an approximation to the true Lagrange multiplier, is computed as an explicit function of the violation

with a user-defined spring constant ks. Increasing the spring constant reduces the violation of active constraints, but at the same time makes the problem stiffer and more difficult for the nonlinear solver.

The discontinuous derivative of the max() operator may also contribute to solver nonconvergence. Therefore, the feature can optionally replace the max() operator with a smoothed ramp() function with a user-defined transition zone size (dt). The transition zone extends on both sides of mi = 0, meaning that the Lagrange multiplier fc can be positive also at locations where mi is slightly negative. In a contact problem, this in practice means that the contact force will start to increase slightly before contact is made between the surfaces.

While the penalty method in many cases provides a fair approximation of the inequality constraint, its solution is inherently dependent on the unphysical spring constant ks. You always end up in a tradeoff between constraint violations and solvability. As a way to circumvent this, the feature optionally implements an augmented Lagrangian update strategy.

When the augmented Lagrangian method is used, the Lagrange multiplier is represented as a dependent variable field, λa, which is updated in a segregated iteration

The correction applied in each iteration, j, is effectively computed using the penalty method and a user-defined spring constant. This iteration can converge such that λa,j+1 = λa,j only when either mi,j < 0 and λa,j = 0, or mi,j = 0. This holds independently of the value of the spring constant ks, which therefore now only controls the convergence rate and not the accuracy of the final solution. A higher value of the spring constant may lead to faster convergence, but also to divergence. Conversely, a lower value generally leads to a slower but more stable convergence.

The Lagrange multiplier fc is applied on the overall system of equations as an additive weak term fc*Nf, where Nf is a constraint Jacobian expression containing the test() or var() operators. The constraint Jacobian in practice controls how the constraint force affects the equations. By default, Nf is set to test(fc), which means that active constraints are enforced in a symmetric way. The weak contribution is integrated using a specified integration order. By default the order is taken from the surrounding physics interface, based on its current discretization order.

When using the augmented Lagrangian method, the Lagrange multiplier field λa must be discretized using a chosen shape function and shape function order. In practice, Lagrange and Discontinuous Lagrange shape functions are the most useful. Constraints are better respected locally when using discontinuous elements or higher shape function orders, but convergence behavior may suffer. The augmented Lagrangian update equation is set up as a weak equation using an integration order matching the Lagrange multiplier shape order.

Solution of the augmented Lagrangian system of equations requires the use of a Segregated stationary solver, where the Lagrange multiplier variable is solved for in the last step in each iteration. When constant discontinuous or first-order Lagrange shape functions are used, this update step may use a Lumped Solver. Otherwise, it should be a standard Segregated Step with default settings. If the augmented Lagrangian method is used in some Weak Inequality Constraint feature, the default solver is automatically modified in this way.

The default solver settings for other variables are not modified. In theory, the augmented Lagrangian iteration expects the remaining variables to first be solved to convergence in each segregated iteration. This means that you may have to change the for these variables from the default to with a sufficient maximum .

Select an option from the list: (the default) or . For either option, enter a , which is to be set less than or equal to 0. For example, entering 2-(u+v) constrains u+v to be larger than or equal to 2.

For , also enter a . The constraint force expression must use the test() or var() operator. For example, write test(u) to enforce the constraint by modifying only the u equation with positive reaction term.

	The sign of the constraint force expression is important. Changing the sign of the constraint force expression does not only change the sign of the Lagrange multiplier as is the case for a standard Weak Constraint. Instead it changes the meaning of the complementary conditions and usually leads to a diverging solution.

From the list, select a method for approximating the inequality constraint: (the default) or . For both methods, specify a value for the to be used when enforcing the constraint. The appropriate value depends both on the scale of the main system of equations (in practice the material model) and on the scale of the constraint expression. It must normally be chosen from experience or using trial and error.

Enable to relax the complementarity condition in a zone around where the constraint becomes active. This may improve convergence. Specify a which may be interpreted as in the same units as the constraint expression. For example, if the is some function g(u) and the is set to 0.1, then the reaction term will start to apply already when g(u) becomes larger than -0.1.

When the selected is , specify an for the Lagrange multiplier. The zero default is usually appropriate, but a positive value can be used to prevent large constraint violations during the first iterations, in particular if the is chosen relatively weak.

These settings affect the numerical integration, and you do normally not need to change them. The checkbox is selected by default, meaning that the settings are taken from the main equation in the interface. If the checkbox is cleared, the following settings become available:

The is a positive integer that specifies the desired accuracy of integration during discretization. Polynomials of at most the given integration order are integrated without systematic errors. For smooth constraint forces, a sufficient integration order is typically twice the order of the shape functions used in the constraint. For example, the default integration order for second-order Lagrange elements is 4.

The constraint forces generated by inequality constraints are in general continuous, but because of the presence of the max() operator they are not always polynomials over each element. Raising the integration order may therefore in some situations improve accuracy.

The setting determines which frame to base the integration on: , , , or . The default frame is the one used for the physics interface.

	This section is available for 2D axisymmetric and 1D axisymmetric models.

By default, the π checkbox is selected to make the Lagrange multiplier represent, for example, the constraint force per area rather than by length and a full revolution. If the checkbox is cleared, the Lagrange multiplier is not multiplied by 2π r where it is used in the weak equation, and therefore represents flux per length and a full revolution.