## Abstract

The work presented here justifies the use of restraints with large positive and negative stiffness values to asymptotically model geometric constraints of a structure in linear structural analysis. Replacing constraints by very stiff restraints improves the versatility of variational methods such as the Rayleigh–Ritz method as the limitation on the choice of functions is removed. Based on recently published theorems on the existence of natural frequencies of systems with artificial restraints of positive and negative stiffness and their convergence towards the frequencies of the corresponding constrained systems, a proof is given to show that the displacement of a constrained structure caused by any action along its direction is approached and bracketed by the displacement of systems with artificial restraints of positive and negative stiffness (positive and negative penalty functions) as the magnitude of stiffness is increased. A procedure for the practical use of positive and negative penalty functions in static analysis of linear structures is also proposed.

## 1. Introduction

In applying the Rayleigh–Ritz method to determine the natural frequencies of continuous systems, it has been a common practice to use restraints with large positive stiffness to asymptotically model geometric constraints (Courant 1943; Gorman 1989; Cheng & Nicolas 1992; Yuan & Dickinson 1992, 1994; Lee & Ng 1994; Amabili & Garziera 1999). This removes a limitation on the choice of shape functions as these functions need not satisfy the constraints individually. This is also achievable using the Lagrangian multiplier method in which the individual functions are permitted to violate the constraints which are satisfied by the series as a whole (Ilanko & Dickinson 1999). However, the Lagrangian multiplier method requires additional programming effort since the final matrix formulation includes two types of equations: the modified minimization equations and the constraint equations. This makes the programming more time consuming than for the regular Rayleigh–Ritz method which results in only minimization equations. Furthermore, in the Lagrangian multiplier method, continuous constraints such as line supports have to be approximated by a series of point constraints. This is not necessary when using the Rayleigh–Ritz method with artificial restraints as the strain energy in continuous restraints can be included in an integral form.

Due to the above reasons, the use of restraints with large stiffness instead of constraints is becoming increasingly common when applying the Rayleigh–Ritz method for structural analysis, particularly for calculating natural frequencies. A more general version of this method, namely the use of penalty functions, is widely used to solve a vast range of constrained variational and optimization problems in applications such as the Rayleigh–Ritz method, Finite Elements and Elements Free Galerkin's Method (Zienkiewicz 1974, 1977; Gavete *et al*. 2000; Pannachet & Askes 2000). Large penalty parameters which correspond to restraints with high stiffness are multiplied by quadratic error functions that represent the resultant violation of the constraints, and added to the function being minimized or optimized. However, the problem of determining a suitable penalty parameter that is large enough to minimize the error and small enough to avoid computational problems, first noted by Courant (1943), has remained unresolved for several decades. The error decreases with the size of the penalty parameter, and in many cases a suitable penalty parameter is found by repeating the calculations for increasing values of the penalty parameter until numerical convergence is reached (Zienkiewicz 1974). Another problem with this method is the absence of any information on the magnitude of error caused by the violation of constraints.

This paper justifies a method to overcome the above problems, proposed in a recent publication (Ilanko 2002*a*) where the Rayleigh–Ritz method was applied to determine the deflection of a propped cantilever. It was shown that by using a combination of positive and negative values for the stiffness of an asymptotic restraint that replaced the prop, any error due to the violation of a geometric constraint may be kept within any desired tolerance. This was only a demonstration of the method and lacked proof. However, justification for using positive and negative stiffness in determining natural frequencies of constrained systems is now available in the form of two theorems (Ilanko 2002*b*). This has also been proven to be true in the determination of critical loads of linear elastic structures (Ilanko 2003). This paper extends the justification to static analysis of linear elastic structures. Three new theorems are presented, proving that the displacements of asymptotic models with artificial restraints having positive and negative stiffness (penalty parameters) asymptotically approach the displacement of a constrained structure as the magnitude of the stiffness approaches infinity, and that if the displacement is caused by a force or moment applied along its direction then the displacement of the asymptotic models also bound the displacement of the constrained structure.

