- Email: [email protected]
Numerical non-linear analysis of secondary buckling in stability problems
Computer methods in applied mechanics and engineerlng ELSEVIER
Numerical
Comput.
non-linear
Methods
Appl.
Mech.
Engrg.
120 (1995)
183-193
analysis of secondary buckling in stability problems Baisheng
Mathematics Department,
Wu
Jilin University. Changchun
Received
19 October
130023, China
1993
Abstract In this paper we discuss the numerical calculation of secondary buckling of structures. This kind of stability problem contains two parameters (one is conservative load parameter, the other is auxiliary one) and satisfies Zz-symmetry. By using the symmetry and introducing a new parameter we construct an augmented system. The system avoids the singularities of usual ones. Adjusting the newly introduced parameter, one can calculate a double buckling load and two corresponding buckling modes, secondary buckling load and secondary buckling mode directly without tracing primary post-buckling path. Here, the use of the augmented system yields a quadratically convergent iteration scheme. We also discuss the implementation of Newton’s method solving the system - a partitioning procedure. With this algorithm the amount of computation is largely saved.
1. Introduction
Rods, plates and shells as members, often appear in engineering structures. These members of structures are usually subjected to axial or membranous forces to give full play to load-bearing capability. Hence, bucklings of the structures have to be considered. In some structures, there are many different buckling modes, under the same buckling load; or artificially imposed approach or equality of buckling loads with different buckling modes, for example, optimization of structures. When these buckling loads are close or equal, the interactions among buckling modes make post-buckling behavior of structures more complex [l], for example, secondary buckling (secondary bifurcation in buckling problem) which means a new buckling (bifurcation) appears in a primary post-buckling path except trivia1 (undeformed configuration) equilibrium path. For the definition of secondary bifurcation, we may also refer to [2] for details. Once a load arrives a buckling one, the stability behavior of a structure will change. This make the investigation of secondary buckling of structures more important. Early in 1975, Bauer et al. [2] proposed ‘multiple eigenvalues lead to secondary bifurcation’ by using a buckling mode1 of rod of two degrees of freedom with two parameters. They also gave a perturbation method to study it. By using this method, a few authors found secondary buckling phenomenon in some structures [3-51. Going over the literature above, one easily discovers that secondary bucklings are generally related to stability problems with &-symmetry and appear in symmetric primary postbuckling paths. Based on the &-symmetry, local and qualitative analyses of secondary bifurcation (buckling) have been given by Shearer [6], Potier-Ferry [7], Golubitsky and Schaeffer [S] by using Liapunov-Schmidt reduction and normal form theory, respectively. Because of diverse and complicated computations of the coefficients of these analyses, their approaches are only valid to simple structures, yet do not apply to numerical calculation. Also, a secondary buckling load might be arbitrarily close to a primary buckling load, it will be difficult to detect the existence of secondary buckling by the continuation method because of trouble tracing in itself primary post-buckling path near the primary buckling load; moreover, even if one ascertains the existence of a secondary buckling point, usual 0045.7825/95/$09.50 0 1995 Elsevier SSDI 004S-7825(94)(~Oo55-R
Science
S.A.
All rights
reserved
184
B. Wu I Comput. Methods Appl. Mech. Engrg.
120 (1995)
183-193
augmented systems [9, lo] computing bifurcation points in non-trivial solution branch near a double eigenvalue are singular or almost singular. Chien [ll] has used a local perturbation skill to deal with the computation of secondary post-buckling paths. However, it is very difficult to choose the perturbation vector d for performing the local perturbation, yet he did not point out how to calculate a secondary buckling point exactly. In this paper, we investigate the numerical calculation of secondary buckling problems satisfying Z,-symmetry and containing two parameters. A Newton-type method for direct calculations of double buckling points and secondary buckling points will be presented. The basis for this procedure is an extension of the non-linear set of equations by constraint conditions which introduces additional information about these points-so-called augmented systems. The discussions of computing a simple limit and bifurcation point of non-linear stability problems by an augmented system may be found in [9, lo]. Note that these techniques lead to a Newton scheme with quadratic convergence behavior. Furthermore, augmented system provides enough information to branch-switching. As compared to existing work, our augmented system introduces a new parameter except for the original two parameters. It brings the double buckling point and secondary buckling point into one system. By adjusting a newly introduced parameter, one can calculate the double buckling load and corresponding two buckling modes; secondary buckling load and secondary buckling mode directly without through tracing the primary post-buckling path. Hence, we can use another augmented system to compute secondary post-buckling path. Two numerical examples in secondary buckling problems will be used to account for the effectiveness of the method.
2. Calculation of secondary buckling via augmented system Consider a general structural stability problem containing two parameter,
viz. (1)
f(& A, CL)=0
where f:R” x R2+Rn is a smooth, column-vector-valued function. The column vector x E R”’ represents the state variables or basic unknowns of the problem, and the scalar A E IF&! is a conservative It is assumed that (1) arises from a model of a load parameter, and p E IF8is an auxiliary one. mechanical system with finitely many degrees of freedom, or from a finite-dimensional approximation (e.g. by finite element method) of a continuous system. Suppose that the total potential energy of the system is V(X, A, /L), v : RI” x R2+ R”
(2)
and satisfies the following Z,-symmetry: SZI,
s2=z,
V(Sx, A, CL)= I+,
x N
matrix S, such that (3a)
A, P)
where I is N x N identity operator. that f in (1) satisfies f(&
there is a N
A, CL)= V(& A, EL).
(3b) By the principle of stationary potential energy and (3b), we obtain (4)
Note that the Z,-symmetry appears to be the most common in physical examples. In the following discussions, we make the hypotheses: (H,). f(0, A, p) = 0 holds for any (A, CL)E R2. This means that for any load, the undeformed state is the trivial equilibrium path. (H2). When p = pO, the buckling load A = A, is a double one, and the corresponding two buckling modes 4,) +2 satisfy, respectively
B. Wu I Comput. Methods Appl. Mech. Engrg.
120 (1995)
183-193
185
The point (x, A, p) = (0, A,, pO) is called double buckling point. Note that these two assumptions have appeared in [6-S] for the purpose of theoretic investigations. f: := f,(O, A,, pO), ‘t’ Throughout this paper, we make notations f,(x, A, CL):= #(x, A, p)lax, represents transposition, etc. By (3a) we have natural decomposition rw”=@%rw~,
(6a)
lR~={xERN,SX=X})
(6b)
where
rw: = {x E RN, sx = -x}
)
(64
and f maps R,N into itself. The rWr and Iwr are called symmetric and antisymmetric subspace, respectively, of RN. a secondary buckling appears in the symmetric primary As pointed out in the Introduction, post-buckling path. Therefore, we first give an augmented system detecting the symmetric primary post-buckling path which bifurcate from the primary buckling point in the trivial equilibrium path as follows
here, E, : R,” x IF8x lR*-+ Rf” x [w, and if c#O
&(x, t A, CL,&I=
C’b)
if e=O . Let E = 0 in (7), then it becomes the augmented system locating the primary buckling point. By using the implicit function theorem, one can prove [12] so long as a, =
+:fxodP, +o >
(8)
then for every p near pO, the symmetric primary post-buckling A = A,, may uniquely be presented by E: I;, = {x = EX”~(E, II), A = h”((~,p) : I&I< Ed)
path near the double buckling load
(9)
where X”,(E,p) E II??,“. Because secondary buckling points are closely related to double buckling points, we first use the following augmented system to calculate the double buckling load and corresponding two buckling modes: [.&CO>A, Y)x,~ o &(x,,
x2,
A,
PU) =
=
(10)
where E, : $” x rWr x Pi*--+ rWr x II: x R*, x1 and x2 represent, respectively, the symmetric and antisymmetric buckling modes. It is easy to know 1121that the Jacobi of E, at the point (x1, x2, A, P) = (4,) +2, A,,, p(,) is non-singular provided the determinant
4:f:A4,
4:f:,4
4if4*42
4ifi,42
z. .
So one can use Newton’s method to calculate the point.
(11)
B. Wu I Comput. Methods Appl. Mech. Engrg.
186
120 (1995) 183-193
REMARK
2.1. A non-singular system similar to (10) has been proposed by Werner [13] for the detection of a double singular point of (l), however, it is not in combination with a conservative system. We also examine the computation of secondary buckling and give an efficient implementation of Newton’s method on system (10).
Suppose one has located the double buckling point by augmented system (10). To compute the secondary buckling point that occurred in the symmetric primary post-buckling path above, instead of usual continuation method, we introduce a new parameter S in (lo), and use the augmented system
if s F1@1,~2,
A,
~6)
f
o
(124
=
if6=0’
(124 to directly calculate the secondary buckling load and corresponding buckling mode by continuating 6 from 6 = 0 and (x1, x2, A, P) = (~&,4~, A,, P”). In fact, set S = 0 in (12), then it changes to (10). Using the non-singularity of the Jacobi of (10) at the point (x, , x2, A, p) = (4,) c$~,A,, p”) and the implicit function theorem, we can represent the solutions of (11) in terms of 6
(X,(O),
X2(0),
h(O),
CL(O))
where f,(S) E Rf”, X2(6)~Rr. @,(6),
A@), 4 = 0.
=
(4,427
A,>
14
For any given 6 f0,
( 13b)
let 6 = b(6).
By (12a) and (12c), one has (14)
Using the uniqueness of symmetric primary post-buckling path when p = G(6) (see (9)), we see that (6X,(6), A(6)) is in & and 6X,(S) #O. Again, from (12a) and (12d) one knows that there is a antisymmetric eigenvector X2(8) at the point @X,(S), A(6)) of the r,,. So when P = b(S), (x, A) = (6X, (6 ), h(6 )) is just a secondary buckling point and X;(S ) is the corresponding secondary buckling mode. Augmented systems (10) and (12) are associated with N + 2 unknowns. The numerical effort to solve such systems may again be decreased by using a partitioning method (see Section 3). Once the secondary buckling point and corresponding antisymmetric post-buckling mode have been calculated by (12), we can use another augmented system (fixed 6 # 0):
( 15b) (15c)
B. Wu I Compu~. Methods Appl. Mech. Engrg.
_
I
120 (109.5) 183-193
187
if y=O,
f,(T A, I.@>M
to trace the initial secondary post-buckling branch by continuating y and starting from y = 0 and (x, $, A) = (6X,(6), X*(S), A(6)). When y is appropriately large, one may continuate the load parameter A to compute secondary post-buckling path by existing continuation methods [14].
3. Implementation
of Newton’s
method solving augmented
system (10)
In this section, we take Newton’s method solving augmented system (10) as an example to illustrate its implementation - a partitioning method. Of course, this method will also be valid to augmented system (12) so long as 6 is proper small. By (6), one has N = N, + N,
(16)
where N, = dim Rf”, N, = dim Rr. In general, it will be easy to identify R~(rW~) with R”S(R”a) using the isomorphism Z,: Rr--+ RN,, Z,(x) = x,(Z, : rWt-+ lRNd, Z,(x) = x,). Using this isomorphism, augmented system (10) may be changed into
(17)
where & : RN5 x RN, x R* + RNs x RNa x R2. The incremented leads to the following N + 2 dimensional linear systems:
solution scheme using Newton’s method
A,Ax,+B,AA+C,Ap=-r,
(Isa)
A,Ara+B,AA+C,Ap=-r,
(lgb)
x:AxS=
-r SC
(184
xi5 Ax,, = -r,,
(184
where
A, = Lf;f;(O,A, /4,’
,
A, = I,f,(O> A, P>I,’
,
B, = 4f,,(O, 4 /-+,‘x, > B, = z,f,,(O, A, PV,‘X, >
Cs= 4f,,(% At /-X1x, C, = tJJ0,
7
A, /-4Z,‘x,,
ra = LfJO, A, PVC’X,, ra = LL(O, A, PK~x, , rsr = (x:x, - 1)/2 , rap = (xix, - 1)/2 .
188
B. Wu I Comput. Methods Appl. Mech. Engrg. 120 (1995) 183-193
The tangent matrix in (18) is not fully filled, therefore we may use a special elimination technique which comes out to be a very efficient way to solve (18) within one Newton iteration step. Step 1. Solve
Step 2. Calculate D,=Ah’B,+C, R, = Ah2 B, + Y, Step 3. Solve
Step 4. Calculate the new increments
of AX,, AA:
As can be seen, the matrixes which have to be factorized in this solution process are the coefficient ones in Steps 1 and 3 with dimensions Ns + 1 and N, + 1, respectively. Thus, the cost of using (18), specially for large systems, is decreased considerably (note that N, + N, = N) REMARK 3.1. At the double buckling point, the matrix A, and A, are all singular. But the coefficient matrix of Step 1 is non-singular under the condition (8); furthermore, the coefficient matrix of Step 3 is also non-singular provided the condition (11) holds. The hypotheses in (8) and (11) are common in the theoretical investigations of secondary bifurcation (see [6-g]). Hence, both the two matrixes may be factored without ill-conditioning and our method is quadratically convergent.
To compute the double buckling point approximation for ($,, &, h,, pO) where +, = I,41 >
4,, = Ia42 ’
with Newton’s
method,
one has to have
a suitable
(19)
Note that the double buckling load and two corresponding buckling modes are just the eigenvalues and corresponding eigenvectors of the operator equation linearized from the trivial solution path x = 0 of (1) (see [lOI) Ax = ABx
(20)
where A and B are N x N symmetric matrixes and A positive definite and p is implicitly contained in A and B. For given CL,by using the isomorphisms Z, and I,, we can solve the two smallest buckling loads with, respectively, symmetric and antisymmetric buckling modes from (20) (see e.g. 1151). If the two buckling loads are near, then we may respectively take the isomorphisms of the two corresponding buckling modes, one of the two buckling loads, and I_Las a initial approximation to (4,, 4,, A,, po).
B. Wu I Comput.
Merhods Appl. Mech. Engrg.
189
120 (199-i) 183-lY3
4. Numerical examples If this section, we would like to present two numerical examples which will illustrate the characteristics of the algorithms developed. These examples include computation of double buckling points as well as secondary buckling points. All computations have been performed on a SUNworkstation. EXAMPLE 4.1. Coupled elastic rods as a simple system demonstrating secondary buckling [ 161. Two rigid rods are connected by a springhinge. The bottom rod is constrained by another hinge and spring, and bends in the plane y = 0 through an angle /3. The top rod bends through an angle a in the plane containing the bottom rod and the y-axis. The hinges and springs constrain the rods to buckle in mutually orthogonal planes. The top rod is subjected a vertical compressive load p. The total potential energy of the system is la2+;k,~2+p(l,+l,cosa)cos~
V(a,p,p,k,,k,.l,,fZ)=~k
(21)
where (Y, p = bending angles of top and bottom rods, 1,) I, = rod lengths, k,, k, = spring stiffnesses. Introducing non-dimensional variables y=---k, k, ’
p+, 2
*=E!? k2
one obtains non-dimensional
V=& ’
2
potential energy
V(ol,p,n,~,U)=~~12*+~~*+A(1+~coSa)cos~.
(22)
Define
s=(:, -:>
(23)
and x = (a, p)‘, then the potential energy in (22) satisfies (3a) and (3b) with fixed V. By using stationary potential principle, we derive equilibrium equations of the system ~a - AE.L sin (Ycos p = 0,
(244
p - A( 1 + p cos a!) sin p = 0 .
(24b)
From (23), we also have
I,=(1
01,
In the following becomes
z-1 = 5
0
Z,=(O l),
;
numerical
computation,
I,‘=
001
we let v = 0.25. Then the augmented
(25) system in (12)
ACL sin@4
1. 1 o
25a
_
P-hll+CLco&LI),P *-1 ;*-1
=O,
(26)
and the augmented system (7) (set E = 0) for detecting the primary buckling load related symmetric primary post-buckling branch becomes 0.25a - h/_&cX 0. (Y2-l 1 =
1
At a double buckling point of the system, (A, &, A,, F,,) = (1, 1, f , f). Using augmented
(27) systems
190
B. Wu
Table 1 Dependence
I Comput.
Methods
Appl.
Mech.
Engrg.
120 (1995)
of secondary buckling of coupled rods on p
E
IL(&)
A,,(E)
k(E)
4)
0
0.33333333 0.33336297 0.33345193 0.33407721 0.33634673 0.34601193 0.35390944 0.36437317 0.37798559 0.39554873 0.41818119 0.44747901
0.75000000 0.74993332 0.74973325 0.74833000 0.74328060 0.72251844 0.70639540 0.68610979 0.66140088 0.63203338 0.59782699 0.55868542
0.75000000 0.75013335 0.75053356 0.75334226 0.76347776 0.80575708 0.83947685 0.88336613 0.93963315 1.01168479 1.10498659 1.22882878
0.00000000 0.04000000 0.08000000 0.20000000 0.40000000 0.80000000 1.00000000 1.20000000 1.40000000 I .60000000 1.80000000 2 .oooooooo
0.04 0.08 0.20 0.40 0.80 1.00 1.20 1.40 1.60 1.80 2.00
1113-193
(26) resp. (27) t o continuate 6 resp. p with the methods in Section 3. We get Table 1 which contains the dependences of the primary buckling load A,, secondary buckling load A, and corresponding bending angle (Y, of the top rod on the length ratio p.
EXAMPLE
4.2.A simply supported
compression.
The normal displacement
d4w + h d’w ~+~w+w~=o, dx4
beam on a non-linear elastic foundation subjected to axial w(x) of the beam satisfies the following equation [7]:
(284
o
(28b) where A is proportional to the axial compression and pw + w3 is the reaction of the foundation. It is easy to check that the solution of (28) is just the stationary point of the following potential energy function under the conditions w(0) = w(n) = 0:
(29) We also see that u satisfies the Z,-symmetry U(S,w, A, I*) = VW, A, P)
(30a)
s, w(x) = w(7T-x)
W’b)
where )
x E [O, Tr] .
It is instructive to describe in detail how a discretization of (28) leads to a finite dimensional problem involving symmetry analogous to (28). To fully utilize the symmetry we choose a uniform mesh size, h = nrl(N + l), with N = 2n, II: integer, and nodes x, = ih, i = 1, 2, . , . , N. Starting from (29) and using difference-variation (evaluate the integral with trapezoid formula and use central difference quotient approximations to d*wldX* and dwldx) discretization to (28) leads to finite dimensional problem (1) where x = (w,, w2,. . . , w~)~ is the column vector of the nodes and the components of f are 7
f, =w3 - 4w, - (5 + h4p + h4w;)w,
-;
At-w,
+ 3w,) ,
f; = w,_2 - 4W,_l - 4w,+, + wj+2 + (6 + h4F + h4w;)wI 2SiaN-l,w,=w
Iv+1 = 0,
h2 - 4 A(2w, - wi_2 -
Wi+*)
)
B. Wu I Comput. Methods Appl. Mech. Engrg.
120 (1995)
183-193
191
7
+ (5 + h4p + h4w;)w,
fN = wN_* - 4~,~_,
+ % A(w,_,
- 3w,)
.
Define 0
1
s=
1 :. C1
E[W2nx2n,
(31)
01
Then the S and f above, satisfy (3a) and (4). Now we have I, =;(I,,
E,)+4,
Z,‘=(I,
E,)‘,
r,’
.*. .f
:)Wfl,
)
I, = ; (Z, = (I,
(324
-En)+!,,
Wb)
-En)’
where 1;-6
I,, I, E RnxZn ; Then the augmented
I,‘,
&=(Y
I,’
l
..* ;)ERY
(32~)
E lRZnxn .
system (12) is equivalent to the following system
(33)
By continuating p, the corresponding following augmented system
symmetric primary post-buckling load can been calculated by the
@u-JO~A, PPCW,= o
1
w;w,- 1
I .
(34)
At a double buckling point of the original continuous problem (28), (35) Let N = 80(n = 40). By using the method in Section 3, we get the dependences of the symmetric primary buckling load h,, the secondary buckling load A, and corresponding maximum deflection wraX of the beam on the stiffness of the linear part of the reaction or. as following Table 2. From Table 2, one knows that as the deviation of p from pCJ,the secondary buckling load slightly departs from the symmetric primary buckling load, yet the maximum deflection of the beam increase rapidly.
4. Conclusions An augmented system containing a newly introduced parameter has been given for the numerical calculation of secondary buckling in non-linear stability problems. The implementation of Newton’s method solving the resulted system has been developed by a partitioning method. The method directly yields secondary buckling load and corresponding secondary buckling mode without tracing primary post-buckling branch. The eigenmode associated with secondary buckling point can be used for
B. Wu I Comput. Methods Appl. Mech. Engrg.
192 Table
I20 (1995) 183-193
2
Dependence
of secondary
buckling
of beam on Jo
f
CL(E)
,$A)
A,(E)
Wyy&)
0 0.5 1.0 1.5 2.0 3.0 4.0 5.0 5.5 6.0 6.5 7.0
4.00502046 3.99885508 3 .Y803742Y 3.94962411 3.90668192 3.18469338 3.61569885 3.40159514 3.27841124 3.14499291 3.00180535 2.84938408
5.00728033 5.00111864 4.98262180 4.95185619 4.90889246 4.78684272 4.61776342 4.40355229 4.28030659 4.14682133 4.00356194 3.85106421
5.00728033 5.01037506 5.01966595 5.03517321 5.05693084 5.11940548 5.20766262 5.32255394 5.39034378 5.46526926 5.54754715 5.63742969
0.00000000 0.11108165 0.22211179 0.33303849 0.44380898 0.66466313 0.88421666 1.101Y6390 1.20998301 1.31732600 1.42390320 1.52961495
branch-switching. A consistent linearization of the governing equations leads to quadratic convergence behavior. The two numerical examples have displayed that one can follow the secondary buckling points sufficiently far away from the double buckling points so that the further computation of the points can be carried out by using a continuation procedure based on the augmented system [9]. The algorithm has been applied to coupled elastic rods and a simply supported beam in a non-linear elastic foundation. In both cases the algorithm is very successful in computing double buckling points and secondary buckling points and shows good convergence.
Acknowledgment
This work has been financially supported by the National Science Foundation of China. The author would like to thank Professor T. Kiipper and the Mathematisches Institut der Universitat zu Koln for their hospitality. This paper was finished while the author visited Koln, where he was supported by Volkswagen-Stiftung, and he would like to thank it for its support.
References
111H.
Troger, Application of bifurcation theory to the solution of nonlinear stability problems in mechanical engineering, in: T. Kiipper, H.D. Mittelmann and H. Weber, eds., Numerical Methods for Bifurcation Problems, ISNM 70 (Birkhauser Verlag, 1984). t4 L. Bauer, H.B. Keller and E.L. Reiss, Multiple eigenvalues lead to secondary bifurcation, SIAM Rev. 17 (1975) 101-122. L.J. Putnick and E.L. Reiss, Secondary states of rectangular plates, SIAM J. Appl. Math. 38 (1980) 131 B.J. Matkowsky, 38-51. Secondary buckled state of rectangular sandwich plates, Acta Mechanica Sinica 6 (1990) 141 He Luwu and Cheng Changjun, 237-247. Secondary bifurcation and imperfection sensitivity of cylindrical shell of finite length, Internat. J. [51 M.O. Oyesanya, Non-Linear Mech. 25 (1990) 405-415. SIAM J. Math. Anal. 11 (1980) 365-389. PI M. Shearer, Secondary bifurcation near a double eigenvalues, Multiple bifurcation, symmetry and secondary bifurcation, in: Nonlinear Problems of Analysis in Geometry f71 M. Potier-Ferry, and Mechanics, Research Notes in Math. 46 (Pitman, 1981) 158-167. and D. Schaeffer, Singularities and Groups in Bifurcation Theory, Vol. 1 (Springer-Verlag, 1985). PI M. Golubitsky of symmetry-breaking bifurcation points, SIAM J. Numer. Anal. 21 (1984) 191 B. Werner and A. Spence, The computation 388-399. convergent procedure for the calculation of stability points in finite [lOI P. Wriggers, W. Wagner and C. Miehe, A quadratically element analysis, Comput. Methods Appl. Mech. Engrg. 70 (1988) 329-347. [Ill C.S. Chien, Secondary bifurcation in the buckling problem, J. Comput. Appl. Math. 25 (1989) 277-287. [121 B. Wu, The compution of secondary bifurcation and its applications to buckling of structures, Ph.D. Thesis, Jilin University, Changchun, 1991.
6. Wu / Comput. Methods Appl. Mech. Engrg. [13] B. Werner,
120 (1995) 183-193
193
Regular systems for bifurcation points with underlying symmetries, in: T. Kiipper, H.D. Mittelmann and H. Weber, eds., Numerical Methods for Bifurcation Problems, ISNM 70 (Birkhauser Verlag, 1984) 562-574. [14] H.B. Keller, Lectures on Numerical Methods in Bifurcation Problems (Springer-Verlag, 1987). [15] G.H. Golub and C.F. Van Loan, Matrix Computation, 2nd edition (The Johns Hopkins University Press, Baltimore, Maryland, 1989). [16] S.J. Tavener and T. Mullin, Buckling of coupled elastic rods, Physica D 30 (1988) 382-398.
Numerical
Comput.
non-linear
Methods
Appl.
Mech.
Engrg.
120 (1995)
183-193
analysis of secondary buckling in stability problems Baisheng
Mathematics Department,
Wu
Jilin University. Changchun
Received
19 October
130023, China
1993
Abstract In this paper we discuss the numerical calculation of secondary buckling of structures. This kind of stability problem contains two parameters (one is conservative load parameter, the other is auxiliary one) and satisfies Zz-symmetry. By using the symmetry and introducing a new parameter we construct an augmented system. The system avoids the singularities of usual ones. Adjusting the newly introduced parameter, one can calculate a double buckling load and two corresponding buckling modes, secondary buckling load and secondary buckling mode directly without tracing primary post-buckling path. Here, the use of the augmented system yields a quadratically convergent iteration scheme. We also discuss the implementation of Newton’s method solving the system - a partitioning procedure. With this algorithm the amount of computation is largely saved.
1. Introduction
Rods, plates and shells as members, often appear in engineering structures. These members of structures are usually subjected to axial or membranous forces to give full play to load-bearing capability. Hence, bucklings of the structures have to be considered. In some structures, there are many different buckling modes, under the same buckling load; or artificially imposed approach or equality of buckling loads with different buckling modes, for example, optimization of structures. When these buckling loads are close or equal, the interactions among buckling modes make post-buckling behavior of structures more complex [l], for example, secondary buckling (secondary bifurcation in buckling problem) which means a new buckling (bifurcation) appears in a primary post-buckling path except trivia1 (undeformed configuration) equilibrium path. For the definition of secondary bifurcation, we may also refer to [2] for details. Once a load arrives a buckling one, the stability behavior of a structure will change. This make the investigation of secondary buckling of structures more important. Early in 1975, Bauer et al. [2] proposed ‘multiple eigenvalues lead to secondary bifurcation’ by using a buckling mode1 of rod of two degrees of freedom with two parameters. They also gave a perturbation method to study it. By using this method, a few authors found secondary buckling phenomenon in some structures [3-51. Going over the literature above, one easily discovers that secondary bucklings are generally related to stability problems with &-symmetry and appear in symmetric primary postbuckling paths. Based on the &-symmetry, local and qualitative analyses of secondary bifurcation (buckling) have been given by Shearer [6], Potier-Ferry [7], Golubitsky and Schaeffer [S] by using Liapunov-Schmidt reduction and normal form theory, respectively. Because of diverse and complicated computations of the coefficients of these analyses, their approaches are only valid to simple structures, yet do not apply to numerical calculation. Also, a secondary buckling load might be arbitrarily close to a primary buckling load, it will be difficult to detect the existence of secondary buckling by the continuation method because of trouble tracing in itself primary post-buckling path near the primary buckling load; moreover, even if one ascertains the existence of a secondary buckling point, usual 0045.7825/95/$09.50 0 1995 Elsevier SSDI 004S-7825(94)(~Oo55-R
Science
S.A.
All rights
reserved
184
B. Wu I Comput. Methods Appl. Mech. Engrg.
120 (1995)
183-193
augmented systems [9, lo] computing bifurcation points in non-trivial solution branch near a double eigenvalue are singular or almost singular. Chien [ll] has used a local perturbation skill to deal with the computation of secondary post-buckling paths. However, it is very difficult to choose the perturbation vector d for performing the local perturbation, yet he did not point out how to calculate a secondary buckling point exactly. In this paper, we investigate the numerical calculation of secondary buckling problems satisfying Z,-symmetry and containing two parameters. A Newton-type method for direct calculations of double buckling points and secondary buckling points will be presented. The basis for this procedure is an extension of the non-linear set of equations by constraint conditions which introduces additional information about these points-so-called augmented systems. The discussions of computing a simple limit and bifurcation point of non-linear stability problems by an augmented system may be found in [9, lo]. Note that these techniques lead to a Newton scheme with quadratic convergence behavior. Furthermore, augmented system provides enough information to branch-switching. As compared to existing work, our augmented system introduces a new parameter except for the original two parameters. It brings the double buckling point and secondary buckling point into one system. By adjusting a newly introduced parameter, one can calculate the double buckling load and corresponding two buckling modes; secondary buckling load and secondary buckling mode directly without through tracing the primary post-buckling path. Hence, we can use another augmented system to compute secondary post-buckling path. Two numerical examples in secondary buckling problems will be used to account for the effectiveness of the method.
2. Calculation of secondary buckling via augmented system Consider a general structural stability problem containing two parameter,
viz. (1)
f(& A, CL)=0
where f:R” x R2+Rn is a smooth, column-vector-valued function. The column vector x E R”’ represents the state variables or basic unknowns of the problem, and the scalar A E IF&! is a conservative It is assumed that (1) arises from a model of a load parameter, and p E IF8is an auxiliary one. mechanical system with finitely many degrees of freedom, or from a finite-dimensional approximation (e.g. by finite element method) of a continuous system. Suppose that the total potential energy of the system is V(X, A, /L), v : RI” x R2+ R”
(2)
and satisfies the following Z,-symmetry: SZI,
s2=z,
V(Sx, A, CL)= I+,
x N
matrix S, such that (3a)
A, P)
where I is N x N identity operator. that f in (1) satisfies f(&
there is a N
A, CL)= V(& A, EL).
(3b) By the principle of stationary potential energy and (3b), we obtain (4)
Note that the Z,-symmetry appears to be the most common in physical examples. In the following discussions, we make the hypotheses: (H,). f(0, A, p) = 0 holds for any (A, CL)E R2. This means that for any load, the undeformed state is the trivial equilibrium path. (H2). When p = pO, the buckling load A = A, is a double one, and the corresponding two buckling modes 4,) +2 satisfy, respectively
B. Wu I Comput. Methods Appl. Mech. Engrg.
120 (1995)
183-193
185
The point (x, A, p) = (0, A,, pO) is called double buckling point. Note that these two assumptions have appeared in [6-S] for the purpose of theoretic investigations. f: := f,(O, A,, pO), ‘t’ Throughout this paper, we make notations f,(x, A, CL):= #(x, A, p)lax, represents transposition, etc. By (3a) we have natural decomposition rw”=@%rw~,
(6a)
lR~={xERN,SX=X})
(6b)
where
rw: = {x E RN, sx = -x}
)
(64
and f maps R,N into itself. The rWr and Iwr are called symmetric and antisymmetric subspace, respectively, of RN. a secondary buckling appears in the symmetric primary As pointed out in the Introduction, post-buckling path. Therefore, we first give an augmented system detecting the symmetric primary post-buckling path which bifurcate from the primary buckling point in the trivial equilibrium path as follows
here, E, : R,” x IF8x lR*-+ Rf” x [w, and if c#O
&(x, t A, CL,&I=
C’b)
if e=O . Let E = 0 in (7), then it becomes the augmented system locating the primary buckling point. By using the implicit function theorem, one can prove [12] so long as a, =
+:fxodP, +o >
(8)
then for every p near pO, the symmetric primary post-buckling A = A,, may uniquely be presented by E: I;, = {x = EX”~(E, II), A = h”((~,p) : I&I< Ed)
path near the double buckling load
(9)
where X”,(E,p) E II??,“. Because secondary buckling points are closely related to double buckling points, we first use the following augmented system to calculate the double buckling load and corresponding two buckling modes: [.&CO>A, Y)x,~ o &(x,,
x2,
A,
PU) =
=
(10)
where E, : $” x rWr x Pi*--+ rWr x II: x R*, x1 and x2 represent, respectively, the symmetric and antisymmetric buckling modes. It is easy to know 1121that the Jacobi of E, at the point (x1, x2, A, P) = (4,) +2, A,,, p(,) is non-singular provided the determinant
4:f:A4,
4:f:,4
4if4*42
4ifi,42
z. .
So one can use Newton’s method to calculate the point.
(11)
B. Wu I Comput. Methods Appl. Mech. Engrg.
186
120 (1995) 183-193
REMARK
2.1. A non-singular system similar to (10) has been proposed by Werner [13] for the detection of a double singular point of (l), however, it is not in combination with a conservative system. We also examine the computation of secondary buckling and give an efficient implementation of Newton’s method on system (10).
Suppose one has located the double buckling point by augmented system (10). To compute the secondary buckling point that occurred in the symmetric primary post-buckling path above, instead of usual continuation method, we introduce a new parameter S in (lo), and use the augmented system
if s F1@1,~2,
A,
~6)
f
o
(124
=
if6=0’
(124 to directly calculate the secondary buckling load and corresponding buckling mode by continuating 6 from 6 = 0 and (x1, x2, A, P) = (~&,4~, A,, P”). In fact, set S = 0 in (12), then it changes to (10). Using the non-singularity of the Jacobi of (10) at the point (x, , x2, A, p) = (4,) c$~,A,, p”) and the implicit function theorem, we can represent the solutions of (11) in terms of 6
(X,(O),
X2(0),
h(O),
CL(O))
where f,(S) E Rf”, X2(6)~Rr. @,(6),
A@), 4 = 0.
=
(4,427
A,>
14
For any given 6 f0,
( 13b)
let 6 = b(6).
By (12a) and (12c), one has (14)
Using the uniqueness of symmetric primary post-buckling path when p = G(6) (see (9)), we see that (6X,(6), A(6)) is in & and 6X,(S) #O. Again, from (12a) and (12d) one knows that there is a antisymmetric eigenvector X2(8) at the point @X,(S), A(6)) of the r,,. So when P = b(S), (x, A) = (6X, (6 ), h(6 )) is just a secondary buckling point and X;(S ) is the corresponding secondary buckling mode. Augmented systems (10) and (12) are associated with N + 2 unknowns. The numerical effort to solve such systems may again be decreased by using a partitioning method (see Section 3). Once the secondary buckling point and corresponding antisymmetric post-buckling mode have been calculated by (12), we can use another augmented system (fixed 6 # 0):
( 15b) (15c)
B. Wu I Compu~. Methods Appl. Mech. Engrg.
_
I
120 (109.5) 183-193
187
if y=O,
f,(T A, I.@>M
to trace the initial secondary post-buckling branch by continuating y and starting from y = 0 and (x, $, A) = (6X,(6), X*(S), A(6)). When y is appropriately large, one may continuate the load parameter A to compute secondary post-buckling path by existing continuation methods [14].
3. Implementation
of Newton’s
method solving augmented
system (10)
In this section, we take Newton’s method solving augmented system (10) as an example to illustrate its implementation - a partitioning method. Of course, this method will also be valid to augmented system (12) so long as 6 is proper small. By (6), one has N = N, + N,
(16)
where N, = dim Rf”, N, = dim Rr. In general, it will be easy to identify R~(rW~) with R”S(R”a) using the isomorphism Z,: Rr--+ RN,, Z,(x) = x,(Z, : rWt-+ lRNd, Z,(x) = x,). Using this isomorphism, augmented system (10) may be changed into
(17)
where & : RN5 x RN, x R* + RNs x RNa x R2. The incremented leads to the following N + 2 dimensional linear systems:
solution scheme using Newton’s method
A,Ax,+B,AA+C,Ap=-r,
(Isa)
A,Ara+B,AA+C,Ap=-r,
(lgb)
x:AxS=
-r SC
(184
xi5 Ax,, = -r,,
(184
where
A, = Lf;f;(O,A, /4,’
,
A, = I,f,(O> A, P>I,’
,
B, = 4f,,(O, 4 /-+,‘x, > B, = z,f,,(O, A, PV,‘X, >
Cs= 4f,,(% At /-X1x, C, = tJJ0,
7
A, /-4Z,‘x,,
ra = LfJO, A, PVC’X,, ra = LL(O, A, PK~x, , rsr = (x:x, - 1)/2 , rap = (xix, - 1)/2 .
188
B. Wu I Comput. Methods Appl. Mech. Engrg. 120 (1995) 183-193
The tangent matrix in (18) is not fully filled, therefore we may use a special elimination technique which comes out to be a very efficient way to solve (18) within one Newton iteration step. Step 1. Solve
Step 2. Calculate D,=Ah’B,+C, R, = Ah2 B, + Y, Step 3. Solve
Step 4. Calculate the new increments
of AX,, AA:
As can be seen, the matrixes which have to be factorized in this solution process are the coefficient ones in Steps 1 and 3 with dimensions Ns + 1 and N, + 1, respectively. Thus, the cost of using (18), specially for large systems, is decreased considerably (note that N, + N, = N) REMARK 3.1. At the double buckling point, the matrix A, and A, are all singular. But the coefficient matrix of Step 1 is non-singular under the condition (8); furthermore, the coefficient matrix of Step 3 is also non-singular provided the condition (11) holds. The hypotheses in (8) and (11) are common in the theoretical investigations of secondary bifurcation (see [6-g]). Hence, both the two matrixes may be factored without ill-conditioning and our method is quadratically convergent.
To compute the double buckling point approximation for ($,, &, h,, pO) where +, = I,41 >
4,, = Ia42 ’
with Newton’s
method,
one has to have
a suitable
(19)
Note that the double buckling load and two corresponding buckling modes are just the eigenvalues and corresponding eigenvectors of the operator equation linearized from the trivial solution path x = 0 of (1) (see [lOI) Ax = ABx
(20)
where A and B are N x N symmetric matrixes and A positive definite and p is implicitly contained in A and B. For given CL,by using the isomorphisms Z, and I,, we can solve the two smallest buckling loads with, respectively, symmetric and antisymmetric buckling modes from (20) (see e.g. 1151). If the two buckling loads are near, then we may respectively take the isomorphisms of the two corresponding buckling modes, one of the two buckling loads, and I_Las a initial approximation to (4,, 4,, A,, po).
B. Wu I Comput.
Merhods Appl. Mech. Engrg.
189
120 (199-i) 183-lY3
4. Numerical examples If this section, we would like to present two numerical examples which will illustrate the characteristics of the algorithms developed. These examples include computation of double buckling points as well as secondary buckling points. All computations have been performed on a SUNworkstation. EXAMPLE 4.1. Coupled elastic rods as a simple system demonstrating secondary buckling [ 161. Two rigid rods are connected by a springhinge. The bottom rod is constrained by another hinge and spring, and bends in the plane y = 0 through an angle /3. The top rod bends through an angle a in the plane containing the bottom rod and the y-axis. The hinges and springs constrain the rods to buckle in mutually orthogonal planes. The top rod is subjected a vertical compressive load p. The total potential energy of the system is la2+;k,~2+p(l,+l,cosa)cos~
V(a,p,p,k,,k,.l,,fZ)=~k
(21)
where (Y, p = bending angles of top and bottom rods, 1,) I, = rod lengths, k,, k, = spring stiffnesses. Introducing non-dimensional variables y=---k, k, ’
p+, 2
*=E!? k2
one obtains non-dimensional
V=& ’
2
potential energy
V(ol,p,n,~,U)=~~12*+~~*+A(1+~coSa)cos~.
(22)
Define
s=(:, -:>
(23)
and x = (a, p)‘, then the potential energy in (22) satisfies (3a) and (3b) with fixed V. By using stationary potential principle, we derive equilibrium equations of the system ~a - AE.L sin (Ycos p = 0,
(244
p - A( 1 + p cos a!) sin p = 0 .
(24b)
From (23), we also have
I,=(1
01,
In the following becomes
z-1 = 5
0
Z,=(O l),
;
numerical
computation,
I,‘=
001
we let v = 0.25. Then the augmented
(25) system in (12)
ACL sin@4
1. 1 o
25a
_
P-hll+CLco&LI),P *-1 ;*-1
=O,
(26)
and the augmented system (7) (set E = 0) for detecting the primary buckling load related symmetric primary post-buckling branch becomes 0.25a - h/_&cX 0. (Y2-l 1 =
1
At a double buckling point of the system, (A, &, A,, F,,) = (1, 1, f , f). Using augmented
(27) systems
190
B. Wu
Table 1 Dependence
I Comput.
Methods
Appl.
Mech.
Engrg.
120 (1995)
of secondary buckling of coupled rods on p
E
IL(&)
A,,(E)
k(E)
4)
0
0.33333333 0.33336297 0.33345193 0.33407721 0.33634673 0.34601193 0.35390944 0.36437317 0.37798559 0.39554873 0.41818119 0.44747901
0.75000000 0.74993332 0.74973325 0.74833000 0.74328060 0.72251844 0.70639540 0.68610979 0.66140088 0.63203338 0.59782699 0.55868542
0.75000000 0.75013335 0.75053356 0.75334226 0.76347776 0.80575708 0.83947685 0.88336613 0.93963315 1.01168479 1.10498659 1.22882878
0.00000000 0.04000000 0.08000000 0.20000000 0.40000000 0.80000000 1.00000000 1.20000000 1.40000000 I .60000000 1.80000000 2 .oooooooo
0.04 0.08 0.20 0.40 0.80 1.00 1.20 1.40 1.60 1.80 2.00
1113-193
(26) resp. (27) t o continuate 6 resp. p with the methods in Section 3. We get Table 1 which contains the dependences of the primary buckling load A,, secondary buckling load A, and corresponding bending angle (Y, of the top rod on the length ratio p.
EXAMPLE
4.2.A simply supported
compression.
The normal displacement
d4w + h d’w ~+~w+w~=o, dx4
beam on a non-linear elastic foundation subjected to axial w(x) of the beam satisfies the following equation [7]:
(284
o
(28b) where A is proportional to the axial compression and pw + w3 is the reaction of the foundation. It is easy to check that the solution of (28) is just the stationary point of the following potential energy function under the conditions w(0) = w(n) = 0:
(29) We also see that u satisfies the Z,-symmetry U(S,w, A, I*) = VW, A, P)
(30a)
s, w(x) = w(7T-x)
W’b)
where )
x E [O, Tr] .
It is instructive to describe in detail how a discretization of (28) leads to a finite dimensional problem involving symmetry analogous to (28). To fully utilize the symmetry we choose a uniform mesh size, h = nrl(N + l), with N = 2n, II: integer, and nodes x, = ih, i = 1, 2, . , . , N. Starting from (29) and using difference-variation (evaluate the integral with trapezoid formula and use central difference quotient approximations to d*wldX* and dwldx) discretization to (28) leads to finite dimensional problem (1) where x = (w,, w2,. . . , w~)~ is the column vector of the nodes and the components of f are 7
f, =w3 - 4w, - (5 + h4p + h4w;)w,
-;
At-w,
+ 3w,) ,
f; = w,_2 - 4W,_l - 4w,+, + wj+2 + (6 + h4F + h4w;)wI 2SiaN-l,w,=w
Iv+1 = 0,
h2 - 4 A(2w, - wi_2 -
Wi+*)
)
B. Wu I Comput. Methods Appl. Mech. Engrg.
120 (1995)
183-193
191
7
+ (5 + h4p + h4w;)w,
fN = wN_* - 4~,~_,
+ % A(w,_,
- 3w,)
.
Define 0
1
s=
1 :. C1
E[W2nx2n,
(31)
01
Then the S and f above, satisfy (3a) and (4). Now we have I, =;(I,,
E,)+4,
Z,‘=(I,
E,)‘,
r,’
.*. .f
:)Wfl,
)
I, = ; (Z, = (I,
(324
-En)+!,,
Wb)
-En)’
where 1;-6
I,, I, E RnxZn ; Then the augmented
I,‘,
&=(Y
I,’
l
..* ;)ERY
(32~)
E lRZnxn .
system (12) is equivalent to the following system
(33)
By continuating p, the corresponding following augmented system
symmetric primary post-buckling load can been calculated by the
@u-JO~A, PPCW,= o
1
w;w,- 1
I .
(34)
At a double buckling point of the original continuous problem (28), (35) Let N = 80(n = 40). By using the method in Section 3, we get the dependences of the symmetric primary buckling load h,, the secondary buckling load A, and corresponding maximum deflection wraX of the beam on the stiffness of the linear part of the reaction or. as following Table 2. From Table 2, one knows that as the deviation of p from pCJ,the secondary buckling load slightly departs from the symmetric primary buckling load, yet the maximum deflection of the beam increase rapidly.
4. Conclusions An augmented system containing a newly introduced parameter has been given for the numerical calculation of secondary buckling in non-linear stability problems. The implementation of Newton’s method solving the resulted system has been developed by a partitioning method. The method directly yields secondary buckling load and corresponding secondary buckling mode without tracing primary post-buckling branch. The eigenmode associated with secondary buckling point can be used for
B. Wu I Comput. Methods Appl. Mech. Engrg.
192 Table
I20 (1995) 183-193
2
Dependence
of secondary
buckling
of beam on Jo
f
CL(E)
,$A)
A,(E)
Wyy&)
0 0.5 1.0 1.5 2.0 3.0 4.0 5.0 5.5 6.0 6.5 7.0
4.00502046 3.99885508 3 .Y803742Y 3.94962411 3.90668192 3.18469338 3.61569885 3.40159514 3.27841124 3.14499291 3.00180535 2.84938408
5.00728033 5.00111864 4.98262180 4.95185619 4.90889246 4.78684272 4.61776342 4.40355229 4.28030659 4.14682133 4.00356194 3.85106421
5.00728033 5.01037506 5.01966595 5.03517321 5.05693084 5.11940548 5.20766262 5.32255394 5.39034378 5.46526926 5.54754715 5.63742969
0.00000000 0.11108165 0.22211179 0.33303849 0.44380898 0.66466313 0.88421666 1.101Y6390 1.20998301 1.31732600 1.42390320 1.52961495
branch-switching. A consistent linearization of the governing equations leads to quadratic convergence behavior. The two numerical examples have displayed that one can follow the secondary buckling points sufficiently far away from the double buckling points so that the further computation of the points can be carried out by using a continuation procedure based on the augmented system [9]. The algorithm has been applied to coupled elastic rods and a simply supported beam in a non-linear elastic foundation. In both cases the algorithm is very successful in computing double buckling points and secondary buckling points and shows good convergence.
Acknowledgment
This work has been financially supported by the National Science Foundation of China. The author would like to thank Professor T. Kiipper and the Mathematisches Institut der Universitat zu Koln for their hospitality. This paper was finished while the author visited Koln, where he was supported by Volkswagen-Stiftung, and he would like to thank it for its support.
References
111H.
Troger, Application of bifurcation theory to the solution of nonlinear stability problems in mechanical engineering, in: T. Kiipper, H.D. Mittelmann and H. Weber, eds., Numerical Methods for Bifurcation Problems, ISNM 70 (Birkhauser Verlag, 1984). t4 L. Bauer, H.B. Keller and E.L. Reiss, Multiple eigenvalues lead to secondary bifurcation, SIAM Rev. 17 (1975) 101-122. L.J. Putnick and E.L. Reiss, Secondary states of rectangular plates, SIAM J. Appl. Math. 38 (1980) 131 B.J. Matkowsky, 38-51. Secondary buckled state of rectangular sandwich plates, Acta Mechanica Sinica 6 (1990) 141 He Luwu and Cheng Changjun, 237-247. Secondary bifurcation and imperfection sensitivity of cylindrical shell of finite length, Internat. J. [51 M.O. Oyesanya, Non-Linear Mech. 25 (1990) 405-415. SIAM J. Math. Anal. 11 (1980) 365-389. PI M. Shearer, Secondary bifurcation near a double eigenvalues, Multiple bifurcation, symmetry and secondary bifurcation, in: Nonlinear Problems of Analysis in Geometry f71 M. Potier-Ferry, and Mechanics, Research Notes in Math. 46 (Pitman, 1981) 158-167. and D. Schaeffer, Singularities and Groups in Bifurcation Theory, Vol. 1 (Springer-Verlag, 1985). PI M. Golubitsky of symmetry-breaking bifurcation points, SIAM J. Numer. Anal. 21 (1984) 191 B. Werner and A. Spence, The computation 388-399. convergent procedure for the calculation of stability points in finite [lOI P. Wriggers, W. Wagner and C. Miehe, A quadratically element analysis, Comput. Methods Appl. Mech. Engrg. 70 (1988) 329-347. [Ill C.S. Chien, Secondary bifurcation in the buckling problem, J. Comput. Appl. Math. 25 (1989) 277-287. [121 B. Wu, The compution of secondary bifurcation and its applications to buckling of structures, Ph.D. Thesis, Jilin University, Changchun, 1991.
6. Wu / Comput. Methods Appl. Mech. Engrg. [13] B. Werner,
120 (1995) 183-193
193
Regular systems for bifurcation points with underlying symmetries, in: T. Kiipper, H.D. Mittelmann and H. Weber, eds., Numerical Methods for Bifurcation Problems, ISNM 70 (Birkhauser Verlag, 1984) 562-574. [14] H.B. Keller, Lectures on Numerical Methods in Bifurcation Problems (Springer-Verlag, 1987). [15] G.H. Golub and C.F. Van Loan, Matrix Computation, 2nd edition (The Johns Hopkins University Press, Baltimore, Maryland, 1989). [16] S.J. Tavener and T. Mullin, Buckling of coupled elastic rods, Physica D 30 (1988) 382-398.