Vibration Method in Stability Analysis of Planar Constrained Elastica

The primary goal of the research in constrained elastica is to understand the behavior of a thin elastic strip under end thrust when it is subject to lateral constraints. It finds applications in a variety of practical problems, such as in compliant foil journal bearing, corrugated fiberboard, deep drilling, structural core sandwich panel, sheet forming, nonwoven fabrics manufacturing, and stent deployment procedure. By assuming small deformation, Feodosyev (1977) included the problem of a buckled beam constrained by a pair of parallel walls as an exercise for a university strength and material course. Vaillette and Adams (1983) derived a critical axial compressive force an infinitely long constrained elastica can support. Adams and Benson (1986) studied the post-buckling behavior of an elastic plate in a rigid channel. Chateau and Nguyen (1991) considered the effect of dry friction on the buckling of a constrained elastica. Adan et al. (1994) showed that when a column with initial imperfection positioned at a distance from a plane wall is subject to compression, contact zones may develop leading to buckling mode transition. Domokos et al. (1997), Holms et al. (1999), and Chai (1998, 2002) investigated the planar buckling patterns of an elastica constrained inside a pair of parallel plane walls. It was observed that both point contact and line contact with the constraint walls are possible. Kuru et al. (2000) studied the buckling behavior of drilling pipes in directional wells. Roman and Pocheau (1999, 2002) used an elastica model to investigate the post-buckling response of bilaterally constrained thin plates subject to a prescribed height reduction. Chen and Li (2007) and Lu and Chen (2008) studied the deformation of a planar elastica inside a circular channel with clearance. Denoel and Detournay (2011) proposed an Eulerian formulation of the constrained elastica problem.


