## Abstract

Penalty functions can be used to add constraints to systems of equations. In computational mechanics, stiffness-type penalties are the common choice. However, in dynamic applications stiffness penalties have the disadvantage that they tend to *decrease* the critical time step in conditionally stable time integration schemes, leading to increased CPU times for simulations. In contrast, inertia penalties *increase* the critical time step. In this paper, we suggest the simultaneous use of stiffness penalties and inertia penalties, which is denoted as the *bipenalty method*. We demonstrate that the accuracy of the bipenalty method is at least as good as (and usually better than) using either stiffness penalties or inertia penalties. Furthermore, for a number of finite elements (bar elements, beam elements and square plane stress/plane strain elements) we have derived ratios of the two penalty parameters such that their combined effect on the critical time step is neutral. The bipenalty method is thus superior to using stiffness penalties, because the decrease in critical time step can be avoided. The bipenalty method is also superior to using inertia penalties, since the constraints are realized with higher accuracy.

## 1. Introduction

The penalty method is a popular technique to impose constraints, such as Dirichlet conditions, tyings or contact conditions, to a system of equations. In its common format, a penalty function is an addition to the potential energy of the system that consists of the constraint that is to be enforced together with a user-defined number known as the penalty parameter. The interpretation of such penalty parameters is that of (very) stiff springs, and we will refer to these penalties as *stiffness penalties*. If the penalty parameter is infinitely large, the constraint is satisfied exactly; however, in practical applications finite magnitudes of the penalty parameter must be used, implying that the constraint is only enforced approximatively. Thus, there is an inherent dilemma in using penalty functions: for *modelling accuracy* the penalty parameter must be taken as large as possible, but for *numerical accuracy* (related to the machine precision of computers) the penalty parameter must not be too large, see for instance the discussion in Bathe (1996).

There is another drawback of using stiffness penalties, which manifests itself in time domain dynamic applications. As is well known, the speed of sound of a material is proportional to the square root of stiffness divided by mass. The use of stiffness penalties thus leads to an increase of the speed of sound. When a time integration algorithm is used that is only conditionally stable (such as explicit time integration), a critical time step exists such that numerical instabilities may occur if the time step for numerical time integration is chosen larger than this critical time step. The critical time step is inversely proportional to the speed of sound, thus there is an important implication of using stiffness penalties: stiffness penalties lead to a decrease of the critical time step, which implies that numerical simulations require longer CPU times (Belytschko & Neal 1991; Carpenter *et al.* 1991).

In recent years research has been carried out aimed at using lower values of the penalty parameters, while maintaining accuracy of constraint imposition, or aimed at penalty formulations that have a less adverse effect on the critical time step. Ilanko and coworkers (Ilanko & Dickinson 1999; Ilanko 2002*a*,*b*, 2003, 2005*a*; Askes & Ilanko 2006; Askes *et al.* 2008) proposed the use of positive as well as negative values for the penalty parameters and they showed that the exact solution is bracketed by the results obtained with either positive or negative penalties. This result can be used to interpolate between results obtained with moderately large (but not very large) positive and negative penalty parameters. More recently, it was suggested to use inertia penalties instead of stiffness penalties (Ilanko 2005*b*; Williams & Ilanko 2005; Hetherington & Askes 2009). Whereas stiffness penalties can be interpreted as stiff springs, inertia penalties can be regarded as heavy masses; both types of penalties limit the motion of the relevant degree of freedom (d.f.). There is an additional advantage of using inertia penalties instead of stiffness penalties: since the speed of sound is lowered by increasing the mass, the critical time step of a finite element penalized with inertia is *increased* (Hetherington & Askes 2009).

The opposite effects of stiffness penalties and inertia penalties on the critical time step naturally leads to the idea that they can be combined: the motion of a d.f. that is to be constrained can be penalized *simultaneously* with stiffness penalties and with inertia penalties. This approach will be denoted here as the *bipenalty method* and its development is the main aim of this paper. We will focus on two major aspects of the proposed bipenalty method:

—

*Effects on accuracy*. Engineering judgement suggests that the simultaneous use of stiffness penalties and inertia penalties provides similar or better accuracy of constraint imposition than using either of the two penalties with equal magnitude on its own. We will formulate a mathematical proof that underpins this notion.—