## 2. Theoretical derivations

The following notation is used in the derivations. Uppercase letters *A*–*C* are used to name structural systems. Systems with added constraints are denoted by uppercase letters with tilde and have subscripts, which give the number of constraints added. For example, denotes a constrained structure obtained by adding *h* constraints to System *A*. The subscripts are given in brackets (for example, ), if the constraints are not associated with an inertial element possessing mass or moment of inertia. In static analysis whether a structure possesses inertia or not is irrelevant. However, some intermediate steps in the derivations require treatment of an analogous vibratory system and this distinction is useful. Structures with added artificial restraints are denoted by uppercase letters with subscripts which give the number of added restraints. Once again the subscripts may be given in brackets if they are not associated with a mass or inertia.

To explain the derivations an illustrative example will be used but the applicability of the statements are not confined to this example. The example chosen is a propped cantilever with a clamped end and two simple supports (figure 1*a*). This is labelled as it may be considered as having been derived from a cantilever which is labelled System *A* (figure 1*b*) by applying two constraints (*h*=2) that ensure that the displacement at each prop is zero. The mass of the structure may be disregarded in static analysis and therefore the subscript *h* is given in brackets. The purpose of the following derivations is to show that the displacement at any given point in due to any action is approached by the displacement of an asymptotic model *A*_{(h)} in which the *h* constraints are replaced by restraints of large positive and negative stiffness (figure 1*c*). To distinguish asymptotic models with positive and negative stiffness, the labels *A*_{(h+)} and *A*_{(h−)} are used here.

For the following derivations, initially the asymptotic model *A*_{(h)} will be modified by adding *h*+1 inertial elements along the *h* coordinates corresponding to the artificial restraints and the coordinate along which the action is applied, resulting in a vibratory system *B*_{h}. The inertial element may be a particle possessing mass if the restraint is a translational type, or if the applied action is a force. A rigid body possessing a moment of inertia may be used as an inertial element if the restraint is a rotational type or if the applied action is a moment. For the example used as illustration, this modified system is shown in figure 2*a* where three particles possessing mass are attached at coordinates corresponding to the two artificial lateral restraints and the applied force. Let the magnitude of the mass of the particle at the right end corresponding to the applied force *P* be *m*. The use of mass or moment of inertia is only for the purpose of proof and as explained later the actual system need not possess any concentrated mass or inertia at any point or may have distributed mass throughout. In the following section, the displacement corresponding to an action (at the same point in the same sense) will be shown to be bounded and approached by the displacement of an asymptotic model with artificial restraints of positive or negative stiffness.

### (a) Convergence and bounds of displacement due to an action along its coordinate

For finite values of stiffness, System *B*_{h+} is an *h*+1 degree of freedom system. For very small values of negative stiffness, System *B*_{h−} is also an *h*+1 degree of freedom system. However, as explained in Ilanko (2002*b*), as the magnitude of the negative stiffness increases, one by one, the first *h* natural frequencies would become imaginary and for very large negative stiffness values the structure would become a single degree of freedom system. If the inertial elements associated with the artificial restraints were removed, the resulting system *B*_{(h)} is a single degree of freedom system as shown in figure 2*b*. The derivations of the existence and convergence theorems in Ilanko (2002*b*) implicitly assume that the constraints are associated with inertial elements and that the application of each constraint decreases the degrees of freedom by one. Further work (Ilanko submitted) justifies asymptotic modelling even for cases where the constraints do not change the degrees of freedom as in the case of system (figure 2*c*) which may be obtained by setting the stiffness of the artificial restraints in system *B*_{(h)} to very large positive and negative values.

From the existence and convergence theorems (Ilanko 2002*b*), for positive stiffness values, as the stiffness approaches infinity, the highest *h* natural frequencies of System *B*_{h+} become infinite and the remaining natural frequency (lowest) approaches that of the constrained system from below. For system *B*_{(h)} also, similar inequality and convergence relationships for the fundamental natural frequency have been established (Ilanko submitted). For the purpose of the present derivations, let the stiffness values of all artificial restraints be the same, i.e. *k*_{i}=*k* for *i*=1,2,…,*h*.