Introduction
The primary goal of the research in constrained elastica is to understand the behavior of a thin elastic strip under end thrust when it is subject to lateral constraints. It finds applications in a variety of practical problems, such as in compliant foil journal bearing, corrugated fiberboard, deep drilling, structural core sandwich panel, sheet forming, nonwoven fabrics manufacturing, and stent deployment procedure. By assuming small deformation, Feodosyev (1977) included the problem of a buckled beam constrained by a pair of parallel walls as an exercise for a university strength and material course. Vaillette and Adams (1983) derived a critical axial compressive force an infinitely long constrained elastica can support. Adams and Benson (1986) studied the post-buckling behavior of an elastic plate in a rigid channel. Chateau and Nguyen (1991) considered the effect of dry friction on the buckling of a constrained elastica. Adan et al. (1994) showed that when a column with initial imperfection positioned at a distance from a plane wall is subject to compression, contact zones may develop leading to buckling mode transition. Domokos et al. (1997), Holms et al. (1999, and Chai (1998Chai ( , 2002 investigated the planar buckling patterns of an elastica constrained inside a pair of parallel plane walls. It was observed that both point contact and line contact with the constraint walls are possible. Kuru et al. (2000) studied the buckling behavior of drilling pipes in directional wells. Pocheau (1999, 2002) used an elastica model to investigate the post-buckling response of bilaterally constrained thin plates subject to a prescribed height reduction. Chen and Li (2007) and Lu and Chen (2008) studied the deformation of a planar elastica inside a circular channel with clearance. Denoel and Detournay (2011) proposed an Eulerian formulation of the constrained elastica problem.
The emphasis of these studies was placed on the static deformations of the constrained elastica. Very often, multiple equilibria under a specified set of loading condition are possible. Since only stable equilibrium configurations can exist in practice, there is a need to determine the stability of each of these equilibria in order to predict the behavior of the constrained elastica as the external load varies. For an unconstrained elastica, vibration method is commonly used to determine its stability; see Perkins (1990), Patricio et al. (1998), Santillan, et al. (2006), and Chen and Lin (2008). This conventional method, however, becomes useless in the case of constrained elastica.
The difficulty of the conventional vibration method arises from the existence of unilateral constraints. A unilateral constraint is capable of exerting compressive force onto the structure, but not tension. Mathematically, this type of constraints can be represented by a set of inequality equations. This poses challenges in determining the critical states of the loaded structure. In order to overcome this difficulty, the conventional stability analysis needs modification. In this chapter, we introduce a vibration method which is capable of determining the stability of a constrained elastica once the equilibrium configuration is known. The key of solving the vibration problem in constrained elastica is to take into account the sliding between the elastica and the space-fixed unilateral constraint during vibration.
In this chapter, we consider the vibration of an elastica constrained by a space-fixed point constraint. This particular constrained elastica problem is used to demonstrate the vibration method which is suitable to analyze the stability of a structure under unilateral constraint. In Section 2, we describe the studied problem in detail. In Section 3, we describe the static load-deflection relation. In Section 4, we introduce the theoretical formulation of the vibration method. In Section 5 an imperfect system when the point constraint is not at the mid-span is analyzed. In Section 6, several conclusions are summarized. Figure 1 shows an inextensible elastic strip with the right end fully clamped at a point B. On the left hand side there is a straight channel with an opening at point A. The distance between points A and B is L. Part of the strip is allowed to slide without friction and clearance inside the channel. A longitudinal pushing force A F is applied at the left end of the strip inside the channel causing it to buckle in the domain of interest between points A and B. An xy-coordinate system is fixed at point A. A point H fixed at position / 2 x L  and y h  prevents the elastica from deforming freely after the elastica contacts point H.

Problem description
The elastic strip is assumed to be straight and stress-free when A F =0. The effect of gravity is ignored. The strip is uniform in all mechanical properties along its length. The length and the shape of the elastica in the domain of interest vary as the pushing force A F increases. The boundary condition at point A may be called "partially clamped," by which we mean that the strip is allowed to slide freely through the opening A, while the lateral displacement and slope at A are fixed. The dashed and solid curves in Figure 1 represent two typical stages of the elastica deformation when A F increases beyond the buckling load. The dashed curve is a symmetric deformation pattern before the elastica contacts the point constraint.
The solid curve represents an asymmetric deformation after the elastica contacts the point H. Other deformation patterns may also exist, which will be discussed later.

Load-deflection relation
The equilibrium equation at any point (x,y) of the buckled strip between points A and B, as shown in Figure 1, can be written as M are the internal shear force and bending moment, respectively, provided by the partial clamp at A.  (positive when counter-clockwise) is the rotation angle of the strip at point (x,y). EI is the flexural rigidity of the elastic strip. s is the length of the strip measured from point A. For convenience we introduce the following dimensionless parameters (with asterisks):  is the mass per unit length of the elastica. t is time and  is a circular natural frequency, which will be discussed in the dynamic analysis later. After substituting the above relations into Equation (1), and dropping all the superposed asterisks thereafter for simplicity, we obtain the dimensionless equilibrium equation The method of static analysis can be found in . In this section we introduce several deformation patterns of the constrained elastica. All the physical quantities described henceforth are dimensionless.
The length of the elastica being pushed in through the opening is 1 l l    , where l is the dimensionless length of the elastica between points A to B. Figure 2 shows the relation between the edge thrust A F and the length increment l  . The height of the point constraint h is 0.03. The dashed and solid curves in this load-deflection diagram represent unstable and stable configurations, respectively. The method used in determining the stability of the static deformation will be described in detail in Section 4. The symmetric deformation before contact occurs is called deformation (1), whose locus starts at ( A F , l  ) = (1,0) and ends at (0.99668, 0.0022188). The slope of this load-deflection curve is slightly negative. After the middle point C of the elastica touches the point constraint H, the deformation pattern initially remains symmetric, called deformation (2). The load-deflection curve of deformation (2) starts at ( A F , l  )=(0.99668,0.0022188) and ends at (3.97314, 0.0026985). It is noted that the lower part of this load-deflection curve up to (1) (2.03268, 0.0023129) is solid and the upper part is dashed. At the point separating the solid and dashed parts, a symmetry-breaking bifurcation occurs and the elastica evolves to a pair of asymmetric deformations 4(a) and 4(b). As the pushing force continues to increase, it is natural to envision that the elastica may evolve to an "M" shape, i.e., there exist two inflection points in each half of the span, called deformation (3). The load-deflection curve corresponding to deformation (3) starts at ( A F , l  )=(3.97314, 0.0026985) and continues beyond the range of Figure 2. The slope of this curve is slightly negative.
At point ( A F , l  )=(2.03268, 0.0023129), the symmetric deformation (2) described previously may bifurcate to a pair of asymmetric deformations 4(a) and 4(b). Both of (4a) and (4b) have one inflection point on each of the half span separated by the point constraint. Both (4a) and (4b) start at ( A F , l  )=(2.03268, 0.0023129), while (4a) ends at (2.03545, 0.0028996) and (4b) ends at (2.02600 0.0028996). The slopes of curves (4a) and (4b) are positive and negative throughout, respectively. Stability analysis later indicates that deformation (4b) is unstable while (4a) is stable. As l  increases further, the pair of deformations (4a) and (4b) evolve to a pair of asymmetric deformations (5a) and (5b). Deformation (5a) has two inflection points on the left span and one inflection point on the right span. On the other hand, deformation (5b) has two inflection points on the right span and one inflection point on the left span. The slope of load-deflection curve of (5b) is negative throughout. On the other hand, the curve of (5a) is of convex shape with the top being at point ( A F , l  ) =(2.03554, 0.0031711), which is very close to the end of curve of deformation (4a). Deformation (5a) is stable first until the curve reaches the top at ( A F , l  ) =(2.03554, 0.0031711). At this critical point, the elastica will jump to another configuration. The load-deflection curves near the symmetry-breaking point are magnified and shown in the inset of Figure 2.
The theoretical load-deflection curves shown in Figure 2 give us a mental picture how the elastica evolves as the pushing force A F increases. First of all, the elastica remains still when A F is smaller than 1, the Euler buckling load. As soon as A F reaches 1, the elastica jumps to symmetric deformation (2) in contact with the point constraint. As A F increases, a symmetry-breaking bifurcation occurs and the elastica evolves to asymmetric deformation (4a) first and smoothly to (5a). As A F continues to increase up to a certain value, a second jump occurs. Following this jump the elastica will eventually settle to a self-contact configuration. This self-contact configuration requires a length increment l  over 8, which is well beyond the range of Figure 2. The above scenario has been verified experimentally in . In the next section we describe the vibration method used to determine the stability of the static deformations.

Lagrangian and Eulerian descriptions
As mentioned above, the deformation patterns discussed in Section 3 may not necessarily be stable. If the deformation is unstable, then it can not be realized in practice. In order to study the vibration and stability properties of the elastica, we first derive the equations of motion of a small element ds supported by the point constraint, as shown in Figure 3. The geometrical relations between x, y, and  are The balance of moment and forces in the x-and y-directions results in ( , ) x F s t and ( , ) y F s t are the internal forces in the x-and y-directions. The moment-curvature relation of the Euler-Bernoulli beam model is The readers are reminded that the functions x , y ,  , M , x F , and y F in Equations (3)-(8) are dimensionless and are all written explicitly in terms of s and t for clarity. These six equations can be called the Lagrangian version of the governing equations because a material element ds at location s is isolated as the free body. s may be called the Lagrangian coordinate of a point on the elastica. It is noted that s=0, 1 l , and l represent the material points at the left end, the contact point, and the right end, respectively, when the elastica is in equilibrium. During vibration, the elastica may "slide" on the point constraint. As a consequence, the contact point on the elastica may change from s= 1 l to is a small number. This change of contact point is reflected in Equations (6) and (7) H is the Heaviside step function. During vibration, the function ( , ) y F s t may be regarded as the superposition of ( ) ye F s in Equation (10) A variable with subscript "d" represents a small perturbation of its static counterpart with subscript "e." After defining a new variable  as Equation (11) can be rewritten as Apparently, ˆy F and y F are two different functions. It is noted that ( ) ye F  is the static solution as obtained from the static analysis, except that the independent variable s is replaced by  . Similarly, the other perturbed functions may be written as  defined in Equation (12)

Boundary conditions
The exact boundary conditions at the fixed end B are 1 ( , ) The boundary condition at the left end A is more complicated. We denote the material point on the strip right at the opening A of the channel as point A' when the elastica is in equilibrium, as shown in Figure 5(a). Since the strip is under a constant pushing force at the left end, A' will retreat into and protrude out of the channel when the elastica vibrates, as shown in Figures 5(b) and 5(c). We denote this small length of movement as Similarly, we can derive Finally, Equations (43) and (44) may be combined as The three equations (45)-(47) are the linearized boundary conditions at point A.

Constraint equations
When contact occurs, it is required that the elastica always passes through the point constraint. Mathematically, this condition can be written as We also require that the dynamic reactive force must be always normal to the elastica at the point constraint, or mathematically, After using Equations (17), (25)-(26) and neglecting higher-order terms, Equation (52) can be linearized to Equations (50), (51), and (53) are the three constraint equations.

Solution method
In summary, the six linearized differential equations (28)  , and the two dynamic constraint reactions xd R and yd R . It is noted that  in Equations (32)-(33) only appears in the form of 2  . Therefore, if the characteristic value 2  is positive, the corresponding mode is stable with natural frequency  . On the other hand, the equilibrium configuration is unstable if 2  is negative.
A shooting method is used to solve for the characteristic value 2  . Since the linear vibration mode shape is independent of the amplitude, we can set  Figure 6 shows the 2  of the first two modes as functions of the end force A F for deformation (2). The 2  of the first mode becomes negative when A F reaches 2.03268, at which the symmetry-breaking bifurcation occurs. It is noted that the 2  of the second mode becomes negative when A F reaches 3.97314. This happens to be the point at which deformation (3) begins to appear.
The first two mode shapes when A F =1.5 are depicted in Figure 7. The solid and dashed curves represent the static and the vibrating mode shapes of the constrained elastica, respectively. The first mode (a) is asymmetric and the second mode (b) is symmetric. To examine whether sliding at the contact point occurs we examine the values of 0d  and 1d  of each mode. For the asymmetric mode we found that the ratio 0d  / 1d  is approximately 1:400. This means that during vibration the protruding and retreating of the elastica at the channel opening A is negligible compared to the sliding at the contact point. In other words, the elastica length within the domain of interest is almost constant during vibration. This can also be observed from Figure 7(a). For the symmetric mode, we found that the ratio 0d  / 1d  is approximately 2:1. Therefore, sliding at the contact point still occurs. From Figure 7(b) we can observe that the lengths of the elastica on both sides of the point constraint increase (or decrease during the other half of the period) the same amount. Since the elastica is inextensible, the protruding 0d  at the channel opening has to be twice the amount of the sliding 1d  at the contact point. It is noted that a vibration analysis of a constrained structure will cause an erroneous result if sliding at the constraint is neglected. It is noted that the geometric conditions (50)-(51) at the point constraint are exact. Therefore, the vibrating elastica always passes through the point constraint no matter how large the vibration amplitude is. On the other hand, the boundary conditions at points A and B used in the calculation have been linearized from the exact boundary conditions. Therefore, the mode shapes do not necessarily satisfy the exact boundary conditions at points A and B. This may become obvious when the amplitude of vibration is increased dramatically.  negative along the loci of deformations (4b) and (5b), we conclude that deformations (4b) and (5b) are unstable. From these stability analyses, it is concluded that after the symmetrybreaking bifurcation of deformation (2), the elastica branches to deformation 4(a) and continue to deform along locus (5a). After the l  reaches 0.0031711, the elastica will jump to a remote self-contact configuration beyond the range of Figure 2.

Analysis of an imperfect system
The point constraint H in Figure 1 is at the middle between the two ends A and B. In practice, it is very difficult to place the point constraint accurately at the center. Instead, it is almost inevitable that the point constraint may be off the center somewhat. Figure 10  In Figure 11 we describe the change of the load-deflection relation when H  is increased from 0 (ideal case) to 0.01 and 0.05. The height h remains to be 0.03. Focus is placed on how the offset affects the symmetry-breaking bifurcation when deformation (2) branches into asymmetric deformations 4(a) and 4(b) in Figure 2. It is observed that the sharp corner at the bifurcation point degenerates into two smooth curves, called deformations 6(a) and 6(b) in Figure 11. Both deformations 6(a) and 6(b) are asymmetric. The tops of deformation 6(a) and 6(b) are to the left and right, respectively, of the point constraint H. Deformation 6(b) is always unstable. Deformation 6(a) for H  =0.01, on the other hand, is stable before point ( A F , l  )=(1.974, 0.00552). Figure 12 shows the 2 1  along the locus of deformation 6(a) when H  =0.01. It is shown that 2  becomes negative when l  =0.00552. For a larger H  =0.05, deformation 6(b) is stable throughout the range in Figure 11.
Inspecting Figure 11 reveals something unusual about the degeneration of the symmetrybreaking bifurcation due to the offset of the point constraint. For the ideal case with H  =0, part of upper branch (deformations 4(a) and (5a)) is stable until it reaches a maximum. On the other hand, the lower branch (deformations 4(b) and (5b)) is always unstable. When H  increases from 0 to 0.01, part of the lower branch is stable until it reaches a maximum at l  =0.00552. On the other hand, the upper branch is always unstable. It is not clear how the load-deflection curves evolve as H  varies.  In order to answer this question, we plot the load-deflection curves when H  varies with smaller increment. Figure 13(a) shows the load-deflection curves when H  =0, 5 5 10   , and 4 2 10   , respectively. For the case when H  =0, part of the upper branch is stable, while the whole lower branch is unstable, as in the inset of Figure 2. When H  increases by a small amount, the load-deflection curves degenerates into two branches veering away from each other. When H  = 5 5 10   , the stable range on the upper branch shrinks. When H  continues to increase to 4 2 10   , the stable range on the upper branch disappears altogether. For the lower branch, there exists a limit point. The locus with positive slope before the limit point is stable. Figure 13(b) shows another scenario when H  varies from 0 to 5 5 10    , and 4 2 10    . It is observed that for a negative H  , the sharp corner degenerates into two smooth curves crossing each other. The one emerging from the lower part has a stable range which ends at the peak of the curve. The other curve emerging from the top is unstable all the way.

Conclusions
In this chapter we introduce a vibration method which is suitable to analyze the stability of a constrained elastica. A planar elastica constrained by a space-fixed point constraint is used to demonstrate the method. Generally speaking, static analysis allows one to find all the possible equilibrium configurations of a constrained elastica. In order to predict how the elastica behaves in reality, the stability of these equilibrium configurations needs to be determined. The key of the vibration method is to take into account the sliding between the elastica and the unilateral constraint during vibration. In order to accomplish this, Eulerian coordinates are defined to specify the positions of the material points on the elastica. After transforming the governing equations and the boundary conditions from the Lagrangian description to the Eulerian one, the natural frequencies and the vibration mode shapes of the constrained elastica can be calculated. The vibration method is applied to an elastica constrained by a point constraint in this chapter. The same principles can be extended to other similar problems as well, for instance; multiple point constraints  and plane constraints .