*Effects on stability*. Since stiffness penalties and inertia penalties have opposite effects on the critical time step, there are certain combinations of stiffness penalties and inertia penalties for which there is a*zero*net effect on the critical time step. For a number of finite element types it is possible to derive closed-form expressions for the ratio of stiffness penalty parameter over inertia penalty parameter for which the net effect on the critical time step is neutral—such a ratio will be denoted as the ‘critical penalty ratio’.

After reviewing penalty functions based on stiffness and inertia in §2, accuracy and stability of the bipenalty method will be treated in §§3 and 4, respectively. Numerical examples are presented in §5 to illustrate and corroborate the analytical results of the earlier sections. Throughout the paper, the emphasis is on *time domain* analysis; possible extension of the findings to *frequency domain* analysis are beyond the scope of this study.

## 2. Stiffness penalties and inertia penalties

The potential energy of a structure is usually written as
2.1
where **K** the stiffness matrix, **u** is a vector of displacements and **f** the vector of external forces. A constraint is applied to the *n*th d.f. and is its prescribed value. The penalized potential energy is accordingly written as
2.2
where *α*_{s} is a penalty parameter of the stiffness type.

The usual expression for the kinetic energy reads
2.3
where **M** is the mass matrix. The constraint introduced above is written in rate form and added to the kinetic energy as
2.4
where *α*_{m} is an inertial penalty parameter that penalizes the difference in actual velocity and prescribed velocity . The equations of motion follow from
2.5
The components of the mass penalty matrix **M**^{P}, stiffness penalty matrix **K**^{P} and penalty force vector are zero except for the entries related to the *n*th d.f.:
2.6
2.7
2.8
To keep the subsequent derivations simple, in the remainder of the paper we will assume that the prescribed displacement is constant in time, by which . However, varying values of would only change the penalty force vector, which turns out to be of minor importance in our derivations. It is also useful to introduce dimensionless penalty factors *p*_{m} and *p*_{s} that relate the penalty parameters *α*_{m} and *α*_{s} to the corresponding diagonal entries of **M** and **K** via *α*_{m}=*p*_{m}*M*_{nn} and *α*_{s}=*p*_{s}*K*_{nn}, respectively.

## 3. Accuracy of the bipenalty method

For the time integration of the equations of motion (2.5) we will use the Newmark method. Denoting the time-discretized counterparts of the displacements **u**, velocities and accelerations with **d**, **v** and **a**, respectively, the Newmark equations can be written as
3.1
3.2
where *γ* and *β* are user-defined parameters. Expressed in terms of the time-discretized accelerations **a** at time *t*+Δ*t*, the equations of motion read
3.3
Furthermore, we will assume initial values **a**^{0}=**v**^{0}=**0** as well as **d**^{0}=**0** except .

### (a) Convergence of the bipenalty method

Adopting the initialization indicated above, the penalized (*n*th) equation at time *t*=Δ*t* can be written as
3.4
We introduce *α*_{m}=*η*_{m}*p* and *α*_{s}=*η*_{s}*p*, where *η*_{m} and *η*_{s} are user-defined numbers to set the magnitude of inertia penalties relative to the magnitude of stiffness penalties, while *p* is a large number. Thus, convergence is to be verified for the case . Dividing equation (3.4) by *p* and rearranging yields
3.5
by which
3.6
where we have used that all are finite, since they are obtained from equation (3.3) that contains a positive definite system matrix. Thus, convergence to the same (exact) solution is observed irrespective of whether one uses stiffness penalties, inertia penalties or both.

In order to assess the influence of the error in the penalized d.f. on the other d.f., we compare the penalized system of equations given in equation (2.5) with the exact solution obtained via Lagrange multipliers. The latter is obtained from
3.7
where *λ* is the Lagrange multiplier of the constrained *n*th d.f. (that is, *λ* is the reaction force conjugated to *u*_{n}). Furthermore, **A** and *g* are the usual Lagrange multiplier contributions to left-hand-side and right-hand-side: all *A*_{i}=0 except *A*_{n}=1, and . Substracting the first set of equations of expression (3.7) from equation (2.5) yields
3.8
where the error **e** is defined as . For the case of a single constraint on the *n*th d.f., the *n*th equation of expression (3.8) can be simplified as
3.9
Applying the Newmark time integration to the error **e**, equation (3.9) can be written for the first time evaluation *t*=Δ*t* as
3.10
This last expression shows that the error is obtained from
3.11
where the vector of residuals **r** is defined as *r*_{i}=0 for all *i* except *i*=*n*, for which *r*_{n} is defined in equation (3.10). Therefore, the errors of all d.f. *e*_{i} are governed by the same scalar, namely the residual *r*_{n} of the constrained d.f. We can now make the following observations:

— In equation (3.6) it was shown that for , therefore for . In the context of equation (3.10), this can only hold if for .

— Since for , it follows from equation (3.10) that for , which is valid for all

*i*.

Thus, we have shown that the bipenalty method converges for all d.f. Results for subsequent time instants can be obtained with induction.

### (b) Rate of convergence of the bipenalty method

To compare the accuracy of the bipenalty method with that of using either stiffness penalties or inertia penalties, we investigate the rate of convergence of the penalized d.f., that is , for the case . Since it was established in equation (3.6) that all penalty methods converge to the same solution, it can be said that a smaller value of corresponds to a more accurate method.

From equation (3.4) it follows that
3.12
by which
3.13
The magnitude of is governed by the denominator in equation (3.13). Taking the same magnitude for each penalty parameter, it can be seen that the bipenalty method is at least as accurate as (or, in case *β*>0, higher than) using either stiffness penalties or inertia penalties. Again, induction may be used to extend this result to subsequent time steps.

## 4. Stability of the bipenalty method

The critical time step Δ*t*_{crit} that is relevant for conditionally stable time integrators is inversely proportional to the so-called ‘maximum sampling frequency’ *ω*_{max} that can be monitored by the finite element discretization (Hughes 2000), namely
4.1
where the critical sampling ratio *Ω*_{crit} is given for the Newmark time integrator and for a system without physical damping. This maximum sampling frequency is found by solving a dynamic eigenvalue problem for one finite element written as
4.2
and selecting the largest root of *ω*. We aim to establish a relationship between the two penalties *α*_{m} and *α*_{s} such that the maximum sampling frequency is not affected. A penalty ratio *r* is defined in terms of the dimensionless penalty factors as *r*=*p*_{s}/*p*_{m}. For large values of *r*, the stiffness penalty dominates which could lead to an increase in *ω*_{max} and, thus, an undesired decrease in the critical time step. Conversely, for small values of *r* the inertia penalty dominates, which lowers *ω*_{max} and increases the critical time step. We will seek the upper limit of *r* such that the maximum sampling frequency of a bipenalized element is identical to the maximum sampling frequency of an unpenalized element—this upper limit of *r* is called the ‘critical penalty ratio’ and denoted as *r*_{crit}.

The procedure we will follow is summarized as follows:

— Select an element type and specify the mass matrix, which is taken here as either the lumped mass matrix or the consistent mass matrix (these are related to discrete masses and distributed masses, respectively). In this paper, we will present results for a number of finite elements with increasing complexity, namely two-noded bar elements, two-noded beam elements and four-noded square plane stress/plane strain elements.

— Solve equation (4.2) with

*p*_{m}=*p*_{s}=0. This gives the maximum sampling frequency of an unpenalized element.— Substitute the unpenalized value of

*ω*_{max}into equation (4.2) and solve again, but now with*p*_{m}≠0 and*p*_{s}≠0. This gives a relation between*p*_{m}and*p*_{s}, from which*r*_{crit}can be determined.

Readers who wish to skip the technical details may want to refer directly to table 1, where we have summarized all obtained values for the critical penalty ratio.

### (a) Bar elements

The lumped mass matrix, consistent mass matrix and stiffness matrix of a two-noded bar element read
4.3
where *ρ*, *E*, *A* and *h* are the mass density, Young’s modulus, cross-sectional area and element size, respectively. Penalties are imposed by multiplying one of the diagonal terms with 1+*p*_{m} (for the mass matrices) or with 1+*p*_{s} (for the stiffness matrix).

When the lumped mass matrix is taken, the penalized eigenvalue problem reads
4.4
where is the elastic bar velocity (i.e. the speed of sound). The unpenalized eigenvalue problem is retrieved by taking *p*_{m}=*p*_{s}=0, which gives *ω*_{max}=2*c*_{e}/*h*. This value for *ω*_{max} is next substituted into the penalized eigenvalue problem. This renders 2*p*_{m}−*p*_{s}=0 and consequently *r*_{crit}=2.

With a consistent mass matrix, the penalized eigenvalue problem can be written as
4.5
The unpenalized maximum sampling frequency equals . When this value is used for *ω* in equation (4.5) with *p*_{m}≠0 and *p*_{s}≠0, we obtain 4*p*_{m}−*p*_{s}=0 so that *r*_{crit}=4.

### (b) Euler–Bernoulli beam elements

The consistent mass matrix and stiffness matrix of a two-noded Euler–Bernoulli beam element read
4.6
where *I* is the second moment of area of the cross section. Although there is no unique definition of the lumped mass matrix, most authors (Ishihara 1981; Kim 1993; Bathe 1996) use the following expression:
4.7
In deriving the penalized eigenvalue problems, we will distinguish between three cases, namely (i) penalized deflection, (ii) penalized rotation, and (iii) penalized deflection as well as penalized rotation. In terms of modifying the element matrices, this amounts to multiplying the first diagonal component, the second diagonal component, or both, with factors 1+*p*_{m} and 1+*p*_{s}. However, derivations will only be presented for the (most general) third case.

Using a lumped mass matrix and with penalties applied on deflection as well as rotation, the polynomial equation of the eigenvalue problem contains not only linear terms in *p*_{m} and *p*_{s} but also quadratic terms in *p*_{m} and *p*_{s}:
4.8
For an unpenalized element the maximum sampling frequency . With this value substituted back into equation (4.8) we obtain the following equation for the relation between *p*_{m} and *p*_{s}:
4.9
When *p*_{s} is resolved from this last expression, the critical penalty ratio is found as a function of the inertia penalty factor *p*_{m} as
4.10
For we obtain *r*_{crit}↑8 or *r*_{crit}↓2. These two asymptotic values coincide with the critical penalty ratios found for penalized deflection only and penalized rotation only, respectively (the derivations of which are not shown here). In case both deflection and rotation are penalized, the lower value for *r*_{crit} as obtained from equation (4.10) must be taken; since the relevant convergence is from above it is safe to use *r*=2.

The eigenvalue problem with a consistent mass matrix and penalties on deflection as well as rotation of the first node leads to an equation for *ω* as
4.11
With *p*_{m}=*p*_{s}=0 the maximum sampling frequency of an unpenalized element is found as . Substituting this value into equation (4.11) leads to
4.12
Resolving *p*_{s} and dividing the result by *p*_{m} yields
4.13
For it can be seen that *r*_{crit}↑260 or *r*_{crit}↓20. Similar to the beam element with lumped mass, these two asymptotic values coincide with the critical penalty ratios for penalized deflection only and penalized rotation only, respectively (again, these derivations are not presented here for the sake of brevity). For the case where both deflection and rotation are penalized, the value *r*=20 can be used safely.

### (c) Square plane stress/plane strain elements

Finally, we will derive critical penalty ratios for one particular continuum element, namely a four-noded quadrilateral element, with one penalized d.f. If this element is taken to be isoparametric with arbitrary shape, the Jacobian will appear in the derivations. To avoid this, we will assume that the element is square, which is in any case the recommended shape of the element for optimal behaviour. Our results for the critical time steps of this element are in accordance with those presented earlier in Ling & Cherukuri (2002).

Assuming unit thickness in the out-of-plane direction, the lumped mass, consistent mass and stiffness matrices for a four-node square plane stress element are given by
4.14
4.15
and
4.16
where **I**_{8×8} is an 8×8 identity matrix. The plane strain case can be retrieved from the plane stress case by using a fictitious Poisson ratio *ν*′ instead of *ν* and a fictitious Young’s modulus *E*′ instead of *E*, which are given by
4.17
In order to obtain more compact expressions, we define a dimensionless frequency *ζ* as *ζ*=*ω**h*/*c*_{e}.