The theorems in Ilanko (submitted) give,(2.1a)(2.1b)

The notations used in this paper are slightly different from that used by the author in his previous work (Ilanko 2002*b*, submitted) as different systems (*A*_{h+}, *A*_{h−}, *B*_{h−}, etc.) need to be treated together here. The symbols , , , and represent the *i*th natural frequencies of systems *B*_{(h+)}, *B*_{h+}, *B*_{(h−)}, *B*_{h−} and , respectively.

From the existence and convergence theorems (Ilanko 2002*b*), for negative values of stiffness as the magnitude of stiffness approaches infinity, the highest frequency of System *B*_{h−} would approach the natural frequency of the constrained system from above while all other *h* natural frequencies would become imaginary.(2.2a)(2.2b)

In single degree of freedom systems where the constraints and artificial restraints are not associated with any inertial term, there can only be one natural frequency. However, an inequality relationship similar to the above has been shown to hold, if the magnitude of artificial stiffness parameter used is greater than the critical stiffness parameter with the highest magnitude (Ilanko submitted). The critical states correspond to the poles in a frequency squared versus stiffness parameter plot and the number of such critical states can be shown to be equal to the number of artificial restraints *h*. The inequality and convergence equations for systems *B*_{(h−)} and may be written as:(2.2c)(2.2d)

Here, *k*_{crit,h} is the critical stiffness parameter with the highest magnitude.

From equations (2.1*a*) and (2.2*c*),(2.3)

Equation (2.3) states that the natural frequency of System is bounded by the fundamental natural frequency of *B*_{(h+)} and *B*_{(h−)}, and from equations (2.1*b*) and (2.2*d*), it may be seen that as the magnitude of stiffness approaches infinity the difference between these frequencies becomes zero. This means the natural frequency of the constrained system may be delimited to any desired accuracy by calculating the natural frequencies of systems *B*_{(h+)} and *B*_{(h−)}.

From first principles the natural frequency of the single degree of freedom system is(2.4)where is the structural stiffness coefficient of the constrained system defined as the force required to cause a unit displacement along the coordinate of motion of *m*. This stiffness coefficient is independent of the mass, and is the same as that of system . If the motion of the mass is along the applied action, then the displacement along the applied action would be given by(2.5)

Using equation (2.4),(2.6)

Equation (2.6) may also be obtained from d'Alembert's principle. It should be noted here that the ‘mass’ has been introduced only to provide an inertial resistance equivalent to that of the elastic stiffness. Similarly, the displacement of the asymptotic system *A*_{(h)} in figure 1*c* may be related to the natural frequencies of the corresponding vibratory systems. Let the displacement of the systems *A*_{(h+)} and *A*_{(h−)} (see figure 1*c*) along any prescribed action *P* be *δ*^{(+)} and *δ*^{(−)}, respectively. Then,(2.7)and(2.8)

From equations (2.3), (2.6), (2.7) and (2.8),

For all positive *k* and for all negative *k* where |*k*|>|*k*_{crit,h}|,(2.9)

Furthermore, from equations (2.1*b*), (2.2*b*), (2.6), (2.7) and (2.8),(2.10a,b)

The force *P* and mass *m* used in equations (2.4)–(2.9) may be replaced with a moment term and a moment of inertia term, in which case the displacements , *δ*^{(+)} and *δ*^{(−)} would represent rotations instead of translations. Equations (2.9) and (2.10*a*,*b*) may be stated in the form of the following theorems:

*The displacement of a constrained structure subject to any prescribed action,* *along the coordinate of such action, is greater than or equal to the corresponding displacement of an asymptotic model of the structure in which the constraints are replaced with restraints of negative stiffness provided the magnitude of the negative stiffness is greater than the highest critical stiffness,and less than or equal to the corresponding displacement of an asymptotic model of the structure where the constraints are replaced with restraints of positive stiffness. The action may be a force or a moment. If the action is a force then the displacement would be a translation, and if it is a moment the displacement would be a rotation.*

