OpenSees Cloud

OpenSees AMI

What Is a Good Penalty Number?

Original Post - 29 Nov 2023 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

I often see the penalty constraint handler used with seemingly high penalty numbers like the following:

ops.constraints('Penalty',1e18,1e18)

I’m not sure why these specific numbers are used so often in scripts, but I suspect these values were used in an old example and have been passed down. And while the 1e18 values might have worked for the original example, they may not work (lead to non-convergence) for your model because the penalty numbers to use depend on the stiffness as well as the base units of your model.

Any way, the first value, alphaSP, is the penalty number applied to single point constraints defined via the fix command–the fully fixed, the pins, the rollers. The second value, alphaMP, is the penalty number applied to multi-point constraints, e.g., rigid beams and rigid diaphragms.

Consider the simple beam divided into two elements of equal length, one flexible, the other rigid (using the rigidLink command with the -beam option). Further details on the model can be found in a previous post.

When the axial rigidity of element 2 is enforced, the horizontal displacements of nodes 2 and 3 should be equal. Likewise, the rotations of nodes 2 and 3 will be equal when the flexural rigidity of element 2 is enforced.

Multi-Point Constraints

The common case is to use the penalty constraint handler to enforce multi-point constraints. Performing analyses of the beam for successively larger alphaMP values with alphaSP held constant at 1e12 gives the following convergence for the relative displacements and rotations of nodes 2 and 3.

As the penalty number exceeds a few orders of magnitude times the reference axial and flexural stiffness of the beam, the multi-point constraints are enforced. Going beyond 1e19, the linear equation solver fails due to round off error, i.e, poor conditioning.

The script to perform this analysis is very simple.

for r in range(20):
   ops.wipe()

   #
   # Define model
   #

   ops.constraints('Penalty',1e12,10**r)

   #
   # Analyze model, get results
   #

Single-Point Constraints

While the constraint handlers in OpenSees are typically used for multi-point constraints, the handlers must also enforce single-point constraints. When using a constraint handler, the equations associated with single-point constraints are not simply removed from the system of equations like they are when using the Plain handler.

Flipping alpha values around, performing analyses of the beam for successively larger alphaSP values with alphaMP held constant at 1e12, i.e., enforcing the rigid beam but not the support conditions, gives the following convergence for the horizontal and rotational boundary condition at node 1.

Again, as the penalty number exceeds a few orders of magnitude times the reference axial and flexural stiffness of the beam, the boundary conditions are enforced.

Rules of Thumb

When figuring out good penalty numbers for your model, use values that are 4-5 orders of magnitude higher than some characteristic, or reference, stiffness of your model. What that reference stiffness is is somewhat squishy–you don’t need to be exact, just get in the ballpark. Also note that due to differing units of the DOFs in frame models, you should use the larger reference stiffness between translational and rotational DOFs.