The penalized eigenvalue problem for a square plane stress element with lumped mass is found to be
4.18
Taking *p*_{m}=*p*_{s}=0 yields , irrespective of the magnitude of *ν*. We substitute this value for *ζ* into equation (4.18) to obtain 6(*ν*+1)*p*_{m}−(3−*ν*)*p*_{s}=0, hence *r*_{crit}=6(1+*ν*)/(3−*ν*) for plane stress. The equivalent plane strain case is found by using the expression for *ν*′ as given in equation (4.17). As a result, *r*_{crit}=6/(3−4*ν*) for plane strain.

Using the consistent mass matrix, the eigenvalue problem for a square plane stress element leads to a polynomial equation in *ζ* as
4.19
The maximum dimensionless sampling frequency depends on the value of the Poisson ratio. For we obtain , which implies for the penalized system that 8(1+*ν*)*p*_{m}−(3−*ν*)*p*_{s}=0, therefore *r*_{crit}=8(1+*ν*)/(3−*ν*) for plane stress. The counterpart transition value for *ν* in plane strain is . For the maximum dimensionless sampling frequency of the unpenalized system is . For the penalized eigenvalue problem this value leads to 8*p*_{m}−(3−4*ν*)*p*_{s}=0 and therefore *r*_{crit}=8/(3−4*ν*) for plane strain. No solutions for *r*_{crit} could be found for the case in plane stress and in plane strain. In the numerical examples we will use the expressions for *r*_{crit} given above irrespective of *ν*, and verify these are safe for and , respectively.