*As the magnitude of stiffness of artificial restraints approaches infinity the displacement due to any applied action along its direction,* *obtained using the asymptotic model, would approach the displacement of the constrained structure, from above in the case of positive stiffness, and from below in the case of negative stiffness.*

### (b) Convergence of the displacement profile of an asymptotic model

The above theorems only state that the displacement of a constrained structure along the coordinate of an applied action is bounded and approached by the displacement of the asymptotic models as the magnitude of stiffness of the artificial restraints takes very large values. They do not prove that the displacement along any other coordinate is also bounded and approached in this manner. The following proof shows that for linear, elastic constrained structures, the displacement at any point or along any coordinate due to any given load is approached by the displacements of the corresponding asymptotic models with positive and negative stiffness. While numerical results show that for very large stiffness values the displacement profile of a constrained structure may also be bounded by those of the asymptotic models, the author is unable to prove this at present.

Suppose the displacement of the constrained structure along a given coordinate *q*_{j} due to an action corresponding to another coordinate *q*_{i} is required. For illustrative purpose, let the required displacement be a lateral translation at *x*=*x*_{j} as shown in figure 3*a*. For this proof, let be the flexibility coefficient of the constrained system giving the displacement along *q*_{i} due to a unit action along *q*_{j} and the corresponding flexibility coefficients from the asymptotic models with positive and negative stiffness be and .

Applying theorem 2.2 to the situation where the action and displacement are along the same coordinate,(2.11a,b)

Let us introduce two unit masses, *m*_{i}=1 and *m*_{j}=1 along *q*_{i} and *q*_{j}, respectively. From the convergence theorems (Ilanko 2002*b*, submitted), the natural frequencies of the resulting two degree of freedom system *C*_{(h)} (either *C*_{(h+)} or *C*_{(h−)}) shown in figure 3*b* asymptotically approach the natural frequencies of the constrained system shown in figure 3*c* as follows.(2.12a,b)

The direction of approach depends on the sign of the artificial stiffness (Ilanko submitted). The following relationships between the flexibility coefficients and the fundamental natural frequencies of the two degree of freedom constrained system and its asymptotic models may be obtained by rearranging the frequency equations for these cases which are readily obtainable from the equations of motion by following the standard procedure given in textbooks (Tse *et al*. 1978).(2.13)(2.14)

From equations (2.11) and (2.12),(2.15)

Substituting equations (2.13) and (2.14) into the above equation gives:(2.16)

Since the choice of coordinates *q*_{i} and *q*_{j} is arbitrary, we may conclude that for large values of stiffness, the displacement along any coordinate at any point in the asymptotic model approaches the displacement of the constrained system. This means the displacement profile of the constrained structure may be found using the asymptotic models. Using the principle of superposition, this convergence must also hold for any given set of loading. This may be stated in the form of the following theorem:

*The displacement of a linear, elastic, constrained structure at any given point subject to any given loading is approached by the displacement of an asymptotic model in which the constraints are replaced by restraints with very large positive or negative stiffness, as the magnitude of the stiffness becomes infinite.*

The above derivations only assure convergence of the structural displacement profile towards that of the constrained system, and no proof is yet available to confirm the bracketing of the displacement profile of the constrained solution by those obtained for the asymptotic models with positive and negative stiffness values. The latter is only assured in the case of a specific displacement corresponding to the coordinate of the applied action (theorem 2.2).

Due to linearity, since the displacement profile converges, the stress, strain and internal forces and moments at any point or section in the asymptotic models are also expected to approach that of the constrained structure in the same manner. It should be noted that in using a discretization method such as the Rayleigh–Ritz method, any error in the calculated displacement profile will result in greater error in stress, strain and internal actions such as bending moment or shearing force as they involve differentiation of the displacement. However, as explained later, numerical results for an illustrative example show that any error due to asymptotic modelling may be controlled by using a combination of positive and negative penalty parameters.

