- Email: [email protected]
The finite deformation theory for beam, plate and shell Part I. The two-dimensional beam theory
Computer methods in applbd mechanics and engineering Comput. Methods Appl. Mech. Engrg. 146 (1997) 53-63
ELSEVIER
The finite deformation theory for beam, plate and shell Part I. The two-dimensional beam theory Mingrui Li Applied Mechanics,
China Agricultural
University (East Campus),
Received 15 May 1996; revised 17 September
Beijing, 100083, PR China
1996
Abstract
A finite deformation theory is proposed that can describe precisely, the nonlinear geometric behavior of a two-dimensional beam structure. Unlike the existing nonlinear theories, this theory does away with simplifications such as the assumptions of small displacement, small normal or shearing strain, small rotation, as well as small loading step. The finite element method and its formulation of this theory are deduced. Numerical examples performed by the theory proved it valid and computed results obtained by the theory are independent of loading steps.
1. Introduction The existing finite deformation theory can describe the nonlinear geometric behavior of truss, 2-D and 3-D continuum structures pretty well, but it is unsatisfactory for beam, plate and shell structures. All the nonlinear theories of the beam, plate and shell [l-13] have to adopt some kind of simplification, such as small displacement, small normal or shearing strain, small rotation or small loading step in computation. There are the Bernoulli beam and the Kirchhoff plate that neglect the shearing effect; and the Timoshenko beam and the Mindlin plate that include the shearing effect. Then there are several algorithms, such as TL, UL and CR [14-181, which are supposed to solve the nonlinear geometric problem approximately by using certain simplifications. But the different assumptions, different theories, different algorithms or even different loading steps led to different computational results. Except for those that can be compared with analytical solutions, it is very difficult to discern which results are better. Although the powerful algorithm called arc length method [19-231 can trace the load path accurately in buckling and bifurcation problems, it is an automatic incremental loading step method. Thus, when the loading step becomes larger it will violate the assumption of small deformation and may lead to absurd results. The reason for the above phenomena is the lack of an exact finite deformation theory which can describe the complete nonlinear geometric behavior of beam, plate and shell structures exactly. In this paper we present a finite deformation theory for a 2-D beam that are free of the assumptions of small displacement, small rotation, small strain and small loading step. We then derive the corresponding finite element method (FEM) by using total Lagrangian (TL) algorithm. Three numerical test problems and three examples are given to show the correctness of the theory. The general finite deformation theory will be presented in our next paper.
2. The finite deformation theory of a 2-D beam We will adopt the Timoshenko
beam theory, i.e. a plane cross section originally normal to its neutral
0045-7825/97/$17X10 0 1997 Elsevier Science S.A. All rights reserved PII SOO45-7825(96)01213-3
M. Li I Comput. Methods Appl. Mech. Engrg.
54
146 (1997) 53-63
(4 Fig. 1. Configuration
of a beam segment
at different
times.
(a) Initial configuration
at time
t,); (b) current
configuration
at time
t.
axis will always remain a plane during the deformation but because of shearing this section does not remain normal to the neutral axis. The initial configuration of a beam segment at time t, is shown in Fig. l(a). Fig. l(b) shows the current configuration at time t. The x-axis is chosen as the neutral axis in the initial configuration and the y-axis is on the cross section along its principal axis. A line segment A,B, represents an arbitrary cross section at time t, and after the deformation, it becomes the inclined line segment A, B, , its current position at time t. The basic assumption is: the length of A,,B,equals that of A ,B,,and that is the only assumption of this paper. The following corollary is obvious, along the cross section, the Poisson ratio p = 0. It follows then the shearing modulus G and the Young’s module E have the relationship
E=2G
(1)
Although I_L= 0 is adopted in all the linear and nonlinear theories of beams, the relationship between the elastic moduli has not been duly noted. For the convenience of comparing the results with those of others, E and G will both appear in the following derivations. But the relationship E = 2G will be used in the numerical computations. From the basic assumption that the line segment AB will act as a rigid body, its displacement can be taken as the combination of a rigid body translation and a rigid body rotation with a rotation angle 8. The radius vectors of A o,B,,A 1,B, are TA,,= (6 o)r
(2a)
rB0 = (x9 Y)’
(2b)
rAl
=
‘A0
+
‘A0
(W
TBl
=
rlJ0
+
so
=
rA,
+A$,
(24
where uAo= (u, u) r is the displacement of A, located on the neutral line, us0 = (U, V)',is the displacement of B,,an arbitrary point on the cross section at time to. From the basic assumption, the length of A,B, equals that of A,B,, it can be seen from Fig. l(b) that A ,B, = (-y
Substituting
sin 0, y cos 0)’
the above expressions
U=u-ysin0 v=u
-y(l
of uAO, ugo and A ,B, in the expression of rB1, we obtain (3a)
-cos0)
(3b)
Thus, the displacement of an arbitrary point can be expressed in terms of the displacement of a point on the neutral line (rigid body translation) and the rotation angle of the cross section (rigid body rotation). We will now study the deformation with Lagrange’s point of view. Taking the derivatives of Eq. (3) with respect to x and y, we obtain u,X = u,, - ye,, cos 8
(4a)
M. Li I Comput.
Methods
Appl.
Mech.
Engrg.
55
146 (1997) 53-63
U,? = -sin 8 V,, = v,, - y
(4b) sin 0
(4c)
lJy = -1 + case
(4d)
We can describe the finite deformation
of a beam with the Green-Lagrange
E,, + U,, + (U:,, + VT,)/2 = u,, + (u: + 332 E,?. = V, + u.? + UJ,
+ VJ’.,
strain,
+ y’8,,/2 - yO.xres
= yen
(54 (5b)
E,,, = 0
(5c)
Y,, = (1 + u,~) cos 8 + v,, sin 0
(6a)
r en = v,, cos 0 -
(6b)
where
(1 + u,~) sin 0
The geometric meaning of the above expressions can be described as follows. Let us define r, = ‘dsl”ds, where “ds is the length of the line segment at t,, and ‘d.s is the length of the same line segment at time t. Then r, is called the stretch, representing the relative extension of the segment. Let ‘ds = du, we will have ‘d? = (dx + dU)’ + dv’ and r<,= ((1 + u,,)~ + vTJ1”
Pa)
We can also find the tangent vector of the point A, located on the neutral line at time t, r, = (1 + u,,, v,JT
C’b)
Obviously, the length of the vector r, is re. This tangent vector describes the deformation of the neutral line. Let 4 denote the angle between r, and the x-axis, then 8 also denotes the angle between the normal at A, and y-axis as shown in Fig. l(b). We have tan 4 = v,,/(l
+ u.,)
@a)
Therefore sin 4 = v., lr,
(8b)
cos 4 = (I + U,x)lr,
(gc)
Thus, Eq. (6) can be rewritten
as
r PI =r,cos$cos6+r,sin$sin8=r,cos(+-f9)=r,cosy
Pa)
r en = rr
Pb)
sin d, cos 8 - r, cos C#Isin 8 = r, sin($ - 8 ) = rr sin -y
where y = 4 - 8 (see Fig. l(b)) d enotes the angle of the inclination of the cross section from the normal, i.e. the shearing angle. In the Bernoulli beam theory, the shearing effect is neglected, i.e. Y = (b - 0 = 0. But in the finite deformation theory it should be expressed as r t’t*= v,, cos 8 -
(1 + u,,) sin 8 = 0
(10)
We can call Eq. (10) the Bernoulli condition of finite deformation of a beam. But in the Bernoulli beam theory the relationship v., = 0 is also used. This is the result of a series of linearization assumptions: 0-tane-tan~=v,,/(l+u,,)=v,,
56
M. Li I Comput. Methods Appl. Mech. Engrg.
146 (1997) 53-63
Thus, it is no longer suitable to adopt the Bernoulli beam in the finite deformation theory. The above are deduced from our finite deformation theory. If any approximation assumptions on the strain or the rotation are used, one can obtain the existing approximate nonlinear geometric models.
3. The virtual work equation and the FEM formulation The virtual work of an elastic body at equilibrium
can be expressed as
SW(i) = SW(~)
(Ha)
where SW”’ and SW”’ denote, respectively, the internal virtual work of an elastic body at equilibrium. Let P denote the equivalent nodal force subjected to all the external nodal, surface and body forces, and let 6q denote the virtual displacements of nodes, then we have SW”’ = Q’
p
(lib)
Because the Green-Lagrange while the finite deformation SW”’ =
incremental strains are conjugated with the 2nd Piola-Kirchhoff stresses is taking place, the internal virtual work of a beam can be expressed as
(SETx S,, + SE,:, S,,) du
(114
where the Kirchhoff stresses S,,, S,, and the Green strains E,,, EXYsatisfy the constitutive (within the limit of the elastic state) S,, = EE,,
and
relationships
S,, = GExY
(lld)
Here, we only emphasize on the geometric part of the finite deformation dSij = Dijk, dE,,
theory.
If we use We)
instead of (lld) and choose a suitable constitutive relation, then it will fit for any kind of problem without changing any part of the geometric formulation. Note that we use the symbol Exy, which is twice the value of Green strain tensor E,,. Eq. (lld) may be written simply as S, = EE,
(i, j = 1,2)
(110
In this paper we will always take the initial configuration as the reference configuration to derive the FEM formulations, i.e. the TL method. For simplification, the symbol of time in the reference configuration will be omitted. To derive the variation of the Green strains we first write out the variation of Eq. (6), Sres = cos 13Su,, + sin q Su,, + r,, 68
(124
Sren = -sin 8 Su,, + cos 8 Su,, - TesF8
(12b)
Thus, we have SE,, = A,, Se
(134
SEXY= A,, Se
(13b)
e = (u,,, u,,, O,,, 0)’
(I3c)
where
4,
= (I + u,, -yf,
cos 8, u,, - yqx sin 6 Y*O,, - YL,
A,, = (-sin 8, cos 8, 0, -res)
-y%r,,)
(134 (134
M. Li I Comput. Methods Appl. Mech. Engrg.
The interpolation
146 (1997) 53-63
57
functions for u, u and 8 are defined as
u=H,q
(144
v=H,q
(14b)
0 = Hoe
(144
H is defined as Se=H&j
(154
While the types of interpolation H=[H,,H,
,...,
H; ,...,
functions for u, LJand 8 are the same, H will be H,,]
(15b)
and
(t = 1,2, . . . ) n)
(154
where n is the number of nodes of the beam element. Then Eq. (3) can be expressed as 6E,, = A xxH i5q = B,, Sq
(164
SE,, = A,,H
(16b)
Sg = BxY s4
Substituting Eq. (16) in Eq. (11) and also making use of the fact that the nodal virtual displacements can be chosen arbitrarily if the constrains are satisfied, we then have the equilibrium equation F(q) = 1” (B,T,&, + B$%J The Kirchhoff N=
(17)
internal forces can be defined as
I A L dA
(@a)
sxydA
(18b)
Q=/ M =
dy - P = 0
A
IA
Substituting
yS,, dA
Eq. (5) and Eq. (11) in Eq. (18) we obtain
N = EAu,, + EA(u: Q = GA,r,
sin y
M = - EZ0,,r, cos y A, is called the effective
+ 4/2
+ EZO32
(194
(19b) (19c)
shearing area and can be calculated in the usual way. (17) and (19) are the equilibrium equations of the FEM and the cross section, respectively. The first term of the axial force in Eq. (19a) is the same as in the linear theory. The second and third terms arise from the finite deformation of the neutral line. The fourth term is due to bending. The shearing force of Eq. (19b) shows the coupling of the extension (r,) and the shearing (sin y). The bending moment of Eq. (19~) shows the coupling of the curvature (O,X),the extension (r,) and the shearing (cos y ). We can still use a similar relationship M = --EZ*curvature*C, where C is a modification factor due to the extension and shearing, i.e.
M. Li I Comput. Methods Appl. Mech. Engrg.
58
C=
146 (1997) 53-63
cos y
re
The above relationships apply to all situations, regardless of the size of the displacements, angular rotations (within the limit of elasticity).
strains or
4. Solving the FEM equilibrium equations with the Newton’s method The FEM equilibrium equation (17) can be solved by using Newton’s method iteratively. We want to expand Eq. (17) in the form of Taylor’s series, which is F( q + Aq) = F(q) + AF + O(Aq*) = 0 Ignoring the higher-order
terms of Aq and denoting AF as
AF= K,Aq Eq. (17) becomes K, Aq = -F(q)
=P -
iu
(B,T,.S,, + B&J
du
K, is called the tangent stiffness matrix. From Eq. (17). AF can be written as AF =
(AB,‘,S,, + B,T, AS,, + AB,:, SXY+ BTy AS,,) dv
I”
Note that the operator
A only effects the nodal displacements,
(2Oc) so from Eq. (16) we obtain
AB,‘, = HTFx,H Aq
(214
AB,T, = HFx,H 6q
@lb)
where
L
=
r
1
I
0
sym _
-y cos 8 L ye,, sin 8
1 -y sin0
y*
-ye,, cos 0
-yr,n
0
Fxy= i
vm
00
-COS 8
-sin0 0 0
During the integration, I,ydA,
YQk
i,y2dA,
-r,, I
the following four terms will appear i,y-‘d&
i,y4dA
The first integration will be zero if the neutral line passes through the center of the cross section. The second is the moment of inertia. The third will also be zero if the neutral layer, i.e. the xz-plane, is the plane of symmetry of the cross section, otherwise it should be calculated. The fourth term is the fourth-order moment of the cross section. Using the constitutive law of Eq. (lld), we have AS,, = E AE,, = EB,, Aq AS_, = G AE,, = GBxy Aq Substituting
these two in the expression of AF, the tangent stiffness matrix can be written as
K, = K,, + K,, + K, + K,
(234
M. Li I Comput. Methods Appl. Mech. Engrg.
146 (1997) 53-63
5c)
where K,, = E
I”
BTxB,, dv
(23b) (23~)
K,, =
K, =
HTFx,S,,H
dv
WI
H.“F,,S,,.H
dv
P-34
I L’
(23d) and (23 e ) are usually called the geometric stiffness matrix. During the iteration process, whenever an incremental vector Aq”’ is found, the displacement vector can be immediately updated as 9
cr+l) = q”’
+ Aq”’
This is a standard Newton’s method. It is entirely unnecessary to wait for the end of the convergence of each loading step, which is essentially the modified Newton’s method and has a lower convergence rate than the standard one [13]. Using our updating method, the quadratic convergence rate can be observed whenever Newton’s attraction area is entered. If the initial value is far from the final result, the process may diverge due to the escape from Newton’s attraction area. The above theory and the FEM formulations are exactly the same for straight as well as curved beams. The difference lies in the fact that a local coordinate system can be established for a whole straight beam element, whereas for a curved beam, the local coordinate system must be established for each individual point. In practical computation, it is sufficient to establish it at each numerical integration point instead. In short, the coordination transformation matrix can be taken out of the integral for a straight beam, while it must be left inside the integral for a curved beam.
5. Numerical examples Three simple but crucial test problems and three numerical examples are given below. The purpose of the tests is to verify the validity of the proposed theory. The parameters are chosen as follows in all the test problems. A clamped beam has a rectangular cross section, b = 1, h = 6, R = 72, E = 2G = 4. The convergence condition is set to be Y, < E or
r,
where r, = IlAull/ll~Il, is t h e ratio of the Euclid norm of the incremental displacement and total displacement; rr = jlArll/ [[AP(~, is the ratio of the Euclid norm of the unbalanced force and the incremental force in a loading step. The tolerance error E:is set to be l.e-10. In such test problem the beam is discretized to 13 nodes and 25 nodes by using a 3-node, 4-node and 5-node element successively in each discretization. All together there are 6 cases for each test problem. Test 1
A straight beam has length L = 2nR. The boundary condition is set at its left end with u = v = 0, and 13specified such that the beam rotates an angle 8 around its left end. Let it rotate A8 = 7r/2 in each step. After 4 steps the beam will go back to its initial position without any deformation. All six cases passed the test. When we let A0 > 7r/2 in a loading step, it failed due to the divergence. The reason is that the initial value is too far away from Newton’s attraction area. Thus, A8 = n/2 may be called the Newtonian extreme incremental loading in this test problem.
60
M. Li I Comput. Methods Appl. Mech. Engrg.
146 (1997) 53-63
Test2 The same straight beam with length L = 27rR and boundary condition u = u = 8 = 0 at its left end has now a lumped moment applied to its right end. When M = M, = 2nEZ/L, the beam should bend into a closed circle. That means the right end will have theoretical values of u = u = 0 and 8 = 27~. During the computation, the incremental loading is taken to be AM = M,,/4in each loading step, so after 4 steps the beam will bend into a perfect circle. The results are tabulated in Tables 1 and 2. Au/L and Au/L are the ratios of the difference of the computed displacements and the theoretical displacements to the total length of the beam. AO/27r is the ratio of the difference of the computed angle and the theoretical result to the theoretical angle 21r. The results show that the 4-node and 5-node elements have very good performance, whereas the 3-node element performed unsatisfactorily. However, when the beam was discretized into a finer mesh the results improved. Comparing the 3-node and 4-node element, we find that it is better to use a higher-order interpolation function rather than a finer mesh. But when the 4-node and 5-node elements are compared, the higher-order approach is not obviously better. Therefore, the 4-node element with cubic interpolation is recommended for the finite deformation theory of beams. We also tried to take AM > n/2 in a loading step but failed. Thus, AM = MO/4 can be called the Newtonian extreme incremental loading in this test problem. When AM < MC,/4is tried, and we reached the following conclusions: (a)The final results are exactly the same, i.e. the computed results are independent of the loading steps, just as we expected. (b) Smaller loading increments result in faster convergence, however, the total loading steps will be more. Test 3 A curved beam in the shape of a half circle with the radius R is fixed at its left end, i.e. the boundary condition is set to be u = u = t? = 0 at the left end. A lumped moment M is applied to the right end to bend the beam against its curvature. The loading increment is set to be AM = M,/4, where MO = 2EZ/ R. After 2 steps the curved beam will bend into a straight line, after 4 steps it will bend into a half circle in the opposite direction, and after 6 steps it will become a closed circle. All the 25 node cases passed the test. However, only the 4-node and 5-node elements passed the 13 node mesh. The 3-node element of the 13 node case can be bent into a straight line but no further. This means the requirement of the FEM for the finite deformation is higher in both interpolation and mesh discretization. All the above computations used the 5-point Gauss numerical integration. When fewer integration points were used, the computation either diverged or produced bad results. When more integration points were tried, no improvement could be obtained at all. The following three numerical examples are all using the 4-node element only. Example 1 is from Crivelli [18] and others. They obtained the results by using 10 or more loading steps, while we use only one loading step only. And we also give the result by doubling the load and Table 1 The results of a beam with 13 nodes Nodes/elem
Total elements
Au/L
Au/L
AtIl2n
Number of iterations/step
3 nodes 4 nodes 5 nodes
6 4 3
0.06 0.056 2.3e-3
0.06 0.009 1.6e-3
0.5 0.054 2.3-4
14 17 18
Table 2 The results of a beam with 25 nodes Nodes/elem
Elements
AulL
AviL
Ael2n
Number of iterations/step
3 nodes 4 nodes 5 nodes
12 8 6
0.17 9.3e-4 2.le-3
0.11 2.3e-6 1.4e-5
0.19 9.2e-4 2.2e-4
16 18 18
M. Li
Table 3 The tip displacements
I Comput. Methods Appl. Mech. Engrg. 146 (1997) 53-63
61
of a cantilever beam with two transverse loads
Elements
U
V
0
Pl
P2
6 8 16 16
30.78 30.79 30.79 50.12
67.03 67.02 67.03 80.82
1.043 1.043 1.043 1.328
850 850 850 1700
1350 1350 1350 2700
Fig. 2. A cantilever beam subjected to two transverse loads.
using another loading step. The calculation is done by using 6, 8, 16 and more elements successively. Since no further improvements can be made by using more than 16 elements, we believe this is the exact answer according to our theory. The results are listed in Table 3. The two deformed shapes of this example are shown in Fig. 2. Example 2 is from Zienkiewicz by using plane stress elements and whose results are shown in Fig. 19.8 and Fig. 19.9 of [2]. We calculate the problem by discretizing the structure into 16 beam elements. Our results are listed in Table 4. The first 5 loading steps are away from the limit point. It is followed by 2 small loading steps approaching the limit point. Then snap through occurs. The last loading step is in post buckling state. Fig. 3 shows the initial configuration and the deformed shapes at loading steps l-5 and post buckling. Example 3 is created by the author. The structure is four 90” arcs encircled by a closed circle, all of the same radius as shown in Fig. 4. Due to the symmetry in both the vertical and the horizontal directions, only one quarter of the structure with 8 elements is analyzed. The closed circle is compressed by a load from the top and the bottom along the diameter. The results are shown in Table 5 by using 5 loading steps. But the final result is the same when only one loading step is used. It is interesting to note that though the deformation is already very large, passing the center of the circle, no buckling has occurred even with larger loads. The above examples show the efficiency and the correctness of the proposed theory, as well as the fact that the limitation of the Newtonian extreme incremental loading is not essential. Table 4 The top displacements
of a clamped-hinged
arch (16 four-node element)
Ll
V
l9
Loads (PR’/EI)
-6.747 -21.53 -43.11 -56.01 -59.80 -60.44 -60.65 -62.53 5.947
-8.609 -25.70 -58.21 -93.63 -110.1 -112.5 -113.2 -119.8 -206.1
-0.0168 -0.0560 -0.0973 -0.1123 -0.0832 -0.0703 -0.0653 -0.0018 2.783
2 4 6 8 9 9.10 9.12 9.14 (Critical) 10.0
M. Li
62
(b)
(4 Fig. 3. Large
deflection
Fig. 4. A closed circle inscribed quarter of the structure. Table
I Comput. Methods Appl. Mech. Engrg. 146 (1997) 53-63
(4 of a deep clamped-hinged
4 arcs compressed
by a vertical
arch.
(a) Initial
(b) configuration;
load. (a) Initial configuration;
(b) deformed (b) deformed
shapes. configurations
of one
5
The vertical
displacements
u
21.28 51.33 82.51 108.7 126.6
at bottom
and the horizontal
displacements
at right end of a closed
ring with 4 inscribed
arcs
Loads 15.87 28.90 33.67 31.68 26.70
8 16 24 32 40
6. Conclusion
In this paper we proposed a finite deformation theory and the FEM formulation for two dimensional beams. In the limit of the Newtonian extreme incremental loading, all the results can be obtained by using only one loading step. When more steps are used, the results will be exactly the same. This fact is important for all nonlinear geometric problems. The suggested beam element is 4 or 5 node element and the suggested number of Gauss integration points is 5. The finite deformation theory for plates, shells and 3-D beams will be presented in future papers.
References [1] K.J. Bathe, Large displacement analysis of three dimensional beam structures, Int. J. Numer. Methods Engrg. 4 (1979) 961-986. [2] K.J. Bathe, An assessment of current finite element analysis of nonlinear problems in solid mechanics, in: Symp. Numer. Solution of Partial Differential Equations (University of Maryland, 1975). [3] J. Argyris et al., On large displacement, small strain analysis of structures with rotational degrees of freedom, Comput. Methods Appl. Mech. Engrg. 14 (1978) 401-451. [4] G. Narayanan et al., An investigation of geometric nonlinear formulations of 3D beam elements. Int. J. Numer. Methods Engrg. 25 (1990) 643-662. [5] J.L. Meek, Geometric nonlinear analysis of space frames, in: Proc. Second Int. Conf. on Computer Aided Analysis and Design in Civil Engrg. (University of Roorkae, India, 1985) G20-G30. [6] K.M. Hsiao et al., Nonlinear finite element analysis of elastic frames, Comput. Struct. 26 (1987) 693-701. [7] E.N. Dvorkin et al., On the nonlinear formulation for curved Timoshenko beam elements considering large displacement/ rotation increments, Int. J. Numer. Methods Engrg. 26 (1988) 1597-1613. [8] A. Cardona, A beam finite element nonlinear theory with rotation, Int. J. Numer. Methods Engrg. ‘26 (1988) 2403-2438. [9] J.C. Simo, A finite strain beam formulation. The three dimensional dynamic problem, Part I. Comput. Methods AppI. Mech. Engrg. 58 (1986) 79-116. -101 W. Pietraszkiewicz, Geometrically nonlinear theories of thin elastic shells, Adv. Mech. 12 (1989) 51-130. 111 H. Stumpf, Buckling and post-buckling of shells for unrestricted and moderate rotations, in: E.L. Axelrad et al., eds., Flexible Shells (Springer, Berlin, 1984) 76-90.
M. Li I Comput.
Methods Appl. Mech. Engrg.
[12] Y. Basar et al., Finite rotation elements for the nonlinear
146 (1997) 53-63
63
analysis of thin shell structures, Int. .I. Solids Struct. 26 (1990) 83-97. [13] K.J. Bathe, Finite Element Procedures in Engineering Analysis (Prentice Hall, Englewood Cliffs, NJ. 1982). [14] J.L. Meek et al., Geometrically nonlinear analysis of space frames by an increment iterative technique. Comput. Methods Appl. Mech. Engrg. 47 (1984) 261-282. (151 K.M. Hsiao et al., A co-rotational procedure that handles large rotations of spatial beam structures, Comput. Methods Appl. Mech. Engrg. 27 (1987) 769-781. (161 M.A. Crisfield, Co-rotational beam element for two and three dimensional nonlinear analysis (IUTAMIIACN Symp. Vienna, Austria, 1989). (171 M.A. Crisfield, A consistent co-rotational formulation for nonlinear three dimensional beam elements, Comput. Methods Appl. Mech. Engrg. 81 (1990) 131-150. [18] L.A. Crivelli et al., A three dimensional nonlinear Timoshenko beam based on the core congruential formulation. Int. J. Numer. Methods Engrg. 21 (1993) 3647-3673. [19] E. Ricks, The application of Newton’s method to the problem of elastic stability, J. Appl. Mech. 39 (1979) 1060-1066. [ZO] G.A. Wempner. Discrete approximations related to the nonlinear theories of solids, Int. J. Solids Struct. 7 (1971) 1581-1599. [21] M.A. Crisfeld, A fast incremental iterative solution procedure that handles snap-through, Comput. Struct. 13 (1981) 55-62. [22] P.X. Bellini et al., An improved automatic incremental algorithm for the efficient solution of nonlinear finite element equations, Comput. Struct. 26 (1987) 99-110. (231 Mingrui Li, A new arc length method for tracing the complete load path, in: Recent Developments in Finite Element Analysis (Palo Alto, CA, 1994) 171-180.
ELSEVIER
The finite deformation theory for beam, plate and shell Part I. The two-dimensional beam theory Mingrui Li Applied Mechanics,
China Agricultural
University (East Campus),
Received 15 May 1996; revised 17 September
Beijing, 100083, PR China
1996
Abstract
A finite deformation theory is proposed that can describe precisely, the nonlinear geometric behavior of a two-dimensional beam structure. Unlike the existing nonlinear theories, this theory does away with simplifications such as the assumptions of small displacement, small normal or shearing strain, small rotation, as well as small loading step. The finite element method and its formulation of this theory are deduced. Numerical examples performed by the theory proved it valid and computed results obtained by the theory are independent of loading steps.
1. Introduction The existing finite deformation theory can describe the nonlinear geometric behavior of truss, 2-D and 3-D continuum structures pretty well, but it is unsatisfactory for beam, plate and shell structures. All the nonlinear theories of the beam, plate and shell [l-13] have to adopt some kind of simplification, such as small displacement, small normal or shearing strain, small rotation or small loading step in computation. There are the Bernoulli beam and the Kirchhoff plate that neglect the shearing effect; and the Timoshenko beam and the Mindlin plate that include the shearing effect. Then there are several algorithms, such as TL, UL and CR [14-181, which are supposed to solve the nonlinear geometric problem approximately by using certain simplifications. But the different assumptions, different theories, different algorithms or even different loading steps led to different computational results. Except for those that can be compared with analytical solutions, it is very difficult to discern which results are better. Although the powerful algorithm called arc length method [19-231 can trace the load path accurately in buckling and bifurcation problems, it is an automatic incremental loading step method. Thus, when the loading step becomes larger it will violate the assumption of small deformation and may lead to absurd results. The reason for the above phenomena is the lack of an exact finite deformation theory which can describe the complete nonlinear geometric behavior of beam, plate and shell structures exactly. In this paper we present a finite deformation theory for a 2-D beam that are free of the assumptions of small displacement, small rotation, small strain and small loading step. We then derive the corresponding finite element method (FEM) by using total Lagrangian (TL) algorithm. Three numerical test problems and three examples are given to show the correctness of the theory. The general finite deformation theory will be presented in our next paper.
2. The finite deformation theory of a 2-D beam We will adopt the Timoshenko
beam theory, i.e. a plane cross section originally normal to its neutral
0045-7825/97/$17X10 0 1997 Elsevier Science S.A. All rights reserved PII SOO45-7825(96)01213-3
M. Li I Comput. Methods Appl. Mech. Engrg.
54
146 (1997) 53-63
(4 Fig. 1. Configuration
of a beam segment
at different
times.
(a) Initial configuration
at time
t,); (b) current
configuration
at time
t.
axis will always remain a plane during the deformation but because of shearing this section does not remain normal to the neutral axis. The initial configuration of a beam segment at time t, is shown in Fig. l(a). Fig. l(b) shows the current configuration at time t. The x-axis is chosen as the neutral axis in the initial configuration and the y-axis is on the cross section along its principal axis. A line segment A,B, represents an arbitrary cross section at time t, and after the deformation, it becomes the inclined line segment A, B, , its current position at time t. The basic assumption is: the length of A,,B,equals that of A ,B,,and that is the only assumption of this paper. The following corollary is obvious, along the cross section, the Poisson ratio p = 0. It follows then the shearing modulus G and the Young’s module E have the relationship
E=2G
(1)
Although I_L= 0 is adopted in all the linear and nonlinear theories of beams, the relationship between the elastic moduli has not been duly noted. For the convenience of comparing the results with those of others, E and G will both appear in the following derivations. But the relationship E = 2G will be used in the numerical computations. From the basic assumption that the line segment AB will act as a rigid body, its displacement can be taken as the combination of a rigid body translation and a rigid body rotation with a rotation angle 8. The radius vectors of A o,B,,A 1,B, are TA,,= (6 o)r
(2a)
rB0 = (x9 Y)’
(2b)
rAl
=
‘A0
+
‘A0
(W
TBl
=
rlJ0
+
so
=
rA,
+A$,
(24
where uAo= (u, u) r is the displacement of A, located on the neutral line, us0 = (U, V)',is the displacement of B,,an arbitrary point on the cross section at time to. From the basic assumption, the length of A,B, equals that of A,B,, it can be seen from Fig. l(b) that A ,B, = (-y
Substituting
sin 0, y cos 0)’
the above expressions
U=u-ysin0 v=u
-y(l
of uAO, ugo and A ,B, in the expression of rB1, we obtain (3a)
-cos0)
(3b)
Thus, the displacement of an arbitrary point can be expressed in terms of the displacement of a point on the neutral line (rigid body translation) and the rotation angle of the cross section (rigid body rotation). We will now study the deformation with Lagrange’s point of view. Taking the derivatives of Eq. (3) with respect to x and y, we obtain u,X = u,, - ye,, cos 8
(4a)
M. Li I Comput.
Methods
Appl.
Mech.
Engrg.
55
146 (1997) 53-63
U,? = -sin 8 V,, = v,, - y
(4b) sin 0
(4c)
lJy = -1 + case
(4d)
We can describe the finite deformation
of a beam with the Green-Lagrange
E,, + U,, + (U:,, + VT,)/2 = u,, + (u: + 332 E,?. = V, + u.? + UJ,
+ VJ’.,
strain,
+ y’8,,/2 - yO.xres
= yen
(54 (5b)
E,,, = 0
(5c)
Y,, = (1 + u,~) cos 8 + v,, sin 0
(6a)
r en = v,, cos 0 -
(6b)
where
(1 + u,~) sin 0
The geometric meaning of the above expressions can be described as follows. Let us define r, = ‘dsl”ds, where “ds is the length of the line segment at t,, and ‘d.s is the length of the same line segment at time t. Then r, is called the stretch, representing the relative extension of the segment. Let ‘ds = du, we will have ‘d? = (dx + dU)’ + dv’ and r<,= ((1 + u,,)~ + vTJ1”
Pa)
We can also find the tangent vector of the point A, located on the neutral line at time t, r, = (1 + u,,, v,JT
C’b)
Obviously, the length of the vector r, is re. This tangent vector describes the deformation of the neutral line. Let 4 denote the angle between r, and the x-axis, then 8 also denotes the angle between the normal at A, and y-axis as shown in Fig. l(b). We have tan 4 = v,,/(l
+ u.,)
@a)
Therefore sin 4 = v., lr,
(8b)
cos 4 = (I + U,x)lr,
(gc)
Thus, Eq. (6) can be rewritten
as
r PI =r,cos$cos6+r,sin$sin8=r,cos(+-f9)=r,cosy
Pa)
r en = rr
Pb)
sin d, cos 8 - r, cos C#Isin 8 = r, sin($ - 8 ) = rr sin -y
where y = 4 - 8 (see Fig. l(b)) d enotes the angle of the inclination of the cross section from the normal, i.e. the shearing angle. In the Bernoulli beam theory, the shearing effect is neglected, i.e. Y = (b - 0 = 0. But in the finite deformation theory it should be expressed as r t’t*= v,, cos 8 -
(1 + u,,) sin 8 = 0
(10)
We can call Eq. (10) the Bernoulli condition of finite deformation of a beam. But in the Bernoulli beam theory the relationship v., = 0 is also used. This is the result of a series of linearization assumptions: 0-tane-tan~=v,,/(l+u,,)=v,,
56
M. Li I Comput. Methods Appl. Mech. Engrg.
146 (1997) 53-63
Thus, it is no longer suitable to adopt the Bernoulli beam in the finite deformation theory. The above are deduced from our finite deformation theory. If any approximation assumptions on the strain or the rotation are used, one can obtain the existing approximate nonlinear geometric models.
3. The virtual work equation and the FEM formulation The virtual work of an elastic body at equilibrium
can be expressed as
SW(i) = SW(~)
(Ha)
where SW”’ and SW”’ denote, respectively, the internal virtual work of an elastic body at equilibrium. Let P denote the equivalent nodal force subjected to all the external nodal, surface and body forces, and let 6q denote the virtual displacements of nodes, then we have SW”’ = Q’
p
(lib)
Because the Green-Lagrange while the finite deformation SW”’ =
incremental strains are conjugated with the 2nd Piola-Kirchhoff stresses is taking place, the internal virtual work of a beam can be expressed as
(SETx S,, + SE,:, S,,) du
(114
where the Kirchhoff stresses S,,, S,, and the Green strains E,,, EXYsatisfy the constitutive (within the limit of the elastic state) S,, = EE,,
and
relationships
S,, = GExY
(lld)
Here, we only emphasize on the geometric part of the finite deformation dSij = Dijk, dE,,
theory.
If we use We)
instead of (lld) and choose a suitable constitutive relation, then it will fit for any kind of problem without changing any part of the geometric formulation. Note that we use the symbol Exy, which is twice the value of Green strain tensor E,,. Eq. (lld) may be written simply as S, = EE,
(i, j = 1,2)
(110
In this paper we will always take the initial configuration as the reference configuration to derive the FEM formulations, i.e. the TL method. For simplification, the symbol of time in the reference configuration will be omitted. To derive the variation of the Green strains we first write out the variation of Eq. (6), Sres = cos 13Su,, + sin q Su,, + r,, 68
(124
Sren = -sin 8 Su,, + cos 8 Su,, - TesF8
(12b)
Thus, we have SE,, = A,, Se
(134
SEXY= A,, Se
(13b)
e = (u,,, u,,, O,,, 0)’
(I3c)
where
4,
= (I + u,, -yf,
cos 8, u,, - yqx sin 6 Y*O,, - YL,
A,, = (-sin 8, cos 8, 0, -res)
-y%r,,)
(134 (134
M. Li I Comput. Methods Appl. Mech. Engrg.
The interpolation
146 (1997) 53-63
57
functions for u, u and 8 are defined as
u=H,q
(144
v=H,q
(14b)
0 = Hoe
(144
H is defined as Se=H&j
(154
While the types of interpolation H=[H,,H,
,...,
H; ,...,
functions for u, LJand 8 are the same, H will be H,,]
(15b)
and
(t = 1,2, . . . ) n)
(154
where n is the number of nodes of the beam element. Then Eq. (3) can be expressed as 6E,, = A xxH i5q = B,, Sq
(164
SE,, = A,,H
(16b)
Sg = BxY s4
Substituting Eq. (16) in Eq. (11) and also making use of the fact that the nodal virtual displacements can be chosen arbitrarily if the constrains are satisfied, we then have the equilibrium equation F(q) = 1” (B,T,&, + B$%J The Kirchhoff N=
(17)
internal forces can be defined as
I A L dA
(@a)
sxydA
(18b)
Q=/ M =
dy - P = 0
A
IA
Substituting
yS,, dA
Eq. (5) and Eq. (11) in Eq. (18) we obtain
N = EAu,, + EA(u: Q = GA,r,
sin y
M = - EZ0,,r, cos y A, is called the effective
+ 4/2
+ EZO32
(194
(19b) (19c)
shearing area and can be calculated in the usual way. (17) and (19) are the equilibrium equations of the FEM and the cross section, respectively. The first term of the axial force in Eq. (19a) is the same as in the linear theory. The second and third terms arise from the finite deformation of the neutral line. The fourth term is due to bending. The shearing force of Eq. (19b) shows the coupling of the extension (r,) and the shearing (sin y). The bending moment of Eq. (19~) shows the coupling of the curvature (O,X),the extension (r,) and the shearing (cos y ). We can still use a similar relationship M = --EZ*curvature*C, where C is a modification factor due to the extension and shearing, i.e.
M. Li I Comput. Methods Appl. Mech. Engrg.
58
C=
146 (1997) 53-63
cos y
re
The above relationships apply to all situations, regardless of the size of the displacements, angular rotations (within the limit of elasticity).
strains or
4. Solving the FEM equilibrium equations with the Newton’s method The FEM equilibrium equation (17) can be solved by using Newton’s method iteratively. We want to expand Eq. (17) in the form of Taylor’s series, which is F( q + Aq) = F(q) + AF + O(Aq*) = 0 Ignoring the higher-order
terms of Aq and denoting AF as
AF= K,Aq Eq. (17) becomes K, Aq = -F(q)
=P -
iu
(B,T,.S,, + B&J
du
K, is called the tangent stiffness matrix. From Eq. (17). AF can be written as AF =
(AB,‘,S,, + B,T, AS,, + AB,:, SXY+ BTy AS,,) dv
I”
Note that the operator
A only effects the nodal displacements,
(2Oc) so from Eq. (16) we obtain
AB,‘, = HTFx,H Aq
(214
AB,T, = HFx,H 6q
@lb)
where
L
=
r
1
I
0
sym _
-y cos 8 L ye,, sin 8
1 -y sin0
y*
-ye,, cos 0
-yr,n
0
Fxy= i
vm
00
-COS 8
-sin0 0 0
During the integration, I,ydA,
YQk
i,y2dA,
-r,, I
the following four terms will appear i,y-‘d&
i,y4dA
The first integration will be zero if the neutral line passes through the center of the cross section. The second is the moment of inertia. The third will also be zero if the neutral layer, i.e. the xz-plane, is the plane of symmetry of the cross section, otherwise it should be calculated. The fourth term is the fourth-order moment of the cross section. Using the constitutive law of Eq. (lld), we have AS,, = E AE,, = EB,, Aq AS_, = G AE,, = GBxy Aq Substituting
these two in the expression of AF, the tangent stiffness matrix can be written as
K, = K,, + K,, + K, + K,
(234
M. Li I Comput. Methods Appl. Mech. Engrg.
146 (1997) 53-63
5c)
where K,, = E
I”
BTxB,, dv
(23b) (23~)
K,, =
K, =
HTFx,S,,H
dv
WI
H.“F,,S,,.H
dv
P-34
I L’
(23d) and (23 e ) are usually called the geometric stiffness matrix. During the iteration process, whenever an incremental vector Aq”’ is found, the displacement vector can be immediately updated as 9
cr+l) = q”’
+ Aq”’
This is a standard Newton’s method. It is entirely unnecessary to wait for the end of the convergence of each loading step, which is essentially the modified Newton’s method and has a lower convergence rate than the standard one [13]. Using our updating method, the quadratic convergence rate can be observed whenever Newton’s attraction area is entered. If the initial value is far from the final result, the process may diverge due to the escape from Newton’s attraction area. The above theory and the FEM formulations are exactly the same for straight as well as curved beams. The difference lies in the fact that a local coordinate system can be established for a whole straight beam element, whereas for a curved beam, the local coordinate system must be established for each individual point. In practical computation, it is sufficient to establish it at each numerical integration point instead. In short, the coordination transformation matrix can be taken out of the integral for a straight beam, while it must be left inside the integral for a curved beam.
5. Numerical examples Three simple but crucial test problems and three numerical examples are given below. The purpose of the tests is to verify the validity of the proposed theory. The parameters are chosen as follows in all the test problems. A clamped beam has a rectangular cross section, b = 1, h = 6, R = 72, E = 2G = 4. The convergence condition is set to be Y, < E or
r,
where r, = IlAull/ll~Il, is t h e ratio of the Euclid norm of the incremental displacement and total displacement; rr = jlArll/ [[AP(~, is the ratio of the Euclid norm of the unbalanced force and the incremental force in a loading step. The tolerance error E:is set to be l.e-10. In such test problem the beam is discretized to 13 nodes and 25 nodes by using a 3-node, 4-node and 5-node element successively in each discretization. All together there are 6 cases for each test problem. Test 1
A straight beam has length L = 2nR. The boundary condition is set at its left end with u = v = 0, and 13specified such that the beam rotates an angle 8 around its left end. Let it rotate A8 = 7r/2 in each step. After 4 steps the beam will go back to its initial position without any deformation. All six cases passed the test. When we let A0 > 7r/2 in a loading step, it failed due to the divergence. The reason is that the initial value is too far away from Newton’s attraction area. Thus, A8 = n/2 may be called the Newtonian extreme incremental loading in this test problem.
60
M. Li I Comput. Methods Appl. Mech. Engrg.
146 (1997) 53-63
Test2 The same straight beam with length L = 27rR and boundary condition u = u = 8 = 0 at its left end has now a lumped moment applied to its right end. When M = M, = 2nEZ/L, the beam should bend into a closed circle. That means the right end will have theoretical values of u = u = 0 and 8 = 27~. During the computation, the incremental loading is taken to be AM = M,,/4in each loading step, so after 4 steps the beam will bend into a perfect circle. The results are tabulated in Tables 1 and 2. Au/L and Au/L are the ratios of the difference of the computed displacements and the theoretical displacements to the total length of the beam. AO/27r is the ratio of the difference of the computed angle and the theoretical result to the theoretical angle 21r. The results show that the 4-node and 5-node elements have very good performance, whereas the 3-node element performed unsatisfactorily. However, when the beam was discretized into a finer mesh the results improved. Comparing the 3-node and 4-node element, we find that it is better to use a higher-order interpolation function rather than a finer mesh. But when the 4-node and 5-node elements are compared, the higher-order approach is not obviously better. Therefore, the 4-node element with cubic interpolation is recommended for the finite deformation theory of beams. We also tried to take AM > n/2 in a loading step but failed. Thus, AM = MO/4 can be called the Newtonian extreme incremental loading in this test problem. When AM < MC,/4is tried, and we reached the following conclusions: (a)The final results are exactly the same, i.e. the computed results are independent of the loading steps, just as we expected. (b) Smaller loading increments result in faster convergence, however, the total loading steps will be more. Test 3 A curved beam in the shape of a half circle with the radius R is fixed at its left end, i.e. the boundary condition is set to be u = u = t? = 0 at the left end. A lumped moment M is applied to the right end to bend the beam against its curvature. The loading increment is set to be AM = M,/4, where MO = 2EZ/ R. After 2 steps the curved beam will bend into a straight line, after 4 steps it will bend into a half circle in the opposite direction, and after 6 steps it will become a closed circle. All the 25 node cases passed the test. However, only the 4-node and 5-node elements passed the 13 node mesh. The 3-node element of the 13 node case can be bent into a straight line but no further. This means the requirement of the FEM for the finite deformation is higher in both interpolation and mesh discretization. All the above computations used the 5-point Gauss numerical integration. When fewer integration points were used, the computation either diverged or produced bad results. When more integration points were tried, no improvement could be obtained at all. The following three numerical examples are all using the 4-node element only. Example 1 is from Crivelli [18] and others. They obtained the results by using 10 or more loading steps, while we use only one loading step only. And we also give the result by doubling the load and Table 1 The results of a beam with 13 nodes Nodes/elem
Total elements
Au/L
Au/L
AtIl2n
Number of iterations/step
3 nodes 4 nodes 5 nodes
6 4 3
0.06 0.056 2.3e-3
0.06 0.009 1.6e-3
0.5 0.054 2.3-4
14 17 18
Table 2 The results of a beam with 25 nodes Nodes/elem
Elements
AulL
AviL
Ael2n
Number of iterations/step
3 nodes 4 nodes 5 nodes
12 8 6
0.17 9.3e-4 2.le-3
0.11 2.3e-6 1.4e-5
0.19 9.2e-4 2.2e-4
16 18 18
M. Li
Table 3 The tip displacements
I Comput. Methods Appl. Mech. Engrg. 146 (1997) 53-63
61
of a cantilever beam with two transverse loads
Elements
U
V
0
Pl
P2
6 8 16 16
30.78 30.79 30.79 50.12
67.03 67.02 67.03 80.82
1.043 1.043 1.043 1.328
850 850 850 1700
1350 1350 1350 2700
Fig. 2. A cantilever beam subjected to two transverse loads.
using another loading step. The calculation is done by using 6, 8, 16 and more elements successively. Since no further improvements can be made by using more than 16 elements, we believe this is the exact answer according to our theory. The results are listed in Table 3. The two deformed shapes of this example are shown in Fig. 2. Example 2 is from Zienkiewicz by using plane stress elements and whose results are shown in Fig. 19.8 and Fig. 19.9 of [2]. We calculate the problem by discretizing the structure into 16 beam elements. Our results are listed in Table 4. The first 5 loading steps are away from the limit point. It is followed by 2 small loading steps approaching the limit point. Then snap through occurs. The last loading step is in post buckling state. Fig. 3 shows the initial configuration and the deformed shapes at loading steps l-5 and post buckling. Example 3 is created by the author. The structure is four 90” arcs encircled by a closed circle, all of the same radius as shown in Fig. 4. Due to the symmetry in both the vertical and the horizontal directions, only one quarter of the structure with 8 elements is analyzed. The closed circle is compressed by a load from the top and the bottom along the diameter. The results are shown in Table 5 by using 5 loading steps. But the final result is the same when only one loading step is used. It is interesting to note that though the deformation is already very large, passing the center of the circle, no buckling has occurred even with larger loads. The above examples show the efficiency and the correctness of the proposed theory, as well as the fact that the limitation of the Newtonian extreme incremental loading is not essential. Table 4 The top displacements
of a clamped-hinged
arch (16 four-node element)
Ll
V
l9
Loads (PR’/EI)
-6.747 -21.53 -43.11 -56.01 -59.80 -60.44 -60.65 -62.53 5.947
-8.609 -25.70 -58.21 -93.63 -110.1 -112.5 -113.2 -119.8 -206.1
-0.0168 -0.0560 -0.0973 -0.1123 -0.0832 -0.0703 -0.0653 -0.0018 2.783
2 4 6 8 9 9.10 9.12 9.14 (Critical) 10.0
M. Li
62
(b)
(4 Fig. 3. Large
deflection
Fig. 4. A closed circle inscribed quarter of the structure. Table
I Comput. Methods Appl. Mech. Engrg. 146 (1997) 53-63
(4 of a deep clamped-hinged
4 arcs compressed
by a vertical
arch.
(a) Initial
(b) configuration;
load. (a) Initial configuration;
(b) deformed (b) deformed
shapes. configurations
of one
5
The vertical
displacements
u
21.28 51.33 82.51 108.7 126.6
at bottom
and the horizontal
displacements
at right end of a closed
ring with 4 inscribed
arcs
Loads 15.87 28.90 33.67 31.68 26.70
8 16 24 32 40
6. Conclusion
In this paper we proposed a finite deformation theory and the FEM formulation for two dimensional beams. In the limit of the Newtonian extreme incremental loading, all the results can be obtained by using only one loading step. When more steps are used, the results will be exactly the same. This fact is important for all nonlinear geometric problems. The suggested beam element is 4 or 5 node element and the suggested number of Gauss integration points is 5. The finite deformation theory for plates, shells and 3-D beams will be presented in future papers.
References [1] K.J. Bathe, Large displacement analysis of three dimensional beam structures, Int. J. Numer. Methods Engrg. 4 (1979) 961-986. [2] K.J. Bathe, An assessment of current finite element analysis of nonlinear problems in solid mechanics, in: Symp. Numer. Solution of Partial Differential Equations (University of Maryland, 1975). [3] J. Argyris et al., On large displacement, small strain analysis of structures with rotational degrees of freedom, Comput. Methods Appl. Mech. Engrg. 14 (1978) 401-451. [4] G. Narayanan et al., An investigation of geometric nonlinear formulations of 3D beam elements. Int. J. Numer. Methods Engrg. 25 (1990) 643-662. [5] J.L. Meek, Geometric nonlinear analysis of space frames, in: Proc. Second Int. Conf. on Computer Aided Analysis and Design in Civil Engrg. (University of Roorkae, India, 1985) G20-G30. [6] K.M. Hsiao et al., Nonlinear finite element analysis of elastic frames, Comput. Struct. 26 (1987) 693-701. [7] E.N. Dvorkin et al., On the nonlinear formulation for curved Timoshenko beam elements considering large displacement/ rotation increments, Int. J. Numer. Methods Engrg. 26 (1988) 1597-1613. [8] A. Cardona, A beam finite element nonlinear theory with rotation, Int. J. Numer. Methods Engrg. ‘26 (1988) 2403-2438. [9] J.C. Simo, A finite strain beam formulation. The three dimensional dynamic problem, Part I. Comput. Methods AppI. Mech. Engrg. 58 (1986) 79-116. -101 W. Pietraszkiewicz, Geometrically nonlinear theories of thin elastic shells, Adv. Mech. 12 (1989) 51-130. 111 H. Stumpf, Buckling and post-buckling of shells for unrestricted and moderate rotations, in: E.L. Axelrad et al., eds., Flexible Shells (Springer, Berlin, 1984) 76-90.
M. Li I Comput.
Methods Appl. Mech. Engrg.
[12] Y. Basar et al., Finite rotation elements for the nonlinear
146 (1997) 53-63
63
analysis of thin shell structures, Int. .I. Solids Struct. 26 (1990) 83-97. [13] K.J. Bathe, Finite Element Procedures in Engineering Analysis (Prentice Hall, Englewood Cliffs, NJ. 1982). [14] J.L. Meek et al., Geometrically nonlinear analysis of space frames by an increment iterative technique. Comput. Methods Appl. Mech. Engrg. 47 (1984) 261-282. (151 K.M. Hsiao et al., A co-rotational procedure that handles large rotations of spatial beam structures, Comput. Methods Appl. Mech. Engrg. 27 (1987) 769-781. (161 M.A. Crisfield, Co-rotational beam element for two and three dimensional nonlinear analysis (IUTAMIIACN Symp. Vienna, Austria, 1989). (171 M.A. Crisfield, A consistent co-rotational formulation for nonlinear three dimensional beam elements, Comput. Methods Appl. Mech. Engrg. 81 (1990) 131-150. [18] L.A. Crivelli et al., A three dimensional nonlinear Timoshenko beam based on the core congruential formulation. Int. J. Numer. Methods Engrg. 21 (1993) 3647-3673. [19] E. Ricks, The application of Newton’s method to the problem of elastic stability, J. Appl. Mech. 39 (1979) 1060-1066. [ZO] G.A. Wempner. Discrete approximations related to the nonlinear theories of solids, Int. J. Solids Struct. 7 (1971) 1581-1599. [21] M.A. Crisfeld, A fast incremental iterative solution procedure that handles snap-through, Comput. Struct. 13 (1981) 55-62. [22] P.X. Bellini et al., An improved automatic incremental algorithm for the efficient solution of nonlinear finite element equations, Comput. Struct. 26 (1987) 99-110. (231 Mingrui Li, A new arc length method for tracing the complete load path, in: Recent Developments in Finite Element Analysis (Palo Alto, CA, 1994) 171-180.