## 5. Results

In this section we will present a selection of results on stability and accuracy of the bipenalty method. To report on stability and accuracy for all possible combinations of element types, mass matrices, time integration schemes and penalty ratios would be excessive, therefore the discussion is limited to the following issues:

—

*For bar elements*. Illustrate the destabilizing effects of choosing the penalty ratio*r*larger than its critical value. This is done for two combinations of mass matrix and time integration scheme.—

*For bar elements*. Illustrate how the penalty ratio*r*affects the accuracy of constraint imposition. This is again done for two combinations of mass matrix and time integration scheme.—

*For beam elements*. Verify the effects on stability of penalizing deflection, rotation or both. This is done for explicit central difference time integration with a lumped mass matrix only.—

*For square two-dimensional elements*. Demonstrate that the analytically obtained expressions for the critical penalty ratio*r*_{crit}are valid for values of the Poisson ratio . This is done for the central differences time integration scheme combined with a consistent mass matrix.

### (a) Bar problem—stability analysis

The first problem is a one-dimensional bar of length *L*=100 m subjected to a tensile traction *σ*=1 N m^{−2} on the left end. The traction is applied from time *t*=0 s onwards and kept constant for the duration of the simulation up to *t*=150 s. The right end of the bar is fixed, and this support condition is modelled with penalty functions. The material parameters are taken as *E*=1 N m^{−2} and *ρ*=1 kg m^{−3}. A finite element mesh consisting of 100 two-noded bar elements is used. We will show results for two combinations of Newmark parameters and mass matrices, namely (i) the explicit central differences scheme obtained via *β*=0 and with a lumped mass matrix, so that the critical time step Δ*t*_{crit}=1 s, and (ii) the Fox–Goodwin scheme found with and in combination with a consistent mass matrix, so that s. In either case, we apply a time step Δ*t*=Δ*t*_{crit} so that no instabilities are expected while the wave propagates through the bar—any instabilities that arise are thus entirely due to the penalization of the right end of the bar. Our purpose is not to cross compare these two sets of results, but instead to show that our findings on the critical penalty ratio hold for each combination of mass matrix and time integration scheme.