In practical applications it is not necessary to find the natural frequency in order to determine the displacements as these can be done by solving the equations of equilibrium or by minimizing the total potential energy. However, it should be borne in mind that in arriving at the inequality statement in equation (2.9), the magnitude of artificial stiffness needs to be greater than the highest critical stiffness parameter to ensure that all the poles in the frequency squared versus stiffness curve have been passed. In static analysis the displacement *δ*^{(−)} goes through *h* critical states, and in order to ensure that the inequality relationship given by equation (2.2*c*) holds, the magnitude of negative stiffness used must be greater than the highest critical stiffness value as explained in the illustrative example in the next section.

## 3. Illustrative example

Results for a cantilever with one prop have already been reported (Ilanko 2002*a*). To get a more general perspective that will aid in drawing some guidelines for asymptotic modelling, consider the cantilever subject to two artificial restraints (*h*=2) as an illustrative example (figure 1*a*). The artificial restraints of stiffness *k* are applied at distances *b*_{1} and *b*_{2} from the clamped end. A force *P* is applied at the right end where *x*=*L*. The exact solution to this problem was obtained by direct integration of the moment-curvature equation from the Euler–Bernoulli beam bending theory *M*=*EIf*″, where *M* is the bending moment, *EI* is the flexural rigidity and the second derivative of deflection *f* is the curvature. In addition approximate results were generated using the Rayleigh–Ritz method. Since the application of asymptotic modelling is more likely to be useful in a Rayleigh–Ritz formulation, the basic steps in the derivations are presented only for this case. However, for comparison the exact results are also given in tabular and graphical forms.

The following simple power series satisfies the boundary conditions at the clamped end but permits displacements at the restraints.(3.1)

The potential energy of the beam is given by(3.2)

The Rayleigh–Ritz minimization equations, which in this case are statements of minimum potential energy (∂*V*/∂*a*_{i}=0), result in the following equation:(3.3)where(3.3a,b)

The above equation was solved for *b*_{1}=0.25*L* and *b*_{2}=0.5*L*, for various values of the stiffness parameter, and the displacement profile found by substituting {*a*} into equation (3.1). Ten terms (*n*=10) in the series for *f* were found to be sufficient to obtain results that agree with exact results to at least two significant figures.

For convenience the results are presented in terms of non-dimensional parameters, a non-dimensional deflection profile *Δ*(*x*) and a non-dimensional stiffness coefficient of the artificial restraint *α* which are defined as follows:(3.4a)(3.4b)

From hereon this stiffness coefficient will be referred to by its generic name: the penalty parameter. The variation of the non-dimensional tip deflection *Δ*(*L*) with the penalty parameter is shown in figure 4. Some numerical results are given in table 1.

From table 1 and figure 4, it may be seen that the tip deflection for the constrained system is bracketed by the results obtained using asymptotic models with positive and negative penalty parameters of large magnitude. As the magnitude of the penalty parameter increases the difference between the constrained solution and the asymptotic solutions decreases as predicted. There are two poles in figure 4 corresponding to the critical states associated with the penalty parameter in the negative range. The bounding nature of the solution from the positive and negative asymptotic models holds for all positive penalty values and for negative penalty values that are greater in magnitude than the highest critical penalty parameter. Therefore, in order to ensure that the required constrained system solution is bounded by the results obtained from asymptotic models it is necessary to use negative penalty parameters that are greater in magnitude than the highest critical penalty parameter. The critical states occur for negative penalty parameters and the number of critical states would be equal to the number of artificial restraints. These critical values may be obtained numerically by gradually increasing the magnitude of the negative penalty value until *h* poles are passed. Alternatively, they may be found by solving an eigenvalue problem as follows.

Equation (3.3) may be reformulated to separate the effect of the restraints.(3.5)where(3.6)(3.7)

