- Email: [email protected]
Deformation and stability of a spatial elastica under a midpoint force
European Journal of Mechanics A/Solids 54 (2015) 84e93
Contents lists available at ScienceDirect
European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol
Deformation and stability of a spatial elastica under a midpoint force Jen-San Chen*, Re-Ming Chen Department of Mechanical Engineering, National Taiwan University, Taipei, 10617, Taiwan
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 December 2014 Accepted 7 June 2015 Available online 20 June 2015
In this paper we study the deformation and stability of a pinnedepinned buckled beam under the action of a concentrated force at the midpoint. Focus is placed on the snap-through phenomenon, which may take place in a plane or three-dimensional space. We first find the equilibrium configurations by using shooting method. Elastica model is adopted to take into account exact geometry in large deformation. As expected, multiple solutions may exist for a specified set of loading parameters. Vibration method is then employed to determine the stability of the equilibrium solutions. Through these analyses the deformation sequence as the midpoint force increases quasi-statically can be predicted. It is found that the deformation sequence of the elastica is determined by two parameters; (1) the distance between the two ends of the buckled beam, and (2) the bending stiffness ratio of the cross section. Ten different deformation patterns can be identified according to four characteristics; the deformations before, after, and during the jump, and the type of critical point at the jump. A metallic wire with circular cross section is used to verify the predicted deformation sequence. It is concluded that for the specific specimen in the demonstration if one wishes to design an elastica capable of only plane deformation in all range of end distance, then the bending stiffness ratio has to be greater than 28. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Spatial elastica Midpoint force Snap-through
1. Introduction In the design of miniature mechanisms, such as in MEMS, it is impractical to use conventional rigid-body joint pairs due to space constraint. More often, the motion transmission is accomplished by the deformation of a single flexible beam. In particular, the snapthrough phenomenon of a curved beam has been used for microswitches (Maurini et al., 2007; Zang et al., 2007; Krylov et al., 2008; Medina et al., 2012; Ouakad and Younis, 2014). Two kinds of curved beams under midpoint forces can be found in the literature. In the first type, the curved beam is stress free when it is in curved configuration, see Fung and Kaplan (1952), Schreyer (1972), Plaut (1978), Chen et al. (2009), Virgin et al. (2014). In the second type the beam, which is stress free when it is straight, is buckled into a curved shape by edge thrust. The buckled beam is then loaded by a midpoint force, see Seide (1984), Pippard (1990), Patricio et al. (1998), Vangbo (1998), Kublanov and Bottega (1995), Pinto and Goncalves (2000), Cazottes et al. (2009), Chen and Hung (2011, 2012), Chen and Tsao (2013).
* Corresponding author. E-mail address: [email protected] (J.-S. Chen). http://dx.doi.org/10.1016/j.euromechsol.2015.06.006 0997-7538/© 2015 Elsevier Masson SAS. All rights reserved.
The snap-through phenomena observed in these studies are limited in plane deformation. Pippard (1990) noted in his experiment that a sway mode of instability normal to the plane of the strip tends to occur. Intuitively, if the ratio between the width and the thickness of the strip is sufficiently large, the strip should behave like a planar elastica. This raises the question what the minimum ratio is in order for the strip to deform only in a plane. Furthermore, what will happen if three-dimensional deformation does occur? In this paper we pursue this interesting phenomenon by considering an elastic rod capable of out-of-plane deformation. Elastica model is adopted in the analysis. The elastic rod considered in this paper is stress free when it is straight. The two ends of the rod are pinned in space after it is buckled into a curved shape. A force is applied at the midpoint in the plane of the rod. The two principal moments of inertia of the cross section may be different. If the bending stiffness in one of the principal directions is much larger than the other, it is expected that the deformation will be restricted in a plane, as described in the works cited previously. If the two bending stiffness are comparable, on the other hand, spatial deformation may occur. It is the objective of this paper to identify all these deformation patterns. In Section 2 we establish the equations of motion of the loaded rod. In Section 3 we conduct a static analysis to find the equilibrium configurations. In Section 4 we study the vibration characteristics
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
of the loaded elastica. The stability of the loaded elastica can then be determined from the calculated natural frequencies. In Section 5 we present in detail a numerical example showing the loaddeflection diagram and the frequency spectrum. A pinned metallic wire with circular cross section is used to verify the predicted deformation sequence as the midpoint load increases quasistatically. In Section 6 we show a phase diagram using the distance between the two ends and the stiffness ratio between the two principal moments as two parameters. In the diagram, ten different deformation behaviors are identified. In Section 7 several conclusions are summarized.
2. Equations of motion We consider a uniform, inextensible, and unshearable elastic rod with length L0, cross section area A, Young's modulus E, shear modulus G, and mass density per unit volume m. The area moment of inertia in the two principal directions of the cross section are I1 and I2, where I1 I2. The rod is initially straight and stress free. The rod is first buckled in the plane containing the principal direction with I2 by pushing the two ends closer to a distance L*. The two ends are then fixed in space with pin joints, as shown in Fig. 1(a). A space-fixed x*y*z*-coordinate system with origin attached to the left end O is chosen to describe the geometry of the rod. The x*-axis is pointing to the fixed end on the right. The inertial orthonormal frame {ex,ey,ez} is associated with the x*y*z*-coordinate system. The location of a material point on the neutral axis of the deformed rod is denoted by arc length s* measured from end O. We consider the case when the buckled beam is loaded by point force Q* at the midpoint s* ¼ L0/2 in the y* direction. A body-fixed right-handed orthonormal frame {d1,d2,d3} is chosen in such a manner that vector d3 is in the direction of the local tangent of the deformed neutral axis. Vectors d1 and d2 are in the normal cross section of the rod and rotate along with the cross section. Fig. 1(b) shows the case when the cross section is of elliptic
85
shape with semi-major axis b and semi-minor axis a. When the rod is in the unstressed straight state, the frame {d1,d2,d3} coincides with the frame {ex,ey,ez} in such a manner that d3 ≡ ex, d1 ≡ ey, d2 ≡ ez. The vectors di(s*,t*) (i ¼ 1, 2, 3), where t* is time, can be expanded as
di ¼ dix ex þ diy ey þ diz ez
ði ¼ 1; 2; 3Þ
(1)
It is assumed that the cross section of the rod remains plane and normal to the neutral axis after deformation. The rotation of the cross section can be defined by the vector d1. The deformed neutral axis is a space curve defined by position vector R*(s*,t*), which can be related to local tangent d3 as
R*0 ¼ d3
(2)
The evolution of the frame {d1,d2,d3} along the deformed rod is governed by the vector equation * d*0 i ¼ U di
ði ¼ 1; 2; 3Þ
(3)
in Eqs. (2)e(3) represents the derivative with respect to s . U* is the generalized strain vector (van der Heijden et al., 2003), ()0
U* ¼ k*1 d1 þ k*2 d2 þ t* d3
*
(4)
k*1 and k*2 are the curvatures of projections of the neutral axis on the d2 d3 and d1 d3 planes, respectively. t* is composed of the geometric torsion of the neutral axis and the physical twist of the rod. In the case when body fixed triad {d1,d2,d3} coincides with the Frenet trihedron of the neutral axis, then k*1 is zero and t* contains only the geometric torsion. The internal force F*(s*,t*) and internal moment M*(s*,t*) can be written as
F* ¼ F1* d1 þ F2* d2 þ F3* d3 ¼ Fx* ex þ Fy* ey þ Fz* ez
(5)
M* ¼ M1* d1 þ M2* d2 þ M3* d3 ¼ Mx* ex þ My* ey þ Mz* ez
(6)
F1* and F2* are shear forces and F3* is the axial force. M1* and M2* are bending moments in the directions of d1 and d2, respectively. M3* is the twisting moment along the direction of d3. The constitutive equations of the rod can be written as
M1* ¼ EI1 k*1 ; M2* ¼ EI2 k*2 ; M3* ¼ GDt*
(7)
D in the third expression in Eq. (7) depends on the shape of the cross section. If the cross section is elliptic with semi-major axis a 3 3 b and semi-minor axis b, then D ¼ apa 2 þb2 . If the cross section is rectangular with sides a and b, then D ¼ ka3b, where k can be calculated from (Reismann and Pawlik, 1980)
k¼
∞ 1 192 a X 1 ð2n þ 1Þpb 5 tanh 3 3p b n¼0 ð2n þ 1Þ5 2a
(8)
The twist t* is zero in planar deformation, but may be nonzero in spatial deformation. By extending the formulations in (Coleman et al., 1993) and (Goriely and Tabor, 1997) from a rod with circular cross section to noncircular one, we can write the dimensionless governing equations of the spatial elastica under a midpoint force as
R0 ðs; tÞ ¼ d3 ðs; tÞ
Fig. 1. (a) An elastica subject to a midpoint force. (b) Directors {d1,d2,d3}.
0 € tÞ F ðs; tÞ Q dðs 1=2Þey ¼ Rðs;
(9) (10)
86
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
nh i 0 € ðs; tÞ M ðs; tÞ þ d3 ðs; tÞ Fðs; tÞ ¼ a d1 ðs; tÞ d 1 io h € ðs; tÞ þ l d2 ðs; tÞ d 2 d01 ðs; tÞ ¼ Mðs; tÞ d1 ðs; tÞ þ
1g ½Mðs; tÞ d3 ðs; tÞ g
(11)
FðsÞ ¼ F0 þ QHðs 1=2Þey
(21)
H() in Eq. (21) is the Heaviside function. F0 is the internal force at s ¼ 0. The internal moment M(s) in Eqs. (12) and (14) can be obtained from Eq. (16) as
½d1 ðs; tÞ d2 ðs; tÞ (13)
(14)
½d1 ðs; tÞ d3 ðs; tÞ
d() in Eq. (10) is the Dirac delta function. The relations between the dimensionless variables (without asterisks) and the physical ones are *
1 * * * * * M L0 ; R ;s ;L ;a ;b ; M ¼ L0 EI2 L2 GD EI ðF; Q Þ ¼ 0 F* ; Q * ; U ¼ L0 U* ; g ¼ ; l ¼ 1; EI2 EI2 EI2 sffiffiffiffiffiffiffi sffiffiffiffiffiffiffi I mA * 1 EI2 * u ; t¼ 2 a ¼ 22 ; u ¼ L20 t ; dð Þ ¼ L0 d* ð Þ EI2 AL0 L0 mA ðR; s; L; a; bÞ ¼
(20)
For static analysis, the time-dependent variables in Eqs. (10) and (16) vanish. From Eq. (10) one obtains
1l ½Mðs; tÞ d1 ðs; tÞ ¼ Mðs; tÞ d3 ðs; tÞ þ l
M ¼ Mx ; My ; Mz 3. Static analysis
1g ½Mðs; tÞ d3 ðs; tÞ g 1l ½Mðs; tÞ d1 ðs; tÞ ½d3 ðs; tÞ d2 ðs; tÞ þ l
d03 ðs; tÞ
(19)
(12)
½d3 ðs; tÞ d1 ðs; tÞ d02 ðs; tÞ ¼ Mðs; tÞ d2 ðs; tÞ þ
F ¼ Fx ; Fy ; Fz
(15)
u is a dimensionless natural frequency. a is a slenderness parameter. g is a parameter dependent upon material properties and the shape of cross section. For a cross section with equal principal moments of inertia, the bending stiffness ratio l ¼ 1. One may rewrite Eq. (11) by replacing d2 with d3 d1,
MðsÞ ¼ RðsÞ F0 QHðs 1=2ÞRðsÞ ey þ M0
(22)
M0 is the moment at s ¼ 0. One can obtain the static solution by solving Eqs. (9), (12) and (14) for unknowns R(s), d1(s), and d3(s). The boundary conditions at s ¼ 0 can be written as
Rð0Þ ¼ ð0; 0; 0Þ
(23)
d3z ð0Þ ¼ 0
(24)
d1z ð0Þ ¼ 0
(25)
Mz ð0Þ ¼ 0
(26)
The boundary conditions at s ¼ 1 are
Rð1Þ ¼ ðL; 0; 0Þ
(27)
d3z ð1Þ ¼ 0
(28)
d1z ð1Þ ¼ 0
(29)
Mz ð1Þ ¼ 0
(30)
We denote the end angle at s ¼ 0 q0. We then have d3x0 ¼ cosq0, d3y0 ¼ sinq0, d1x0 ¼ sinq0, and d1y0 ¼ cosq0. A shooting method is used to solve for the static deformation of the rod. After making
nh i h € ðs; tÞ þ l € ðs; tÞ d ðs; tÞ d ðs; tÞ þ 2 M0 ðs; tÞ þ d3 ðs; tÞ Fðs; tÞ ¼ a d1 ðs; tÞ d d3 ðs; tÞ d d_ 3 ðs; tÞ d_ 1 ðs; tÞ d3 ðs; tÞ d1 ðs; tÞ 1 3 1 1 io € ðs; tÞ d ðs; tÞ d ðs; tÞ d_ ðs; tÞ d_ ðs; tÞ d ðs; tÞ d ðs; tÞ þ d ðs; tÞ d 3
1
1
3
1
1
3
3
(16)
Eqs. (9), (10), (12), (14), and (16) are five coupled nonlinear equations which can be used to solve for the five vector functions R, F, M, d1, and d3. For ease of writing, we expand all vectors in component forms within the {ex,ey,ez} frame as follows
R ¼ ðx; y; zÞ di ¼ dix ; diy ; diz
(17)
guesses on six parameters F0x, F0y, F0z, M0x, M0y, q0, together with conditions (23)e(26), we can integrate the nine scalar equations in Eqs. (9), (12) and (14) as an initial value problem from s ¼ 0 to 1. The six conditions (27)e(30) are used to check the accuracy of the six guesses. If the result is not satisfactory, a NewtoneRaphson algorithm is adopted for a new set of guesses. After finding R(s), d1(s), and d3(s), the other two functions F and M can be obtained readily.
4. Vibration analysis
ði ¼ 1; 2; 3Þ
(18)
The equilibrium configurations obtained from static analysis may be stable or unstable. In order to predict the behavior of the
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
loaded rod, the stability of these deformations must be determined. In this section we use vibration method to find the natural frequencies of the loaded rod, from which the stability of the elastica can be determined. The static solutions of the rod are denoted Re(s), Fe(s), Me(s), d1e(s), and d3e(s). The end angle at s ¼ 0 is denoted q0e. We superpose small harmonic perturbations onto these static solutions in the forms
Rðs; tÞ ¼ Re ðsÞ þ Rd ðsÞsin ut di ðs; tÞ ¼ die ðsÞ þ did ðsÞsin ut
(31) ði ¼ 1; 3Þ
87
d1y ð0; tÞ ¼ d1ye ð0Þ þ d1yd ð0Þsin ut ¼ cosðq0e þ q0d sin utÞ (46) d3x ð0; tÞ ¼ d3xe ð0Þ þ d3xd ð0Þsin ut ¼ cosðq0e þ q0d sin utÞ (47) d3y ð0; tÞ ¼ d3ye ð0Þ þ d3yd ð0Þsin ut ¼ sinðq0e þ q0d sin utÞ
(48)
After linearizing Eqs. (45)e(48), one obtains
(32)
d1xd ð0Þ ¼ q0d cos q0e
(49)
Fðs; tÞ ¼ Fe ðsÞ þ Fd ðsÞsin ut
(33)
d1yd ð0Þ ¼ q0d sin q0e
(50)
Mðs; tÞ ¼ Me ðsÞ þ Md ðsÞsin ut
(34)
d3xd ð0Þ ¼ q0d sin q0e
(51)
q0 ðtÞ ¼ q0e þ q0d sin ut
(35)
d3yd ð0Þ ¼ q0d cos q0e
(52)
After substituting Eqs. (31)e(34) into Eqs. (9), (10), (12), (14) and (15) and ignoring the higher order terms, one can obtain the linearized equations,
R0d ðsÞ ¼ d3d ðsÞ
(36)
F0d ðsÞ ¼ u2 R d ðsÞ
(37)
M0d ðsÞ ¼ ½Fe ðsÞ d3d ðsÞ þ Fd ðsÞ d3e ðsÞ þ au2 f½d1d ðsÞ d1e ðsÞ þ l½ð½d3e ðsÞ d1e ðsÞ d3d ðsÞÞd1e ðsÞ ð½d3e ðsÞ d1e ðsÞ d1d ðsÞÞd3e ðsÞg (38) d01d ðsÞ ¼ Me ðsÞ d1d ðsÞ þ Md ðsÞ d1e ðsÞ þ
1g f½Me ðsÞ d3e ðsÞ g
ðd3e ðsÞ d1d ðsÞ þ d3d ðsÞ d1e ðsÞÞ þ ½Me ðsÞ d3d ðsÞ þ Md ðsÞ d3e ðsÞðd3e ðsÞ d1e ðsÞÞg (39) d03d ðsÞ ¼ Me ðsÞ d3d ðsÞ þ Md ðsÞ d3e ðsÞ þ
1l f½Me ðsÞ d1e ðsÞ l
ðd1e ðsÞ d3d ðsÞ þ d1d ðsÞ d3e ðsÞÞ þ ½Me ðsÞ d1d ðsÞ þ Md ðsÞ d1e ðsÞðd1e ðsÞ d3e ðsÞÞg (40) Eqs. (36)e(40) are five vector equations which can be used to solve for the frequency and mode shape of the elastica. The boundary conditions at s ¼ 0 and 1 are
Rd ¼ ð0; 0; 0Þ
(41)
d1zd ¼ 0
(42)
d3zd ¼ 0
(43)
Mzd ¼ 0
(44)
It is noted that at the end s ¼ 0, d1x, d1y, d3x, and d3y are function of time during vibration, and can be written as,
d1x ð0; tÞ ¼ d1xe ð0Þ þ d1xd ð0Þsin ut ¼ sinðq0e þ q0d sin utÞ
(45)
The fifteen linearized scalar differential Eqs. (36)e(40) admit nontrivial solutions only when u is equal to an eigenvalue of the system of equations. In the solution method we first set Fxd0 ¼ 1. After guessing 6 variables Fyd0, Fzd0, Mxd0, Myd0, q0d, and u2, one can integrate the differential Eqs. (36)e(40) from s ¼ 0 to s ¼ 1. The assigned Fxd0 ¼ 1, and guessed Fyd0, Fzd0, Mxd0, Myd0, q0d, together with boundary conditions (41)e(44) at s ¼ 0 and (49)e(52) provide the needed initial conditions. The six conditions in (41) and (44) at s ¼ 1 are used to check the accuracy of the 6 guesses. 5. Numerical examples In the calculation we specify material property G/E ¼ 0.376 and
a ¼ 1.767 106. These parameters correspond to the metallic wire we use in the lab demonstration which will be described later. We first consider a case with circular cross section (l ¼ 1) and L ¼ 0.75. Fig. 2(a) shows the midpoint displacement in the negative y-direction, denoted Dh, as a function of midpoint force Q. The solid and dashed curves represent stable and unstable solutions, respectively. The thick and thin lines represent planar (2D) and spatial (3D) deformations, respectively. As Q increases from point O (Q ¼ 0), the elastica starts with planar symmetric deformation (1) until point B. At bifurcation point B, the planar deformation becomes unstable and there emerges a stable spatial symmetric deformation (2). As Q increases up to a limit point S1 (with a vertical tangent), the spatial deformation becomes unstable and the elastica jumps to another planar symmetric deformation (3). The deformation evolution from (1) to (3) described above was depicted in Fig. 2(b) in two viewing directions. Fig. 3 shows the first four u2 of planar deformation (1) as Q increases. The mode shapes corresponding to u1 and u3 are out-ofplane modes, while those corresponding to u2 and u4 are in-plane modes, as shown in Fig. 4 when Q ¼ 0. The solid and dashed lines represent equilibrium configurations and vibrating mode shapes, respectively. The u1 mode has no nodal point. The u2 mode has one nodal point, while u3 and u4 modes have two nodal points. From Fig. 3 we note that the u21 becomes negative at point B (Q ¼ 33.87). This implies that the planar deformation (1) becomes unstable at this bifurcation point. The stability of all other equilibrium positions is determined in the similar manner. In order to see if the above deformation sequence indeed occurs in reality, we build an experimental set-up to mimic the situation in the theory. The elastic rod used in the experiment is made of titaniumenickel alloy with 0.6 mm of diameter. The Young's modulus was measured as 81.3 GPa in an Insight 10 kN tensile test machine by MTS System. The mass density is 8256 kg/m3. The experimental
88
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
Fig. 4. The mode shapes of the first four modes in Fig. 3 when Q ¼ 0. The solid and dashed lines represent equilibrium configurations and vibrating mode shapes, respectively.
Fig. 2. (a) Relation between midpoint height change Dh and force Q whenl ¼ 1 and L ¼ 0.75. (b) Deformation evolution.
set-up is shown in Fig. 5. Both ends of the rod (1) are clamped in small aluminum blocks (2), which are then pinned about axes (3) in frame (5) via roller bearings. The distance of the two pin axes are adjustable by moving the bolts in frame (5) along a linear groove (6). A string (4) is tied onto the midpoint. The other end of the string is attached with a small bucket. Water is flowed into the bucket slowly to mimic the quasi-static point load. The experimental setup is not a very precise mimic of the theoretical model because of the finite size of the end clamps. However, it still can be used to observe the deformation sequence described in Fig. 2. Fig. 6 shows the photographs of deformations (1), (2), and (3) from the front and the side. In this case, the distance between the two ends L* is 20 cm.
Fig. 5. Photograph of the experimental set-up. (1) elastic rod; (2) clamps; (3) pin axis; (4) loading string; (5) adjustable frame; (6) linear groove; (7) fixed frame.
6. Other deformation patterns The deformation pattern described in the previous section is merely one of the many kinds. In this paper we assume that the Fig. 6. Front and side views of the deformation sequence observed in Type I.
Fig. 3. The relation between u2 and Q of the planar deformation (1).
cross section of the beam is of oval shape. If we change the stiffness ratio l and the distance L, we may have ten different deformation patterns when l 1, as shown in the phase diagram Fig. 7(a). The deformation pattern described in the previous section is Type I. Fig. 7(a) is repeated in Fig. 7(b) with different scale. In the following, we describe briefly the behavior of each of these patterns. The analysis procedure is similar to the one in the previous section. For cases when l ¼ 1, the deformation sequences are verified by experiment. Not all cases are accompanied with photographs of the experiment in this paper. The interested reader may refer to the Master Thesis of the second author (Chen, 2014) for more details.
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
89
6.1. Type II In Fig. 8 we pull the two ends in Fig. 2 further away to L ¼ 0.85. In doing so the limit point on the path of 3D deformation (2) ceases to exist. As a consequence the planar symmetric deformation (1) jumps to another planar symmetric deformation (3) at bifurcation point B. The jump is accomplished via a 3D symmetric deformation shown as dashed curve in Fig. 8(b). In the experiment we record the dynamic jump by video camera and play it back slowly. The 3D jump can indeed be observed.
6.2. Type III In Fig. 9 we push the two ends in Fig. 2 closer to L ¼ 0.6. In this case there appear two spatial deformations instead of one. The planar deformation (1) evolves to a stable 3D symmetric deformation (2). At point S1, there emerges another 3D deformation (3) which is asymmetric with respect to the midpoint. After S1 both spatial deformations become unstable and the elastica jumps to planar symmetric deformation (4). Point S1 is a sub-critical bifurcation point. In Fig. 9(b) deformation (2) jumps to (4) via 3D asymmetric deformation (3), as shown by the dashed curve.
Fig. 7. (a) Deformation types in the lL plane; (b) in different scale.
Fig. 8. (a) Load-deflection diagram and (b) deformation evolution of type II (l ¼ 1, L ¼ 0.85).
Fig. 9. (a) Load-deflection diagram and (b) deformation evolution of type III (l ¼ 1, L ¼ 0.6).
90
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
Fig. 11. Two stable deformations (4) and (5) of type IV observed in lab.
Fig. 10. (a) Load-deflection diagram and (b) deformation evolution of type IV (l ¼ 1, L ¼ 0.45).
6.3. Type IV In Fig. 10 we push the two ends even closer to L ¼ 0.45. The planar symmetric deformation (1) evolves to a stable 3D symmetric deformation (2) at B. At point S there emerges a 3D asymmetric deformation (3). The elastica jumps at S, similar to the case in Type III. However, in Type IV there exists one additional 3D deformation (5) which is anti-symmetric with respect to the midpoint and twisted. Deformation (5) becomes stable after it merges with deformation (3) at Point C. As a consequence, there are two stable positions available after the jump, i.e., deformations (4) and (5). It is not an easy matter to predict which way to go after the jump without analyzing the potential energy structure, which is beyond the scope of this paper. In the experiment we found that the elastica jumps at point S from 3D symmetric deformation (2) to deformation (5) via 3D asymmetric deformation (3). Fig. 10(b) shows the deformation sequence. Although we are unable to reach deformation (4) by increasing Q quasi-statically, we can produce this position by hand. The photographs in Fig. 11 show that these two positions indeed exist in reality. In other words, they are both stable. It is noted that on the path of deformation (3) there exists a small stable segment between A1 and A2. Both A1 and A2 are limit
Fig. 12. (a) Load-deflection diagram and (b) deformation evolution of type V (l ¼ 1, L ¼ 0.3).
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
91
Fig. 14. (a) Load-deflection diagram and (b) deformation evolution of type VII (l ¼ 1, L ¼ 0.2).
Fig. 13. (a) Load-deflection diagram and (b) deformation evolution of type VI (l ¼ 1, L ¼ 0.25).
points. Because the range of Q is so narrow that we are unable to produce it in the experiment. 6.4. Type V Fig. 12 shows the case when we push the two ends further closer to L ¼ 0.3. At point B the planar symmetric deformation (1) bifurcates to a stable planar asymmetric deformation (2). At point S there emerges a 3D asymmetric deformation (3). There exists a small stable segment between A1 and A2 on the path of deformation (3), similar to Type IV in Fig. 10. Deformation (3) merges to an independent path of 3D deformation (5) which is anti-symmetric and twisted. Deformation (5) is stable beyond point C. The elastica jumps at point B from planar asymmetric deformation (2) to 3D twisted deformation (5) via 3D asymmetric deformation (3). 6.5. Type VI Fig. 13 shows the case when L ¼ 0.25. There is no spatial deformation in this case. The path of deformation (2) has a limit point S, as shown in Fig. 13(a). Deformation (2) is stable between B and S. Therefore, the planar symmetric deformation (1) bifurcates to planar asymmetric deformation (2) and then jumps at point S to symmetric deformation (3), as shown in Fig. 13(b). The elastica
crosses itself at the midpoint in deformation (3). It is noted that during snapping, the elastica moves sideway in the plane and crosses the pin joint. Due to the finite size of the pin joint in the experiment, the elastica will be blocked by the pin joint following the jump. 6.6. Type VII Fig. 14 shows the case when L ¼ 0.2. At point B the planar symmetric deformation (1) bifurcates to a planar asymmetric deformation (2). Both (1) and (2) are unstable beyond B. Deformation (2) merges at C to a planar symmetric deformation (3), which is stable beyond C. The elastica jumps at point B from deformation (1) to (3) via asymmetric deformation (2), as shown in Fig. 14(b). 6.7. Type VIII Up to this point, we limit our discussion to the case when l ¼ 1. If we set l equal to 2 and start from small L, the elastica behaves from Type VII to VI, and V. If we increase L to 0.6, we encounter a new type of deformation behavior, Type VIII, as shown in Fig. 15. If we compare the load-deflection diagram with the one in Fig. 9 (type III, for l ¼ 1 and L ¼ 0.6), we found that they are similar in structure, except that the 3D symmetric deformation (2) in Type III is replaced by a 2D asymmetric deformation (2) in Type VIII. Apparently, the elastica behaves more like a planar elastica as l increases. As shown in Fig. 15, deformation (3) starts at S1 and ends at S2 on the path of planar asymmetric deformation (2). The deformation sequence is depicted in Fig. 15(b), in which planar deformation (1) bifurcates to
92
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
planar asymmetric deformation (2) and then jumps to planar symmetric deformation (4) via a 3D asymmetric deformation (3). 6.8. Type IX If we increase L from 0.6 to 0.7, the spatial deformation disappears altogether. In Fig. 16 the planar symmetric deformation (1) bifurcates to planar asymmetric deformation (2) and then jumps to planar symmetric deformation (3). The difference between Types VIII and IX is that the elastica in Type VIII jumps via a spatial deformation while the elastica in Type IX jumps via a planar asymmetric deformation. 6.9. Type X For the case when L ¼ 0.9, the bifurcation point B occurs before the limit point, as shown in Fig. 17(a). Therefore, the planar symmetric deformation (1) jumps directly to planar symmetric deformation (3) via a planar asymmetric deformation (2), as depicted in Fig. 17(b). It is noted that the behaviors in Types IX and X are purely planar and have been studied previously by Chen and Hung (2011) with planar elastica model. 7. Conclusions
Fig. 15. (a) Load-deflection diagram and (b) deformation evolution of type VIII (l ¼ 2, L ¼ 0.6).
Fig. 16. (a) Load-deflection diagram and (b) deformation evolution of type IX (l ¼ 2, L ¼ 0.7).
In this paper we study the behavior of a spatial elastica under a midpoint force. The rod under question is buckled into a curved shape and then the two ends are fixed in space with pin joints. Both static and dynamic analyses are conducted. Focus is placed on the snapping phenomenon. The behavior of the elastica depends on two parameters; the distance between the two ends L and the bending stiffness ratio of the cross section l. A total of ten deformation types are identified. Their behaviors are summarized in Table 1, in which we use four characteristics to categorize the deformation patterns: (1) the deformation before jump, (2) the deformation after jump, (3) the type of critical point when jump
Fig. 17. (a) Load-deflection diagram and (b) deformation evolution of type X (l ¼ 2, L ¼ 0.9).
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
93
Table 1 Classification of deformation behaviors. Type
Deformation before jump
Deformation after jump
Critical point
Jump via
I II III IV V VI VII VIII IX X
3D 2D 3D 3D 2D 2D 2D 2D 2D 2D
2D 2D 2D 3D 3D 2D 2D 2D 2D 2D
Limit Bifurcation Bifurcation Bifurcation Bifurcation Limit Bifurcation Bifurcation Limit Bifurcation
3D 3D 3D 3D 3D 2D 2D 3D 2D 2D
symmetric symmetric symmetric symmetric asymmetric asymmetric symmetric asymmetric asymmetric symmetric
symmetric symmetric symmetric twisted twisted symmetric (self-crossing) symmetric (self-crossing) symmetric symmetric symmetric
occurs, and (4) the deformation during the jump. It can be seen that no two deformation patterns are the same in all four characteristics. It is noted that types VI, VII, IX, and X exhibits only planar deformations. Therefore, these deformations can be obtained by using planar elastica model. By observing the phase diagram in Fig. 7 we can see that for the case when l is greater than 28, only these four types of deformation can exist in all ranges of L. We therefore conclude that if one wishes to design an elastica capable of only plane motion in all ranges of L, then l should be greater than 28. It is noted that a shallow planar buckled beam may buckle from a symmetric deformation to another symmetric deformation symmetrically (Chen and Hung, 2012). This occurs when the buckled beam is extensible and flat (L ¼ 1 5.6p2a). With the a in this article, this occurs when L > 0.9999. This type of snapping behavior will not be found in Fig. 6 because the rod in this paper is assumed to be inextensible here. The method presented in this paper is meant for a pinnedepinned elastica. However, the similar method can be applied to other types of boundary conditions, such as clampedeclamped spatial elastica. It is noted that the above conclusions are based on the numerical results on a specific specimen with G/E ¼ 0.376 and slenderness parameter a ¼ 1.767 106. In addition, the cross section is assumed to be of oval shape. If any of these are changed, some of the conclusions may be modified. For instance, if the cross section is changed to rectangular shape, the phase diagram look similar topologically, but the threshold l above which only plane deformation is possible will be changed to 37. From Eq. (11) we see that the slenderness parameter a affects only the inertial terms. Therefore, it does not affect the static equilibrium. In practical applications, the G/E ratios of most metals are in the range between 0.3 and 0.4. Although we did not test all the possible G/E ratios, it is expected that the overall structure of the phase diagram should remain the same in this range. References Cazottes, P., Fernandes, A., Pouget, J., Hafez, M., 2009. Bistable buckled beam: modeling of actuating force and experimental validations. ASME J. Mech. Des. 131, 101001. Chen, J.-S., Ro, W.-C., Lin, J.-S., 2009. Exact static and dynamic critical loads of a shallow arch under a point force at the midpoint. Int. J. Nonlin. Mech. 44, 66e70.
symmetric symmetric asymmetric asymmetric asymmetric asymmetric asymmetric asymmetric asymmetric asymmetric
Chen, J.-S., Hung, S.Y., 2011. Snapping of an elastica under various loading mechanisms. Eur. J. Mech. A Solid 30, 525e531. Chen, J.-S., Hung, S.-Y., 2012. Exact snapping loads of a buckled beam under a midpoint force. Appl. Math. Modell. 36, 1776e1782. Chen, J.-S., Tsao, H.-W., 2013. Static snapping load of a hinged extensible elastica. Appl. Math. Modell. 37, 8401e8408. Chen, R.-M., 2014. Deformation and Vibration of a Spatial Elastic under a Midpoint Force. Master Thesis. Department of Mechanical Engineering, National Taiwan University. in Chinese. Coleman, B.D., Dill, E.H., Lembo, M., Lu, Z., Tobias, I., 1993. On the dynamics of rods in the theory of Kirchhoff and Clebsch. Arch. Ration. Mech. Anal. 121, 339e359. Fung, Y.C., Kaplan, A., 1952. Buckling of Low Arches or Curved Beams of Small Curvature. NACA Technical Note 2840. Goriely, A., Tabor, M., 1997. Nonlinear dynamics of filaments, I. Dynamical instability. Phys. D 105, 20e44. Krylov, S., Ilic, B.R., Schreiber, D., Seretensky, S., Craighead, H., 2008. The pull-in behavior of electrostatically actuated bistable microstructures. J. Michromech. Microeng. 18, 055026. Kublanov, L.B., Bottega, W.J., 1995. On pressing a buckled film. Appl. Math. Modell. 19, 499e507. Maurini, C., Pouget, J., Vidoli, S., 2007. Distributed piezoelectric actuation of a bistable buckled beam. Eur. J. Mech. A Solid 26, 837e853. Medina, L., Gilat, R., Krylov, S., 2012. Symmetry breaking in an initially curved micro beam loaded by a distributed electrostatic force. Int. J. Solids Struct. 49, 1864e1876. Ouakad, H.M., Younis, M.I., 2014. On using the dynamic snap-through motion of MEMS initially curved microbeams for filtering applications. J. Sound Vib. 333, 555e568. Patricio, P., Adda-Bedia, M., Ben Amar, M., 1998. An elastica problem: instabilities of an elastic arch. Phys. D 124, 285e295. Pinto, O.C., Goncalves, P.B., 2000. Nonlinear control of buckled beams under step loading. Mech. Syst. Signal Process. 14, 967e985. Pippard, A.B., 1990. The elastic arch and its modes of instability. Eur. J. Phys. 11, 359e365. Plaut, R.H., 1978. Influence of load position on the stability of shallow arches. J. Appl. Math. Phys. (ZAMP) 30, 548e552. Reismann, H., Pawlik, P.S., 1980. Elasticity: Theory and Applications. John Wiley & Sons, New York. Schreyer, H.L., 1972. The effect of imperfections on the buckling load of shallow circular arches. ASME J. Appl. Mech. 445e450. Seide, P., 1984. Dynamic stability of laterally loaded buckled beams. J. Eng. Mech. 110, 1556e1569. Vangbo, M., 1998. An analytical analysis of a compressed bistable buckled beam. Sens. Actuat. A 69, 212e216. van der Heijden, G.H.M., Neukirch, S., Goss, V.G.A., Champneys, Thompson, J.M.T., 2003. Instability and self-contact phenomena in the writhing of clamped rods. Int. J. Mech. Sci. 45, 161e196. Virgin, L.N., Wiebe, R., Spottswood, S.M., Eason, T.G., 2014. Sensitivity in the structural behavior of shallow arches. Int. J. Nonlin. Mech. 58, 212e221. Zang, Y., Wang, Y., Li, Z., Huang, Y., Li, D., 2007. Snap-through and pull-in instabilities of arch-shaped beam under an electrostatic loading. J. Mechromech. Syst. 16, 684e693.
Contents lists available at ScienceDirect
European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol
Deformation and stability of a spatial elastica under a midpoint force Jen-San Chen*, Re-Ming Chen Department of Mechanical Engineering, National Taiwan University, Taipei, 10617, Taiwan
a r t i c l e i n f o
a b s t r a c t
Article history: Received 23 December 2014 Accepted 7 June 2015 Available online 20 June 2015
In this paper we study the deformation and stability of a pinnedepinned buckled beam under the action of a concentrated force at the midpoint. Focus is placed on the snap-through phenomenon, which may take place in a plane or three-dimensional space. We first find the equilibrium configurations by using shooting method. Elastica model is adopted to take into account exact geometry in large deformation. As expected, multiple solutions may exist for a specified set of loading parameters. Vibration method is then employed to determine the stability of the equilibrium solutions. Through these analyses the deformation sequence as the midpoint force increases quasi-statically can be predicted. It is found that the deformation sequence of the elastica is determined by two parameters; (1) the distance between the two ends of the buckled beam, and (2) the bending stiffness ratio of the cross section. Ten different deformation patterns can be identified according to four characteristics; the deformations before, after, and during the jump, and the type of critical point at the jump. A metallic wire with circular cross section is used to verify the predicted deformation sequence. It is concluded that for the specific specimen in the demonstration if one wishes to design an elastica capable of only plane deformation in all range of end distance, then the bending stiffness ratio has to be greater than 28. © 2015 Elsevier Masson SAS. All rights reserved.
Keywords: Spatial elastica Midpoint force Snap-through
1. Introduction In the design of miniature mechanisms, such as in MEMS, it is impractical to use conventional rigid-body joint pairs due to space constraint. More often, the motion transmission is accomplished by the deformation of a single flexible beam. In particular, the snapthrough phenomenon of a curved beam has been used for microswitches (Maurini et al., 2007; Zang et al., 2007; Krylov et al., 2008; Medina et al., 2012; Ouakad and Younis, 2014). Two kinds of curved beams under midpoint forces can be found in the literature. In the first type, the curved beam is stress free when it is in curved configuration, see Fung and Kaplan (1952), Schreyer (1972), Plaut (1978), Chen et al. (2009), Virgin et al. (2014). In the second type the beam, which is stress free when it is straight, is buckled into a curved shape by edge thrust. The buckled beam is then loaded by a midpoint force, see Seide (1984), Pippard (1990), Patricio et al. (1998), Vangbo (1998), Kublanov and Bottega (1995), Pinto and Goncalves (2000), Cazottes et al. (2009), Chen and Hung (2011, 2012), Chen and Tsao (2013).
* Corresponding author. E-mail address: [email protected] (J.-S. Chen). http://dx.doi.org/10.1016/j.euromechsol.2015.06.006 0997-7538/© 2015 Elsevier Masson SAS. All rights reserved.
The snap-through phenomena observed in these studies are limited in plane deformation. Pippard (1990) noted in his experiment that a sway mode of instability normal to the plane of the strip tends to occur. Intuitively, if the ratio between the width and the thickness of the strip is sufficiently large, the strip should behave like a planar elastica. This raises the question what the minimum ratio is in order for the strip to deform only in a plane. Furthermore, what will happen if three-dimensional deformation does occur? In this paper we pursue this interesting phenomenon by considering an elastic rod capable of out-of-plane deformation. Elastica model is adopted in the analysis. The elastic rod considered in this paper is stress free when it is straight. The two ends of the rod are pinned in space after it is buckled into a curved shape. A force is applied at the midpoint in the plane of the rod. The two principal moments of inertia of the cross section may be different. If the bending stiffness in one of the principal directions is much larger than the other, it is expected that the deformation will be restricted in a plane, as described in the works cited previously. If the two bending stiffness are comparable, on the other hand, spatial deformation may occur. It is the objective of this paper to identify all these deformation patterns. In Section 2 we establish the equations of motion of the loaded rod. In Section 3 we conduct a static analysis to find the equilibrium configurations. In Section 4 we study the vibration characteristics
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
of the loaded elastica. The stability of the loaded elastica can then be determined from the calculated natural frequencies. In Section 5 we present in detail a numerical example showing the loaddeflection diagram and the frequency spectrum. A pinned metallic wire with circular cross section is used to verify the predicted deformation sequence as the midpoint load increases quasistatically. In Section 6 we show a phase diagram using the distance between the two ends and the stiffness ratio between the two principal moments as two parameters. In the diagram, ten different deformation behaviors are identified. In Section 7 several conclusions are summarized.
2. Equations of motion We consider a uniform, inextensible, and unshearable elastic rod with length L0, cross section area A, Young's modulus E, shear modulus G, and mass density per unit volume m. The area moment of inertia in the two principal directions of the cross section are I1 and I2, where I1 I2. The rod is initially straight and stress free. The rod is first buckled in the plane containing the principal direction with I2 by pushing the two ends closer to a distance L*. The two ends are then fixed in space with pin joints, as shown in Fig. 1(a). A space-fixed x*y*z*-coordinate system with origin attached to the left end O is chosen to describe the geometry of the rod. The x*-axis is pointing to the fixed end on the right. The inertial orthonormal frame {ex,ey,ez} is associated with the x*y*z*-coordinate system. The location of a material point on the neutral axis of the deformed rod is denoted by arc length s* measured from end O. We consider the case when the buckled beam is loaded by point force Q* at the midpoint s* ¼ L0/2 in the y* direction. A body-fixed right-handed orthonormal frame {d1,d2,d3} is chosen in such a manner that vector d3 is in the direction of the local tangent of the deformed neutral axis. Vectors d1 and d2 are in the normal cross section of the rod and rotate along with the cross section. Fig. 1(b) shows the case when the cross section is of elliptic
85
shape with semi-major axis b and semi-minor axis a. When the rod is in the unstressed straight state, the frame {d1,d2,d3} coincides with the frame {ex,ey,ez} in such a manner that d3 ≡ ex, d1 ≡ ey, d2 ≡ ez. The vectors di(s*,t*) (i ¼ 1, 2, 3), where t* is time, can be expanded as
di ¼ dix ex þ diy ey þ diz ez
ði ¼ 1; 2; 3Þ
(1)
It is assumed that the cross section of the rod remains plane and normal to the neutral axis after deformation. The rotation of the cross section can be defined by the vector d1. The deformed neutral axis is a space curve defined by position vector R*(s*,t*), which can be related to local tangent d3 as
R*0 ¼ d3
(2)
The evolution of the frame {d1,d2,d3} along the deformed rod is governed by the vector equation * d*0 i ¼ U di
ði ¼ 1; 2; 3Þ
(3)
in Eqs. (2)e(3) represents the derivative with respect to s . U* is the generalized strain vector (van der Heijden et al., 2003), ()0
U* ¼ k*1 d1 þ k*2 d2 þ t* d3
*
(4)
k*1 and k*2 are the curvatures of projections of the neutral axis on the d2 d3 and d1 d3 planes, respectively. t* is composed of the geometric torsion of the neutral axis and the physical twist of the rod. In the case when body fixed triad {d1,d2,d3} coincides with the Frenet trihedron of the neutral axis, then k*1 is zero and t* contains only the geometric torsion. The internal force F*(s*,t*) and internal moment M*(s*,t*) can be written as
F* ¼ F1* d1 þ F2* d2 þ F3* d3 ¼ Fx* ex þ Fy* ey þ Fz* ez
(5)
M* ¼ M1* d1 þ M2* d2 þ M3* d3 ¼ Mx* ex þ My* ey þ Mz* ez
(6)
F1* and F2* are shear forces and F3* is the axial force. M1* and M2* are bending moments in the directions of d1 and d2, respectively. M3* is the twisting moment along the direction of d3. The constitutive equations of the rod can be written as
M1* ¼ EI1 k*1 ; M2* ¼ EI2 k*2 ; M3* ¼ GDt*
(7)
D in the third expression in Eq. (7) depends on the shape of the cross section. If the cross section is elliptic with semi-major axis a 3 3 b and semi-minor axis b, then D ¼ apa 2 þb2 . If the cross section is rectangular with sides a and b, then D ¼ ka3b, where k can be calculated from (Reismann and Pawlik, 1980)
k¼
∞ 1 192 a X 1 ð2n þ 1Þpb 5 tanh 3 3p b n¼0 ð2n þ 1Þ5 2a
(8)
The twist t* is zero in planar deformation, but may be nonzero in spatial deformation. By extending the formulations in (Coleman et al., 1993) and (Goriely and Tabor, 1997) from a rod with circular cross section to noncircular one, we can write the dimensionless governing equations of the spatial elastica under a midpoint force as
R0 ðs; tÞ ¼ d3 ðs; tÞ
Fig. 1. (a) An elastica subject to a midpoint force. (b) Directors {d1,d2,d3}.
0 € tÞ F ðs; tÞ Q dðs 1=2Þey ¼ Rðs;
(9) (10)
86
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
nh i 0 € ðs; tÞ M ðs; tÞ þ d3 ðs; tÞ Fðs; tÞ ¼ a d1 ðs; tÞ d 1 io h € ðs; tÞ þ l d2 ðs; tÞ d 2 d01 ðs; tÞ ¼ Mðs; tÞ d1 ðs; tÞ þ
1g ½Mðs; tÞ d3 ðs; tÞ g
(11)
FðsÞ ¼ F0 þ QHðs 1=2Þey
(21)
H() in Eq. (21) is the Heaviside function. F0 is the internal force at s ¼ 0. The internal moment M(s) in Eqs. (12) and (14) can be obtained from Eq. (16) as
½d1 ðs; tÞ d2 ðs; tÞ (13)
(14)
½d1 ðs; tÞ d3 ðs; tÞ
d() in Eq. (10) is the Dirac delta function. The relations between the dimensionless variables (without asterisks) and the physical ones are *
1 * * * * * M L0 ; R ;s ;L ;a ;b ; M ¼ L0 EI2 L2 GD EI ðF; Q Þ ¼ 0 F* ; Q * ; U ¼ L0 U* ; g ¼ ; l ¼ 1; EI2 EI2 EI2 sffiffiffiffiffiffiffi sffiffiffiffiffiffiffi I mA * 1 EI2 * u ; t¼ 2 a ¼ 22 ; u ¼ L20 t ; dð Þ ¼ L0 d* ð Þ EI2 AL0 L0 mA ðR; s; L; a; bÞ ¼
(20)
For static analysis, the time-dependent variables in Eqs. (10) and (16) vanish. From Eq. (10) one obtains
1l ½Mðs; tÞ d1 ðs; tÞ ¼ Mðs; tÞ d3 ðs; tÞ þ l
M ¼ Mx ; My ; Mz 3. Static analysis
1g ½Mðs; tÞ d3 ðs; tÞ g 1l ½Mðs; tÞ d1 ðs; tÞ ½d3 ðs; tÞ d2 ðs; tÞ þ l
d03 ðs; tÞ
(19)
(12)
½d3 ðs; tÞ d1 ðs; tÞ d02 ðs; tÞ ¼ Mðs; tÞ d2 ðs; tÞ þ
F ¼ Fx ; Fy ; Fz
(15)
u is a dimensionless natural frequency. a is a slenderness parameter. g is a parameter dependent upon material properties and the shape of cross section. For a cross section with equal principal moments of inertia, the bending stiffness ratio l ¼ 1. One may rewrite Eq. (11) by replacing d2 with d3 d1,
MðsÞ ¼ RðsÞ F0 QHðs 1=2ÞRðsÞ ey þ M0
(22)
M0 is the moment at s ¼ 0. One can obtain the static solution by solving Eqs. (9), (12) and (14) for unknowns R(s), d1(s), and d3(s). The boundary conditions at s ¼ 0 can be written as
Rð0Þ ¼ ð0; 0; 0Þ
(23)
d3z ð0Þ ¼ 0
(24)
d1z ð0Þ ¼ 0
(25)
Mz ð0Þ ¼ 0
(26)
The boundary conditions at s ¼ 1 are
Rð1Þ ¼ ðL; 0; 0Þ
(27)
d3z ð1Þ ¼ 0
(28)
d1z ð1Þ ¼ 0
(29)
Mz ð1Þ ¼ 0
(30)
We denote the end angle at s ¼ 0 q0. We then have d3x0 ¼ cosq0, d3y0 ¼ sinq0, d1x0 ¼ sinq0, and d1y0 ¼ cosq0. A shooting method is used to solve for the static deformation of the rod. After making
nh i h € ðs; tÞ þ l € ðs; tÞ d ðs; tÞ d ðs; tÞ þ 2 M0 ðs; tÞ þ d3 ðs; tÞ Fðs; tÞ ¼ a d1 ðs; tÞ d d3 ðs; tÞ d d_ 3 ðs; tÞ d_ 1 ðs; tÞ d3 ðs; tÞ d1 ðs; tÞ 1 3 1 1 io € ðs; tÞ d ðs; tÞ d ðs; tÞ d_ ðs; tÞ d_ ðs; tÞ d ðs; tÞ d ðs; tÞ þ d ðs; tÞ d 3
1
1
3
1
1
3
3
(16)
Eqs. (9), (10), (12), (14), and (16) are five coupled nonlinear equations which can be used to solve for the five vector functions R, F, M, d1, and d3. For ease of writing, we expand all vectors in component forms within the {ex,ey,ez} frame as follows
R ¼ ðx; y; zÞ di ¼ dix ; diy ; diz
(17)
guesses on six parameters F0x, F0y, F0z, M0x, M0y, q0, together with conditions (23)e(26), we can integrate the nine scalar equations in Eqs. (9), (12) and (14) as an initial value problem from s ¼ 0 to 1. The six conditions (27)e(30) are used to check the accuracy of the six guesses. If the result is not satisfactory, a NewtoneRaphson algorithm is adopted for a new set of guesses. After finding R(s), d1(s), and d3(s), the other two functions F and M can be obtained readily.
4. Vibration analysis
ði ¼ 1; 2; 3Þ
(18)
The equilibrium configurations obtained from static analysis may be stable or unstable. In order to predict the behavior of the
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
loaded rod, the stability of these deformations must be determined. In this section we use vibration method to find the natural frequencies of the loaded rod, from which the stability of the elastica can be determined. The static solutions of the rod are denoted Re(s), Fe(s), Me(s), d1e(s), and d3e(s). The end angle at s ¼ 0 is denoted q0e. We superpose small harmonic perturbations onto these static solutions in the forms
Rðs; tÞ ¼ Re ðsÞ þ Rd ðsÞsin ut di ðs; tÞ ¼ die ðsÞ þ did ðsÞsin ut
(31) ði ¼ 1; 3Þ
87
d1y ð0; tÞ ¼ d1ye ð0Þ þ d1yd ð0Þsin ut ¼ cosðq0e þ q0d sin utÞ (46) d3x ð0; tÞ ¼ d3xe ð0Þ þ d3xd ð0Þsin ut ¼ cosðq0e þ q0d sin utÞ (47) d3y ð0; tÞ ¼ d3ye ð0Þ þ d3yd ð0Þsin ut ¼ sinðq0e þ q0d sin utÞ
(48)
After linearizing Eqs. (45)e(48), one obtains
(32)
d1xd ð0Þ ¼ q0d cos q0e
(49)
Fðs; tÞ ¼ Fe ðsÞ þ Fd ðsÞsin ut
(33)
d1yd ð0Þ ¼ q0d sin q0e
(50)
Mðs; tÞ ¼ Me ðsÞ þ Md ðsÞsin ut
(34)
d3xd ð0Þ ¼ q0d sin q0e
(51)
q0 ðtÞ ¼ q0e þ q0d sin ut
(35)
d3yd ð0Þ ¼ q0d cos q0e
(52)
After substituting Eqs. (31)e(34) into Eqs. (9), (10), (12), (14) and (15) and ignoring the higher order terms, one can obtain the linearized equations,
R0d ðsÞ ¼ d3d ðsÞ
(36)
F0d ðsÞ ¼ u2 R d ðsÞ
(37)
M0d ðsÞ ¼ ½Fe ðsÞ d3d ðsÞ þ Fd ðsÞ d3e ðsÞ þ au2 f½d1d ðsÞ d1e ðsÞ þ l½ð½d3e ðsÞ d1e ðsÞ d3d ðsÞÞd1e ðsÞ ð½d3e ðsÞ d1e ðsÞ d1d ðsÞÞd3e ðsÞg (38) d01d ðsÞ ¼ Me ðsÞ d1d ðsÞ þ Md ðsÞ d1e ðsÞ þ
1g f½Me ðsÞ d3e ðsÞ g
ðd3e ðsÞ d1d ðsÞ þ d3d ðsÞ d1e ðsÞÞ þ ½Me ðsÞ d3d ðsÞ þ Md ðsÞ d3e ðsÞðd3e ðsÞ d1e ðsÞÞg (39) d03d ðsÞ ¼ Me ðsÞ d3d ðsÞ þ Md ðsÞ d3e ðsÞ þ
1l f½Me ðsÞ d1e ðsÞ l
ðd1e ðsÞ d3d ðsÞ þ d1d ðsÞ d3e ðsÞÞ þ ½Me ðsÞ d1d ðsÞ þ Md ðsÞ d1e ðsÞðd1e ðsÞ d3e ðsÞÞg (40) Eqs. (36)e(40) are five vector equations which can be used to solve for the frequency and mode shape of the elastica. The boundary conditions at s ¼ 0 and 1 are
Rd ¼ ð0; 0; 0Þ
(41)
d1zd ¼ 0
(42)
d3zd ¼ 0
(43)
Mzd ¼ 0
(44)
It is noted that at the end s ¼ 0, d1x, d1y, d3x, and d3y are function of time during vibration, and can be written as,
d1x ð0; tÞ ¼ d1xe ð0Þ þ d1xd ð0Þsin ut ¼ sinðq0e þ q0d sin utÞ
(45)
The fifteen linearized scalar differential Eqs. (36)e(40) admit nontrivial solutions only when u is equal to an eigenvalue of the system of equations. In the solution method we first set Fxd0 ¼ 1. After guessing 6 variables Fyd0, Fzd0, Mxd0, Myd0, q0d, and u2, one can integrate the differential Eqs. (36)e(40) from s ¼ 0 to s ¼ 1. The assigned Fxd0 ¼ 1, and guessed Fyd0, Fzd0, Mxd0, Myd0, q0d, together with boundary conditions (41)e(44) at s ¼ 0 and (49)e(52) provide the needed initial conditions. The six conditions in (41) and (44) at s ¼ 1 are used to check the accuracy of the 6 guesses. 5. Numerical examples In the calculation we specify material property G/E ¼ 0.376 and
a ¼ 1.767 106. These parameters correspond to the metallic wire we use in the lab demonstration which will be described later. We first consider a case with circular cross section (l ¼ 1) and L ¼ 0.75. Fig. 2(a) shows the midpoint displacement in the negative y-direction, denoted Dh, as a function of midpoint force Q. The solid and dashed curves represent stable and unstable solutions, respectively. The thick and thin lines represent planar (2D) and spatial (3D) deformations, respectively. As Q increases from point O (Q ¼ 0), the elastica starts with planar symmetric deformation (1) until point B. At bifurcation point B, the planar deformation becomes unstable and there emerges a stable spatial symmetric deformation (2). As Q increases up to a limit point S1 (with a vertical tangent), the spatial deformation becomes unstable and the elastica jumps to another planar symmetric deformation (3). The deformation evolution from (1) to (3) described above was depicted in Fig. 2(b) in two viewing directions. Fig. 3 shows the first four u2 of planar deformation (1) as Q increases. The mode shapes corresponding to u1 and u3 are out-ofplane modes, while those corresponding to u2 and u4 are in-plane modes, as shown in Fig. 4 when Q ¼ 0. The solid and dashed lines represent equilibrium configurations and vibrating mode shapes, respectively. The u1 mode has no nodal point. The u2 mode has one nodal point, while u3 and u4 modes have two nodal points. From Fig. 3 we note that the u21 becomes negative at point B (Q ¼ 33.87). This implies that the planar deformation (1) becomes unstable at this bifurcation point. The stability of all other equilibrium positions is determined in the similar manner. In order to see if the above deformation sequence indeed occurs in reality, we build an experimental set-up to mimic the situation in the theory. The elastic rod used in the experiment is made of titaniumenickel alloy with 0.6 mm of diameter. The Young's modulus was measured as 81.3 GPa in an Insight 10 kN tensile test machine by MTS System. The mass density is 8256 kg/m3. The experimental
88
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
Fig. 4. The mode shapes of the first four modes in Fig. 3 when Q ¼ 0. The solid and dashed lines represent equilibrium configurations and vibrating mode shapes, respectively.
Fig. 2. (a) Relation between midpoint height change Dh and force Q whenl ¼ 1 and L ¼ 0.75. (b) Deformation evolution.
set-up is shown in Fig. 5. Both ends of the rod (1) are clamped in small aluminum blocks (2), which are then pinned about axes (3) in frame (5) via roller bearings. The distance of the two pin axes are adjustable by moving the bolts in frame (5) along a linear groove (6). A string (4) is tied onto the midpoint. The other end of the string is attached with a small bucket. Water is flowed into the bucket slowly to mimic the quasi-static point load. The experimental setup is not a very precise mimic of the theoretical model because of the finite size of the end clamps. However, it still can be used to observe the deformation sequence described in Fig. 2. Fig. 6 shows the photographs of deformations (1), (2), and (3) from the front and the side. In this case, the distance between the two ends L* is 20 cm.
Fig. 5. Photograph of the experimental set-up. (1) elastic rod; (2) clamps; (3) pin axis; (4) loading string; (5) adjustable frame; (6) linear groove; (7) fixed frame.
6. Other deformation patterns The deformation pattern described in the previous section is merely one of the many kinds. In this paper we assume that the Fig. 6. Front and side views of the deformation sequence observed in Type I.
Fig. 3. The relation between u2 and Q of the planar deformation (1).
cross section of the beam is of oval shape. If we change the stiffness ratio l and the distance L, we may have ten different deformation patterns when l 1, as shown in the phase diagram Fig. 7(a). The deformation pattern described in the previous section is Type I. Fig. 7(a) is repeated in Fig. 7(b) with different scale. In the following, we describe briefly the behavior of each of these patterns. The analysis procedure is similar to the one in the previous section. For cases when l ¼ 1, the deformation sequences are verified by experiment. Not all cases are accompanied with photographs of the experiment in this paper. The interested reader may refer to the Master Thesis of the second author (Chen, 2014) for more details.
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
89
6.1. Type II In Fig. 8 we pull the two ends in Fig. 2 further away to L ¼ 0.85. In doing so the limit point on the path of 3D deformation (2) ceases to exist. As a consequence the planar symmetric deformation (1) jumps to another planar symmetric deformation (3) at bifurcation point B. The jump is accomplished via a 3D symmetric deformation shown as dashed curve in Fig. 8(b). In the experiment we record the dynamic jump by video camera and play it back slowly. The 3D jump can indeed be observed.
6.2. Type III In Fig. 9 we push the two ends in Fig. 2 closer to L ¼ 0.6. In this case there appear two spatial deformations instead of one. The planar deformation (1) evolves to a stable 3D symmetric deformation (2). At point S1, there emerges another 3D deformation (3) which is asymmetric with respect to the midpoint. After S1 both spatial deformations become unstable and the elastica jumps to planar symmetric deformation (4). Point S1 is a sub-critical bifurcation point. In Fig. 9(b) deformation (2) jumps to (4) via 3D asymmetric deformation (3), as shown by the dashed curve.
Fig. 7. (a) Deformation types in the lL plane; (b) in different scale.
Fig. 8. (a) Load-deflection diagram and (b) deformation evolution of type II (l ¼ 1, L ¼ 0.85).
Fig. 9. (a) Load-deflection diagram and (b) deformation evolution of type III (l ¼ 1, L ¼ 0.6).
90
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
Fig. 11. Two stable deformations (4) and (5) of type IV observed in lab.
Fig. 10. (a) Load-deflection diagram and (b) deformation evolution of type IV (l ¼ 1, L ¼ 0.45).
6.3. Type IV In Fig. 10 we push the two ends even closer to L ¼ 0.45. The planar symmetric deformation (1) evolves to a stable 3D symmetric deformation (2) at B. At point S there emerges a 3D asymmetric deformation (3). The elastica jumps at S, similar to the case in Type III. However, in Type IV there exists one additional 3D deformation (5) which is anti-symmetric with respect to the midpoint and twisted. Deformation (5) becomes stable after it merges with deformation (3) at Point C. As a consequence, there are two stable positions available after the jump, i.e., deformations (4) and (5). It is not an easy matter to predict which way to go after the jump without analyzing the potential energy structure, which is beyond the scope of this paper. In the experiment we found that the elastica jumps at point S from 3D symmetric deformation (2) to deformation (5) via 3D asymmetric deformation (3). Fig. 10(b) shows the deformation sequence. Although we are unable to reach deformation (4) by increasing Q quasi-statically, we can produce this position by hand. The photographs in Fig. 11 show that these two positions indeed exist in reality. In other words, they are both stable. It is noted that on the path of deformation (3) there exists a small stable segment between A1 and A2. Both A1 and A2 are limit
Fig. 12. (a) Load-deflection diagram and (b) deformation evolution of type V (l ¼ 1, L ¼ 0.3).
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
91
Fig. 14. (a) Load-deflection diagram and (b) deformation evolution of type VII (l ¼ 1, L ¼ 0.2).
Fig. 13. (a) Load-deflection diagram and (b) deformation evolution of type VI (l ¼ 1, L ¼ 0.25).
points. Because the range of Q is so narrow that we are unable to produce it in the experiment. 6.4. Type V Fig. 12 shows the case when we push the two ends further closer to L ¼ 0.3. At point B the planar symmetric deformation (1) bifurcates to a stable planar asymmetric deformation (2). At point S there emerges a 3D asymmetric deformation (3). There exists a small stable segment between A1 and A2 on the path of deformation (3), similar to Type IV in Fig. 10. Deformation (3) merges to an independent path of 3D deformation (5) which is anti-symmetric and twisted. Deformation (5) is stable beyond point C. The elastica jumps at point B from planar asymmetric deformation (2) to 3D twisted deformation (5) via 3D asymmetric deformation (3). 6.5. Type VI Fig. 13 shows the case when L ¼ 0.25. There is no spatial deformation in this case. The path of deformation (2) has a limit point S, as shown in Fig. 13(a). Deformation (2) is stable between B and S. Therefore, the planar symmetric deformation (1) bifurcates to planar asymmetric deformation (2) and then jumps at point S to symmetric deformation (3), as shown in Fig. 13(b). The elastica
crosses itself at the midpoint in deformation (3). It is noted that during snapping, the elastica moves sideway in the plane and crosses the pin joint. Due to the finite size of the pin joint in the experiment, the elastica will be blocked by the pin joint following the jump. 6.6. Type VII Fig. 14 shows the case when L ¼ 0.2. At point B the planar symmetric deformation (1) bifurcates to a planar asymmetric deformation (2). Both (1) and (2) are unstable beyond B. Deformation (2) merges at C to a planar symmetric deformation (3), which is stable beyond C. The elastica jumps at point B from deformation (1) to (3) via asymmetric deformation (2), as shown in Fig. 14(b). 6.7. Type VIII Up to this point, we limit our discussion to the case when l ¼ 1. If we set l equal to 2 and start from small L, the elastica behaves from Type VII to VI, and V. If we increase L to 0.6, we encounter a new type of deformation behavior, Type VIII, as shown in Fig. 15. If we compare the load-deflection diagram with the one in Fig. 9 (type III, for l ¼ 1 and L ¼ 0.6), we found that they are similar in structure, except that the 3D symmetric deformation (2) in Type III is replaced by a 2D asymmetric deformation (2) in Type VIII. Apparently, the elastica behaves more like a planar elastica as l increases. As shown in Fig. 15, deformation (3) starts at S1 and ends at S2 on the path of planar asymmetric deformation (2). The deformation sequence is depicted in Fig. 15(b), in which planar deformation (1) bifurcates to
92
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
planar asymmetric deformation (2) and then jumps to planar symmetric deformation (4) via a 3D asymmetric deformation (3). 6.8. Type IX If we increase L from 0.6 to 0.7, the spatial deformation disappears altogether. In Fig. 16 the planar symmetric deformation (1) bifurcates to planar asymmetric deformation (2) and then jumps to planar symmetric deformation (3). The difference between Types VIII and IX is that the elastica in Type VIII jumps via a spatial deformation while the elastica in Type IX jumps via a planar asymmetric deformation. 6.9. Type X For the case when L ¼ 0.9, the bifurcation point B occurs before the limit point, as shown in Fig. 17(a). Therefore, the planar symmetric deformation (1) jumps directly to planar symmetric deformation (3) via a planar asymmetric deformation (2), as depicted in Fig. 17(b). It is noted that the behaviors in Types IX and X are purely planar and have been studied previously by Chen and Hung (2011) with planar elastica model. 7. Conclusions
Fig. 15. (a) Load-deflection diagram and (b) deformation evolution of type VIII (l ¼ 2, L ¼ 0.6).
Fig. 16. (a) Load-deflection diagram and (b) deformation evolution of type IX (l ¼ 2, L ¼ 0.7).
In this paper we study the behavior of a spatial elastica under a midpoint force. The rod under question is buckled into a curved shape and then the two ends are fixed in space with pin joints. Both static and dynamic analyses are conducted. Focus is placed on the snapping phenomenon. The behavior of the elastica depends on two parameters; the distance between the two ends L and the bending stiffness ratio of the cross section l. A total of ten deformation types are identified. Their behaviors are summarized in Table 1, in which we use four characteristics to categorize the deformation patterns: (1) the deformation before jump, (2) the deformation after jump, (3) the type of critical point when jump
Fig. 17. (a) Load-deflection diagram and (b) deformation evolution of type X (l ¼ 2, L ¼ 0.9).
J.-S. Chen, R.-M. Chen / European Journal of Mechanics A/Solids 54 (2015) 84e93
93
Table 1 Classification of deformation behaviors. Type
Deformation before jump
Deformation after jump
Critical point
Jump via
I II III IV V VI VII VIII IX X
3D 2D 3D 3D 2D 2D 2D 2D 2D 2D
2D 2D 2D 3D 3D 2D 2D 2D 2D 2D
Limit Bifurcation Bifurcation Bifurcation Bifurcation Limit Bifurcation Bifurcation Limit Bifurcation
3D 3D 3D 3D 3D 2D 2D 3D 2D 2D
symmetric symmetric symmetric symmetric asymmetric asymmetric symmetric asymmetric asymmetric symmetric
symmetric symmetric symmetric twisted twisted symmetric (self-crossing) symmetric (self-crossing) symmetric symmetric symmetric
occurs, and (4) the deformation during the jump. It can be seen that no two deformation patterns are the same in all four characteristics. It is noted that types VI, VII, IX, and X exhibits only planar deformations. Therefore, these deformations can be obtained by using planar elastica model. By observing the phase diagram in Fig. 7 we can see that for the case when l is greater than 28, only these four types of deformation can exist in all ranges of L. We therefore conclude that if one wishes to design an elastica capable of only plane motion in all ranges of L, then l should be greater than 28. It is noted that a shallow planar buckled beam may buckle from a symmetric deformation to another symmetric deformation symmetrically (Chen and Hung, 2012). This occurs when the buckled beam is extensible and flat (L ¼ 1 5.6p2a). With the a in this article, this occurs when L > 0.9999. This type of snapping behavior will not be found in Fig. 6 because the rod in this paper is assumed to be inextensible here. The method presented in this paper is meant for a pinnedepinned elastica. However, the similar method can be applied to other types of boundary conditions, such as clampedeclamped spatial elastica. It is noted that the above conclusions are based on the numerical results on a specific specimen with G/E ¼ 0.376 and slenderness parameter a ¼ 1.767 106. In addition, the cross section is assumed to be of oval shape. If any of these are changed, some of the conclusions may be modified. For instance, if the cross section is changed to rectangular shape, the phase diagram look similar topologically, but the threshold l above which only plane deformation is possible will be changed to 37. From Eq. (11) we see that the slenderness parameter a affects only the inertial terms. Therefore, it does not affect the static equilibrium. In practical applications, the G/E ratios of most metals are in the range between 0.3 and 0.4. Although we did not test all the possible G/E ratios, it is expected that the overall structure of the phase diagram should remain the same in this range. References Cazottes, P., Fernandes, A., Pouget, J., Hafez, M., 2009. Bistable buckled beam: modeling of actuating force and experimental validations. ASME J. Mech. Des. 131, 101001. Chen, J.-S., Ro, W.-C., Lin, J.-S., 2009. Exact static and dynamic critical loads of a shallow arch under a point force at the midpoint. Int. J. Nonlin. Mech. 44, 66e70.
symmetric symmetric asymmetric asymmetric asymmetric asymmetric asymmetric asymmetric asymmetric asymmetric
Chen, J.-S., Hung, S.Y., 2011. Snapping of an elastica under various loading mechanisms. Eur. J. Mech. A Solid 30, 525e531. Chen, J.-S., Hung, S.-Y., 2012. Exact snapping loads of a buckled beam under a midpoint force. Appl. Math. Modell. 36, 1776e1782. Chen, J.-S., Tsao, H.-W., 2013. Static snapping load of a hinged extensible elastica. Appl. Math. Modell. 37, 8401e8408. Chen, R.-M., 2014. Deformation and Vibration of a Spatial Elastic under a Midpoint Force. Master Thesis. Department of Mechanical Engineering, National Taiwan University. in Chinese. Coleman, B.D., Dill, E.H., Lembo, M., Lu, Z., Tobias, I., 1993. On the dynamics of rods in the theory of Kirchhoff and Clebsch. Arch. Ration. Mech. Anal. 121, 339e359. Fung, Y.C., Kaplan, A., 1952. Buckling of Low Arches or Curved Beams of Small Curvature. NACA Technical Note 2840. Goriely, A., Tabor, M., 1997. Nonlinear dynamics of filaments, I. Dynamical instability. Phys. D 105, 20e44. Krylov, S., Ilic, B.R., Schreiber, D., Seretensky, S., Craighead, H., 2008. The pull-in behavior of electrostatically actuated bistable microstructures. J. Michromech. Microeng. 18, 055026. Kublanov, L.B., Bottega, W.J., 1995. On pressing a buckled film. Appl. Math. Modell. 19, 499e507. Maurini, C., Pouget, J., Vidoli, S., 2007. Distributed piezoelectric actuation of a bistable buckled beam. Eur. J. Mech. A Solid 26, 837e853. Medina, L., Gilat, R., Krylov, S., 2012. Symmetry breaking in an initially curved micro beam loaded by a distributed electrostatic force. Int. J. Solids Struct. 49, 1864e1876. Ouakad, H.M., Younis, M.I., 2014. On using the dynamic snap-through motion of MEMS initially curved microbeams for filtering applications. J. Sound Vib. 333, 555e568. Patricio, P., Adda-Bedia, M., Ben Amar, M., 1998. An elastica problem: instabilities of an elastic arch. Phys. D 124, 285e295. Pinto, O.C., Goncalves, P.B., 2000. Nonlinear control of buckled beams under step loading. Mech. Syst. Signal Process. 14, 967e985. Pippard, A.B., 1990. The elastic arch and its modes of instability. Eur. J. Phys. 11, 359e365. Plaut, R.H., 1978. Influence of load position on the stability of shallow arches. J. Appl. Math. Phys. (ZAMP) 30, 548e552. Reismann, H., Pawlik, P.S., 1980. Elasticity: Theory and Applications. John Wiley & Sons, New York. Schreyer, H.L., 1972. The effect of imperfections on the buckling load of shallow circular arches. ASME J. Appl. Mech. 445e450. Seide, P., 1984. Dynamic stability of laterally loaded buckled beams. J. Eng. Mech. 110, 1556e1569. Vangbo, M., 1998. An analytical analysis of a compressed bistable buckled beam. Sens. Actuat. A 69, 212e216. van der Heijden, G.H.M., Neukirch, S., Goss, V.G.A., Champneys, Thompson, J.M.T., 2003. Instability and self-contact phenomena in the writhing of clamped rods. Int. J. Mech. Sci. 45, 161e196. Virgin, L.N., Wiebe, R., Spottswood, S.M., Eason, T.G., 2014. Sensitivity in the structural behavior of shallow arches. Int. J. Nonlin. Mech. 58, 212e221. Zang, Y., Wang, Y., Li, Z., Huang, Y., Li, D., 2007. Snap-through and pull-in instabilities of arch-shaped beam under an electrostatic loading. J. Mechromech. Syst. 16, 684e693.