Figures 1 (using an inertia penalty factor *p*_{m}=10^{2}) and 2 (using an inertia penalty factor *p*_{m}=10^{4}) show the strain profiles across the bar at time *t*=150 s for the various combinations of Newmark parameters and mass matrices, while also stable and unstable values for the ratio *r*=*p*_{s}/*p*_{m} were taken. At time *t*=150 s the wave front has propagated through the entire bar, been reflected on the penalized right end, and returned halfway along the bar. It can be verified that values of *r*=*r*_{crit} result in stable wave reflection on the penalized boundary, whereas values of *r*>*r*_{crit} lead to destabilizing wave reflection. The particular unstable values of *r* in figures 1 and 2 have been selected such that the instabilities are visible at time *t*=150 s but not grown out of proportion as yet. Nevertheless, it is emphasized that as long as Δ*t*=Δ*t*_{crit} and *r*=*p*_{s}/*p*_{m}>*r*_{crit}, all instabilities are initiated when the wave front reaches the penalized boundary. This observation holds irrespective of the magnitude of the penalty factors *p*_{m} and *p*_{s}.

### (b) Bar problem—accuracy analysis

Next, the accuracy aspects of the bipenalty method are studied for the same bar problem as in §5*a*. The inertia penalty factor *p*_{m}=10^{4} throughout and different stiffness penalty factors *p*_{s} are taken. According to equation (3.13), this may affect the accuracy with which the constraint is modelled. The same combinations of mass matrix and time integration scheme are used, and Δ*t*=Δ*t*_{crit} in all cases.

Figure 3 shows the time history plots of displacement (figure 3*a*–*c*) and acceleration (figure 3*d*–*f*) of the constrained d.f. for the case where lumped mass and central difference time integration (*β*=0, ) is used. Three values for the penalty ratio *r* are taken, namely *r*=0 (figure 3*a*,*d*), (figure 3*b*,*e*) and *r*=2=*r*_{crit} (figure 3*c*,*f*). For this case, the Newmark parameter *β*=0 and equation (3.13) predicts that all three values of *r* should lead to the same convergence rates for the constrained acceleration. This seems to be verified in figure 3, in which the maximum error in the acceleration is about −4×10^{−4} m s^{−2}, irrespective of the value for *r*. However, the error in the acceleration is monotonic in case *r*=0, which leads to strongly accumulating error in the displacement, see figure 3*a*. On the other hand, taking *r*>0 leads to oscillating errors in the acceleration, the result of which is that the displacement error does not accumulate but instead remains bounded. These oscillations in the acceleration error are of higher frequency for larger *r* (i.e. for larger values of the stiffness penalty), which also reduces the magnitude of the displacement error.

Figure 4 shows similar results for the combination of consistent mass with Fox–Goodwin time integration (, ). The three values of the penalty ratio are here taken as, from top to bottom, *r*=0, and *r*=4=*r*_{crit}. Since now , there is a slight effect of *r* on the acceleration error, cf. equation (3.13): larger values for *r* imply larger values for the stiffness penalty and lower error in the acceleration. The relation between *r* and the displacement error is similar to that in figure 3. Zero stiffness penalty (*r*=0) leads to significant accumulation of the displacement error, whereas finite stiffness penalty (*r*>0) gives bounded displacement error. Moreover, the larger the stiffness penalty, the smaller the maximum error in the displacement.

### (c) Beam problem—stability analysis

Next, a horizontal Euler–Bernoulli beam of length *L*=50 m is considered. The cross-sectional area *A*=1 m^{2}, the second moment of area m^{4}, the Young’s modulus *E*=1 N m^{−2} and the mass density *ρ*=1 kg m^{−3}. An upward force *F*=0.01/Δ*t* N is applied at the left end of the beam during the first time step, and support conditions are applied via penalty functions at the right end of the beam—see below for specific details on the support conditions. The beam is subdivided into 50 equal-sized elements. The lumped mass matrix is used. Time integration is carried out with the central differences scheme and Δ*t*=Δ*t*_{crit}. The inertia penalty factor *p*_{m}=10^{4}.

In §4*b* critical penalty ratios were derived for different support conditions, namely *r*_{crit}=8 for a penalized deflection, *r*_{crit}=2 for a penalized rotation and for the combination of penalized deflection with penalized rotation. However, we will simply assume *r*_{crit}=2 for the latter case, since in practice *p*_{m}≫1. In figure 5 displacement profiles across the beam are plotted at time *t*=500 s. The effectiveness of using the bipenalty method on imposing the constraints can be verified qualitatively in figure 5*a*–*c*. In figure 5*d*–*f* the occurrence of instabilities can be verified in case the penalty ratio is chosen just 0.1 per cent larger than its critical value—this indicates that the derived values for the critical penalty ratio *r*_{crit} are crisp, even the approximative value of *r*_{crit} for the combination penalized deflection/ penalized rotation.

### (d) Plane strain problem—stability analysis

The final example uses square plane strain elements. A bar-type structure of height *H*=1 m and length *L*=10 m is modelled with 10 square four-noded finite elements. The Young’s modulus *E*=1 N m^{−2}, the mass density *ρ*=1 kg m^{−3} and a range of Poisson ratios were used (see below). Two leftward forces of 0.01 N each are applied on the two leftmost nodes from time *t*=0 onwards. The horizontal displacement of the top-right node is constrained using the bipenalty method. The consistent mass matrix is used with central difference time integration and Δ*t*=Δ*t*_{crit}.

In §4*c* an expression for the critical penalty factor *r*_{crit} was derived that is valid for a certain range of values for the Poisson ratio, namely in plane strain. As mentioned in §4*c*, an expression for *r*_{crit} could not be found for . It is therefore of interest to verify whether the expression for could be used as well for . Figure 6 shows the deformed mesh at time *t*=50 s, whereby the displacements have not been magnified. Various values were used for the Poisson ratio, namely *ν*∈{0.1,0.2,0.3,0.4}. Figure 6*a*–*d* shows results obtained with using a critical penalty ratio, whereas the results in figure 6*e*–*h* were found with a larger value for the penalty ratio, that is *r*=1.004*r*_{crit}. Taking *r*>*r*_{crit} leads to instabilities in case *ν*=0.3 and, more spectacularly, *ν*=0.4, which are within the range of the derived expression for *r*_{crit} as given in §4*c*. Taking *r*>*r*_{crit} does not seem to lead to instabilities for values of Poisson’s ratio , which is an indication that the suggested expression for *r*_{crit} is safe for .

## 6. Discussion and concluding remarks

In this paper we have introduced the bipenalty method, which consists of the simultaneous use of stiffness-type penalties and inertia-type penalties in the enforcement of constraints. The bipenalty method has significant advantages over the use of either stiffness penalties or inertia penalties. Compared with stiffness penalties, the bipenalty method can be tuned such that the critical time step of conditionally stable time integration schemes is not affected—this can lead to significant savings in computing time. Compared with inertia penalties, the bipenalty method is more accurate when penalty parameters of equal magnitude are taken.

The performance of the bipenalty method is assessed in terms of the ratio of stiffness penalty over inertia penalty, which is termed the ‘penalty ratio’. For reasons of numerical stability in conditionally stable time integration schemes, the penalty ratio is bounded from above. Closed form expressions for the so-called ‘critical penalty ratios’ have been derived for a range of element types, namely bar elements, beam elements and square plane stress/plane strain elements, including a distinction between lumped mass matrices and consistent mass matrices. These values are summarized in table 1.

For bar elements and beam elements, all the derived expressions for the critical penalty ratio have been confirmed by numerical experimentation (although not all individual results were presented). We have found that results are always stable if the penalty ratio is smaller than the critical penalty ratio, and results destabilize if the penalty ratio is larger than the critical penalty ratio. For the square plane stress/plane strain elements, two open issues must be mentioned:

— When using a consistent mass matrix, expressions for the critical penalty ratio could not be found for all values of the Poisson ratio. However, numerical experimentation shows that this does not seem to be a problem; the expressions for the critical penalty ratio are safe for the whole range of tested Poisson’s ratios.

— The theoretical derivations and the numerical simualtions were carried out with a single penalized d.f. However, in practice it is common that a single finite element has multiple constrained d.f. Preliminary derivations have not led to meaningful expressions for the critical penalty ratio, and further research is required to resolve this.

To summarize, the bipenalty method is superior to stiffness penalties (for reasons of stability; moreover, since time step restrictions can be avoided in the bipenalty method CPU times tend to be smaller than with stiffness penalties) and superior to inertia penalties (for reasons of accuracy). Although the term critical penalty ratio refers to the *stability* aspects of the bipenalty method, it was also found that the most *accurate* results are obtained if the penalty ratio is as close as possible to the critical penalty ratio—thus, the critical penalty ratio for stability is also the optimum penalty ratio for accuracy.

## Footnotes

- Received July 6, 2009.
- Accepted November 11, 2009.

- © 2009 The Royal Society