The system is critical if(3.8)

This is a standard eigenvalue problem that could be readily solved. For the present problem, the Rayleigh–Ritz method with 10 terms (*n*=10) gives the critical values of *α* as −21.82 and −970.64 whereas the corresponding exact values obtained by direct integration method are −21.82 and −965.61, respectively.

Results for the deflection profile obtained using both the Rayleigh–Ritz method and the exact direct integration method are presented in figure 5 for *α*=−10 000, −2000, 2000 and 10 000 together with exact results for the constrained system (*α*=∞). Some numerical results are given in table 2. These results show that as the magnitude of the penalty parameter is increased the deflection profiles from the asymptotic models approach that of the constrained system, irrespective of the sign of the stiffness and the method used. The discrepancy between the Rayleigh–Ritz and exact results is not noticeable in the figure.

The deflection-penalty parameter plots were generated for the deflection at every 20th of the span (i.e. for *x*=0, 0.05*L*, 0.1*L*,…,*L*). For 0.5*L*≤*x*≤*L*, the plots were similar to figure 4. For 0<*x*<0.5*L*, a different type of plot was obtained. A typical plot for this range is shown in figure 6 which was obtained for *x*=0.2*L*. The characteristic difference between these types of plots is that for 0.5*L*≤*x*≤*L* the slope of the deflection-penalty parameter curve is always negative, except at the poles where it is undefined. This is only guaranteed for the displacement corresponding to the load (in the present example *Δ*(*L*)), because from equations (2.7) and (2.8) this displacement is inversely proportional to the square of the frequency which can only increase with stiffness which is represented by the penalty parameter (Rayleigh 1945).

In figure 6, the deflection-penalty parameter curve has two stationary points and changes slope twice. An interesting point to note is that one of the stationary points is in the positive penalty range. Checking the accuracy of a penalized solution in the proximity of this point by doing a simple numerical convergence test, as has been the long standing practice in penalty function applications (Zienkiewicz 1974, 1977; Gavete *et al*. 2000; Pannachet & Askes 2000) may lead to incorrect results. For penalty values of 100, 200, 300, 400 and 500, the displacement parameters for *x*=0.2*L* were found to be 0.009, −0.0006, −0.0010, −0.0012 and −0.0012, respectively. These results appear to indicate (incorrectly) that the displacement parameter has converged at −0.0012. However, with further increase in the value of the penalty parameter, the displacement parameter starts to increase and finally converges to 0.0003 for penalty parameters above 10^{5}. A reassuring confirmation is that the same result is obtained for a penalty parameter of −10^{5} and also from direct integration method for the fully constrained structure. The convergence of the solution for positive and negative penalty parameters is from opposite directions as shown in figure 6. This means the deflection of the constrained system is not only approached by the deflection of the asymptotic models as predicted by theorem 2.3, but also bracketed for very large values of the penalty parameter. This was observed to hold also for the deflection at all other points for which results were generated. The author is unable to find a simple explanation or proof for this. Even if such a proof is derivable, the requirement that the magnitude of the negative penalty parameter is greater than the highest critical stiffness is not a sufficient condition to assure bracketing, except when determining the deflection due to a load along its direction, as can be seen from figure 6.

The distribution of the bending moment in the propped cantilever and its asymptotic models were also calculated using both, the exact method and the Rayleigh–Ritz method. Numerical results for the non-dimensional bending moment defined by *ζ*=*M*/(*PL*) are presented in table 3. It may be noted that for large magnitudes of penalty parameter, the exact results for the asymptotic models converge to the constrained solution from opposite directions. Although the agreement between the bending moments obtained using the Rayleigh–Ritz method and the exact method is not good, the results for positively and negatively restrained systems approach each other for large magnitudes of the penalty parameter. The significant discrepancy between the Rayleigh–Ritz solution and the exact solution may be attributed to the presence of the discontinuities in the shear force due to the restraints (or constraints) at *x*=0.25*L* and *x*=0.5*L*. It is possible that the 10 polynomial terms used are not sufficient to allow for these discontinuities.

When using penalty functions or artificial restraints to enforce geometric constraints in static analysis of structures, the procedure in the following section is proposed. In general this is likely to help to control any error due to the violation of the constraints and therefore give more reliable results. In calculating the displacement due to an action along its direction, these steps will ensure that such errors are kept within any desired tolerance.

## 4. Procedure for applying asymptotic modelling in linear structural analysis

Determine the required displacement or stress due to the applied loading on an asymptotic model where all constraints are replaced with artificial restraints of any positive stiffness value (penalty parameter).

Find the critical penalty parameter with the highest magnitude

*α*_{c,h}by solving the eigenvalue problem of the type given by equation (3.8) or by tracing the effect of the negative stiffness (penalty parameter).Use a negative penalty parameter of magnitude greater than

*α*_{c,h}and determine the required displacement or stress and compare it with the result obtained in step 1.Repeat steps 1 and 3 for slightly different magnitudes of the penalty parameter to see if positive and negative values lead to convergence from opposite sides. This step is not necessary if the required displacement is due to a single action (force or moment) along its direction as bounding is assured for this case from theorem 2.2.

If convergence is from different sides, and the difference between the results using the positive and negative penalty values is smaller than the desired tolerance, the computation may be terminated.

Repeat the above steps for increasing magnitudes of positive and negative penalty values until the difference between the results for the positive and negative penalty values reduces to an acceptable tolerance and convergence from opposite directions is observed.

To explain the steps, extracts from the illustrative example are given below.

If the non-dimensional displacement of the propped cantilever in our example at

*x*=0.2*L*is sought to within 0.0005, for*α*=300,*Δ*(0.2*L*)=−0.0010.The critical penalty parameter with the highest magnitude is about −970.

For

*α*=−2000,*Δ*(0.2*L*)=0.0035.For the problem at hand, this is necessary since the load is applied at

*x*=*L*, and the displacement is sought at*x*=0.2*L*. For*α*=−3000,*Δ*(0.2*L*)=0.0019, and for*α*=400,*Δ*(0.2*L*)=−0.0012. Using the results from steps 1, 3 and 4, the slope of the displacement-penalty parameters is positive for negative penalty values selected and negative for the positive penalty values. The approach is not from opposite sides.The computation may not be terminated because the approach is not from opposite sides and the difference between the results is greater than the tolerance.

Repeating the calculations for

*α*=10^{4}and 10^{5}give the displacement parameter as 0.0026 and 0.0028, respectively, giving a positive slope. For*α*=−10^{4}and −10^{5}the corresponding results are 0.0063 and 0.0032, respectively. The slope is now again positive. The constrained solution is now bracketed between 0.0028 and 0.0032. The difference between these values is less than the tolerance. Calculations may now be terminated.

## 5. Concluding remarks

It has been shown that for any linear elastic structure subject to one or more geometric constraints, the displacements due to any given loading may be found by using an asymptotic model where the constraints are replaced with elastic restraints of very large positive or negative stiffness (penalty parameter). It has been proven that the displacement profile of the asymptotic models would approach that of the constrained structure, as the magnitude of the penalty parameter approaches infinity. This approach will be from opposite directions if the displacement is due to a force or moment applied along its direction. A procedure to ensure that the error due to the violation of the constraints is kept within any desired tolerance has been presented.

## Acknowledgments

The work presented here commenced during the author's study leave at the Department of Mechanics and Ocean Engineering, TUHH, Germany. His visiting Professorship was funded jointly by the host institute and the German Academic Exchange Programme (DAAD) for which he is grateful to both organizations. He also wishes to thank his colleagues, Emeritus Professor Harry McCallion for his comments and encouragement, Dr John Smaill for his help in preparing the final version of the diagrams, and the anonymous referees for their suggestions which have helped to enhance this article.

## Footnotes

- Received September 17, 2004.
- Accepted May 18, 2005.

- © 2005 The Royal Society