- Email: [email protected]
Elastoplastic buckling analysis of thin-walled structures
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.1 (1-15)
Aerospace Science and Technology ••• (••••) •••–•••
1
Contents lists available at ScienceDirect
2
67 68
3
Aerospace Science and Technology
4 5
69 70 71
6
72
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13 14
Elastoplastic buckling analysis of thin-walled structures
80 81
16
82
a r t i c l e
i n f o
83
a b s t r a c t
18 19 20 21
84
Article history: Received 21 July 2014 Accepted 4 March 2015 Available online xxxx
22 23 24 25 26 27
78 79
Eugenio Ruocco 1
15 17
77
Keywords: Elastic–plastic buckling Kantorovich technique Buckling of stiffened plates Ramberg–Osgood law Buckling paradox
28 29
The study is concerned with the elastoplastic buckling of thin-walled beams and stiffened plates, subjected to in-plane, uniformly distributed, uniaxial and biaxial load. The ruling differential equations have been solved analytically by using the Kantorovich technique and the obtained displacement field has been employed in a general procedure that, by using the framework derived by the finite element method, is able to analyze the elastoplastic buckling behavior of prismatic beams and stiffened plates with arbitrary cross-section. The inelastic effect is modeled through a stress–strain law of the Ramberg– Osgood type, and both the incremental deformation theory and the J 2 flow theory are here considered. The reliability of the numerical procedure is illustrated for rectangular plates, and the contradicting results obtained by using the two plastic theories are discussed in detail. Finally, the performance of the method is illustrated through the analysis of a C-section and five different closed section columns. © 2015 Published by Elsevier Masson SAS.
85 86 87 88 89 90 91 92 93 94 95
30
96
31
97
32
98
33 34
99
1. Introduction
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
Due to their relevance in the design of thin structures, the critical behavior of plates, both flat and stiffened, has been extensively studied in the past years, as documented by the number of text [1,2] and papers [3,4] on the issue. More recently, a number of authors have also investigated the effect of elastoplastic constitutive models on the prediction of the buckling load of thin structures, generally by employing models based on both the incremental J 2 deformation theory (DT) and the J 2 flow theory (FT). Starting from the pioneering works of – among others – Handelmann and Prager [5], Bijlaard [6] and Bleich [7] in the 1950s, different aspects of buckling behavior in elastoplastic range have been analyzed, and collected in papers dedicated to the analytical investigation of buckling [8] and postbuckling [9] of simple structures. Some of these works examine the causes of the “plastic buckling paradox”, as it was first revealed in a series of theoretical investigations on the buckling of a simply supported flat plate of infinite length, through both the FT [5] and the DT [6]. A comparison of the results shows that “the buckling load obtained from the “more rational” FT model always far exceeds the experimental observation, whereas the results obtained adopting the “simplest”, total strain DT, predict buckling loads lower than FT model and in any case in closest agreement with those obtained experimentally. Moreover, the differences between the two models increase as the level of plasticity
60 61 62 63 64 65 66
1
E-mail address: [email protected]. Tel.: +39 081 5010266.
http://dx.doi.org/10.1016/j.ast.2015.03.001 1270-9638/© 2015 Published by Elsevier Masson SAS.
increases, giving in some cases results with different order of magnitude” [10]. As reported in [8], such aspects drastically affect the behavior and the load carrying capacity of flat plates subject to uniaxial or biaxial in-plane load, regardless of the boundary conditions. For instance, Durban and Zuckerman found, for specified compression/tension ratios, an optimal loading path for the DT that has no correspondence in FT. The works of Wang et al., which analyze the elasto-plastic buckling of thin [11] and thick [12] flat plates by differential quadrature method, confirm the results obtained in [8] and show that the paradox also involves flat plates of different geometry, such as triangular or elliptical, with different boundary conditions. A comparison of incremental DT and flow rule has been recently proposed in [13], where a new algorithm, based on a Generalized Differential Quadrature (GDQ) discretization technique to solve the out-of-plane plate equation, has been presented. The aim of this work is to present a new analytical formulation to define the critical behavior of both flat and stiffened plates, as well as of thin walled columns with open or closed cross sections, and characterized by non-linear elastoplastic behavior. To this end a systematic procedure based on a finite element algorithm capable of accounting for both local and global phenomena has been developed. As they are analytical, the obtained results do not depend on the discretization adopted, and reliable solutions can be achieved with the minimum number of elements required to represent the structural geometry, with the double effect of stabilizing the numerical procedure and accelerating the computing time required for the eigenvalue and eigenvector analysis.
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
2
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.2 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
Fig. 1. Geometry and load condition of a rectangular plate.
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
78
For validation purposes, the proposed procedure has been applied to flat plates, and the results compared to the analytical and/or numerical solutions available in literature. In such a preliminary analysis particular attention is paid to the evaluation of the instantaneous elastic moduli, and to the differences between them deriving from the two models adopted. The analysis indicates that the (correct) use of the constant value of Lamè modulus G in FT shows, as consequence, an irregular behavior of the corresponding generalized elastoplastic modulus which, when the ratio between the in-plane loads is greater than a critical value, in plastic range increases its stiffness, in contrast to what happens to the other modulus. Having adopted the secant modulus G s , the corresponding modulus calculated in the DT does not present the same anomaly. Observing that just in correspondence of such critical ratio the buckling path corresponding to the two theories begins to diverge, it is possible to assert that the irregular behavior of the elastic modulus related to the G-coefficient is a cause, or at least a concomitant cause, of the plastic buckling paradox. Another aspect considered in this work is the adequacy of the von Kármán hypothesis in the negligible contribution of secondorder strain terms related to in-plane displacement. In a previous work devoted to the influence of such terms on the elastic buckling of shell and both flat and stiffened plates [14], the authors demonstrate how, when the critical mode involves comparable in-plane and out-of-plane displacement, the omission of the non-linear strain terms related to in-plane displacements can considerably overestimate the critical load. If such cases never happen for flat plates, for which the von Kármán simplification is thus correct, it is quite common for stiffened plates and prismatic beams, when buckling mode involves torsional, flexo-torsional or global flexural buckling displacement. In the present work is – to the author’s knowledge for the first time – the influence of such nonlinear terms on the elasto-plastic buckling of stiffened plates is presented and discussed. In order to capture the inelastic effect, the constitutive behavior is here modeled through a stress–strain law of the Ramberg– Osgood type that, with the basic equations ruling the elastoplastic buckling of Kirchhoff plates subjected to uniaxial or biaxial compression and under both the FT and DT hypothesis, are introduced in the following section. Section 3 reports an analytical solution based on the Kantorovich technique, and its generalization in a FEM-like procedure for the buckling analysis of structures ascribable to plates rigidly connected together along their edges. The same section includes the numerical procedure necessary to solve the eigenproblem associated to the corresponding elastoplastic stiffness matrix. In Section 4 parametric analyses related to a single plate with different boundary and load conditions, a simply supported open C-section and five different thin-walled closed section beams are reported and discussed. The last part of the paper contains some conclusive considerations.
2. Governing equations
79
Consider a rectangular isotropic plate of dimension (a, b) and uniform thickness h, schematically represented, together with the local reference system adopted in this work, in Fig. 1. The plate is subjected to biaxial in-plane compressive load nx = −ξ Ph and n y = −ηPh, with (ξ, η) fixed ratio parameters, such that (ξ = η = 1) describes equibiaxial compression, (ξ = 1, η = 0) uniaxial compression in x direction, and so on. In plane stress conditions the constitutive behavior can be defined as follows:
81
⎡
⎤
80
⎡
⎤⎡
⎤
σ˙ x αxx αxy αxz ε˙ x ⎣ σ˙ y ⎦ = E ⎣ αxy α y y α yz ⎦ ⎣ ε˙ y ⎦ τ˙xy αxz α yz αzz γ˙xy
σ
σe = P ξ − ξ η + η
2
(3)
σe2
=
ρ
⎣
91 93 94 95 96 97 98 100 102 103 105 107 108 109 110
(4)
112
αi j can
114
113 115 116 117 118
c y y c zz − c 2yz
c xz c yz − c xy c zz c xx c zz − c 2xz
symm
119
⎤
0 0
⎦
c xx c y y − c 2xy
120
(5)
121 122 123 124
where:
c xx = 1 − 3 1 −
c yy = 1 − 3 1 − c xy = −
89
111
12
αxx αxy αxz ⎣ αxy α y y α yz ⎦ αxz α yz αzz 1
88
106
⎤
⎡
87
104
S i j S kl ε˙ kl
For both the FT and DT models the instantaneous moduli be expressed in the form:
⎡
86
101
In Eqs. (2) and (3) (G , λ, G s , λs ) are the elastic and the secant Lamè coefficients, G t is the tangent shear modulus, S i j are the stress deviator components and σe is the equivalent stress that, for the load conditions considered, assumes the value: 2
85
99
(2)
2 e
and the deformation theory of plasticity with the Hencky constitutive relation:
84
92
S i j S kl ε˙ kl
σ˙ i j = 2G s ε˙ i j + λs δi j ε˙ kk − 3(G s − G t )
83
90
(1)
where E is the elastic Young modulus and αi j are instantaneous moduli depending on the plasticity theory considered to model material behavior. Here two plasticity theories are considered, namely the incremental, or flow, theory of plasticity, with the Prandtl–Reuss constitutive equation:
σ˙ i j = 2G ε˙ i j + λδi j ε˙ kk − 3(G − G t )
82
1 2
Et E¯ Et E¯
s2y
4 s2x 4
1 − (1 − 2ν )
125
,
126
c xz = 0,
127 128
, Et E
c yz = 0,
−3 1−
Et E¯
129
sx s y 2
130
131
,
132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.3 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
3
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
Fig. 2. Ramberg–Osgood stress–strain relation.
34 35 36
c zz =
37 38 39 40 41 42 43 44 45 46 47 48
ρ
101
Et
− 1 σx sx = = −ξ ξ 2 − ξ η + η2 2 σe − 1 σy sy = = −η ξ 2 − ξ η + η2 2 σe
55
E¯ = E ;
G¯ = G =
E
58
in flow theory, whereas assume the secant values
E¯ = E s ;
G¯ = G s =
61 62 63 64 65 66
E 2(1 + ν ) + 3( EE − 1)
(9)
ε=
σe E
+k
σ0 σe E σ0
Et
= 1 + nk
σe σ0
σe =1+k Es σ0 E
n−1
(10)
has been adopted to model the stress–strain elastoplastic curve. In Eq. (10) σ0 represents the nominal yield stress, σe is the equivalent stress as defined in Eq. (4) and (k, n) are constant describing the shape of the stress–strain elastoplastic curve, as shown in Fig. 2. In such a hypothesis the ratio of the elastic modulus E and the
106
n−1
107
(11)
2 1
ε˙ y = s˙ y, y +
2
111 112 113 114 115 116 117 118 119 120 122 123
(12)
s˙ 2y ,x + s˙ 2z,x
110
121
The non-linear expressions for the effective strain components contributing to the strain energy can be put in the following form:
ε˙ x = s˙ x,x +
108 109
whose dependent on the (k, n) coefficients are reported, in terms of σe /σ0 ratio, in Fig. 3, where it is possible to observe how, for low values of σe /σ0 , both E s and E t assume the same value of the elastic constant E, then both decrease when σe /σ0 achieve a limit value, with a velocity depending on the assigned values of (k, n). The critical behavior of the plate has been analyzed within the framework of the Kirchhoff–Love kinematic model, hence the dis˙ (x, y , z, t ) = [˙sx s˙ y s˙ z ] T can be expressed placement field velocity U in terms of the three components U˙ = U˙ (x, y , t ), V˙ = V˙ (x, y , t ) ˙ (x, y , t ) of the middle surface as: ˙ =W and W
1
103 105
˙ s˙ z = W
n
102 104
˙ ,y s˙ y = V˙ − z · W
in deformation theory. In this paper the Ramberg–Osgood elastoplastic equation
E
˙ ,x s˙ x = U˙ − z · W
s
59 60
(8)
2(1 + ν )
56 57
(7)
and the constitutive terms ( E¯ , G¯ ), that are the elastic values
53 54
(6)
are constant depending on both the stress state and the incremental theory considered, through the values assumed by
51 52
corresponding secant and tangent values ( E s , E t ) can be obtained in the forms:
G¯ E = c zz c xx c y y − c 2xy Et
49 50
100
s˙ 2y , y + s˙ 2z, y
124 125 126 127 128 129 130
131 132
JID:AESCTE
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.4 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
4
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
Fig. 3. E t / E, and E s / E curves varying the Ramberg–Osgood coefficients n, k.
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
99
γ˙xy = s˙ x, y + s˙ y,x + s˙ y,x s˙ y, y + s˙ z,x s˙ z, y γ˙xz = s˙ x,z + s˙ z,x + s˙ y,x s˙ y,z + s˙ z,x s˙ z,z γ˙ yz = s˙ y,z + s˙ z, y + s˙ y, y s˙ y,z + s˙ z, y s˙ z,z
56
˙ )= δΩ(U˙ , V˙ , W
59
ΩU =
64 65 66
(14)
1 2
+
h/2
ΩV =
2
σx h
106
(17)
˙ ,xx − h2 /12 W ˙ ,xxy y + ξ Ph W ˙ , y y − h2 /12 W ˙ ,yyyy = 0 +ηPh W
107 108 109 110
111 112 113 114
(18)
115
˙ can be obtained. for the out-of-plane displacement W It is worth noting that, because in Eq. (13) the non-linear contribution of the displacement v˙ has been considered, both the Eqs. (17) and (18) are slightly different to the classical equation reported in literature, with new terms usually neglected under the von Kármán hypothesis.
116
(15)
(16)
Eqs. (17) and (18) can be analytically solved by means of the Kantorovich technique with which, through a separation of (x, y ) variables, a system of partial differential equations can be transformed in a set of ordinary equations for which a closed form solution exists. If simply supported, the boundary conditions on the edges x = ±b/2 are automatically satisfied by the following displacement field:
117 118 119 120 121 122
s2y ,x
+
s2z,x
+ σyh
s2y , y
+
s2z, y
124 125
−h/2
105
for the in-plane eigenmode components (U˙ , V˙ ), and:
αxx W˙ ,xxxx + 2(αxy + 2αzz ) W˙ ,xxy y + α y y W˙ , y y y y
102 104
− h(ξ P V˙ ,xx + η P V˙ , y y ) = 0
Eh3 /12
101
123
σ˙ x ε˙ x + σ˙ x ε˙ x + τ˙xy γ˙xy + τ˙xz γ˙xz + τ˙ yz γ˙ yz dz
1
Eh αxx U˙ ,xx + (αxy + αzz ) V˙ ,xy + αzz U˙ , y y = 0 Eh α y y V˙ , y y + (αxy + αzz )U˙ ,xy + αzz V˙ ,xx
100
103
3. The Kantorovich method
represents the strain energy functional for the plate:
62 63
δ(ΩU + Ω V )dΓ = 0
where:
60 61
Γ
57 58
(13)
In Eq. (13) and in the following, the comma before subscripts indicates partial derivative, e.g. s˙ i , j is equivalent to ∂ s˙ i /∂ x j . Unlike the classical von Kármán strain–displacement hypothesis, which considers only the non-linear contribution deriving from the out-ofplane displacement s˙ z , in Eq. (13) the additional non-linear terms related to the in-plane component s˙ y are also considered. Even if negligible for the buckling analysis of flat plate, such a contribution can considerably affect the critical response of stiffened plates and prismatic beams, when buckling mode involves flexo-torsional or global flexural buckling displacement (Figs. 4A, 4B), remaining negligible in case of local flexural mode (Fig. 4C). The equilibrium state is thus sought by invoking the principle of minimum potential energy:
54 55
Kirchhoff’s kinematical model represented in Eq. (12), the minimization of the deformation energy with respect to arbitrary vari˙ ) returns the following local equations: ation of (U˙ , V˙ , W
is the potential energy for the plate subjected to uniform in-plane compressive stress and δ is the variational symbol. Substituting the
126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.5 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
Fig. 4. Pictorial scheme of global flexural (A), flexo-torsional (B) and local flexural (C) buckling mode.
14
80
15 16 17 18 19 20 21 22
81
U˙ (x, y ) = u˙ ( y ) cos
mπ x
82
b
83
mπ x V˙ (x, y ) = v˙ ( y ) sin
˙ (x, y ) = w ˙ ( y ) sin W
84
b mπ x
85 86
(19)
b
87
˙ ( y ) are unknown functions, and ˙ =w where u˙ = u˙ ( y ), v˙ = v˙ ( y ), w m is the assigned number of half-waves in x-direction related to the buckling mode. Indicating with λ = b/mπ the half-wave length in buckling mode, in virtue of the displacement field (19), the system (17) becomes the following linear, ordinary differential equation system related to the in-plane displacements (u˙ , v˙ ):
88
95
30
λ2 αzz u˙ , y y + λ(αxy + αzz ) v˙ , y − αxx u˙ = 0
31
λ2 (α y y − ηPh/ A ) v˙ , y y − λ(αxy + αzz )u˙ , y
23 24 25 26 27 28 29
32 33 34
− (αzz − α Ph/ A ) v˙ = 0
36 37 38 39
˙ ,yyyy + (α y y − ηκ P ) w
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
92 93 94
98
(20)
η h − ξ κ λ2 P
˙ =0 ˙ , y y + λ4 αxx D − ξ hλ2 P w − 2λ2 D (αxy + 2αzz ) w
(21)
99
Fig. 5. (a) Assembled thin-walled plates, geometry and local reference system. (b) Example 1, Buckling load factor κ versus α /β aspect ratio for simply supported plates with three different thicknesses parameter ζ .
which admits a closed form solution in the form:
u˙ = λ11 c 1 e ρ1 y + λ12 c 2 e ρ2 y + λ13 c 3 e ρ3 y + λ14 c 4 e ρ4 y
T˙ (x, y ) = t˙( y ) cos
v˙ = λ21 c 1 e ρ1 y + λ22 c 2 e ρ2 y + λ23 c 3 e ρ3 y + λ24 c 4 e ρ4 y
˙ = δ11 d1 e ω1 y + δ12 d2 e ω2 y + δ13 d3 e ω3 y + δ14 d4 e ω4 y w
(22)
depending on the known constant values (λi j , δ1 j ) and (ρ j , ω j ), which can be derived from the eigenvectors and eigenvalues matrices associated with the coefficients of (20) and (21), and by the unknown vector x = [c i di ] T , whose components can be determined by imposing the proper in-plane and out-of-plane boundary conditions on the edges y = ±a/2. Such boundary conditions, if collected in a single, square matrix in the form:
L·x=0
(23)
allow the analysis of single, flat, rectangular plates through the eigenvalue analysis of the matrix L with which it is possible to derive, for the simplest case, also a closed form solution [8]. The stability analysis of thin-walled structures with an arbitrary cross-section, as reported in Fig. 5, require the additional definition ˙ , y and of a local stiffness matrix k such that, indicating with ϑ˙ = w ˙ = [u˙ v˙ w ˙ ϑ˙ ] y=±a/2 the nodal displacecollecting in the vectors u ˙ ] y=±a/2 the corresponding dual forces, ment, and in p = [t˙ n˙ p˙ m it is possible to describe the nodal force–displacement relationship in the form
p=k·u
(24)
100 101 102 103
Under the displacement hypothesis (19), the corresponding Kirchhoff forces can be expressed as:
65 66
91
97
40 41
90
96
˙ whereas Eq. (18) becomes an ordinary differential equation in w:
35
89
˙ (x, y ) = n˙ ( y ) sin N P˙ (x, y ) = p˙ ( y ) sin
107
b mπ x
108 109
b mπ x
110 111 112
b mπ x
113
(25)
b
n˙ = A
αzz u˙ , y + αzz λ−1 v˙
p˙ = D
117 118
119
α y y w˙ , y y y − λ−2 αxy + 4λ−2 αzz − ηPh w˙ , y
116
α y y v˙ , y − αxy λ−1 u˙
114 115
˙ ( y ) are: ˙ =m where t˙ = t˙( y ), n˙ = n˙ ( y ), p˙ = p˙ ( y ) and m t˙ = A
105 106
mπ x
˙ (x, y ) = m ˙ ( y ) sin M
104
−2
120
121 122
(26)
123
In Eq. (26) A = E · h/(1 − ν 2 ) and D = E · h3 /12(1 − ν 2 ) representing respectively the extensional and the flexural plate stiffness. ˙ ) reported Substituting in Eq. (26) the displacement field (u˙ , v˙ , w in Eq. (22) it is possible to express the relationship between the nodal force and corresponding displacement, in the form reported in Eq. (24) [14]. Finally, it is possible to obtain the global stiffness matrix K in a FEM-like procedure, and the related plastic buckling parameter
125
˙ =D m
α y y w˙ , y y − αxy λ
˙ w
124 126 127 128 129 130 131 132
JID:AESCTE
AID:3273 /FLA
6
[m5G; v1.149; Prn:9/03/2015; 15:29] P.6 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
Fig. 6. Example 1, Buckling load factor
33
κ versus ξ/η load ratio for simply supported plates with three different thicknesses parameter ζ .
99
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
100
P as the lowest eigenvalue of K, and the corresponding buckling mode as the related eigenvector. However, for the non-linearity of the matrix K( P ), both because P is contained through (λi j , δ1 j ) in the general solution (22) and because it is contained, through σe , in the elastoplastic constitutive coefficients ( E t , E s ), the eigenvalue solution requires a double, nested iterative strategy based on a trial and error procedure. Assuming the Shanley concept of continued loading, i.e. no unloading occurs at the instant of the plastic buckling, the outer iteration starts by assuming a value of P = P 0 small enough to not cause plastic deformation and evaluating the related instantaneous moduli αi j ( P 0 ). With such values the inner iteration solves the corresponding eigenproblem, in which the critical load P 1 represents the lower eigenvalue. If ( P 1 − P 0 )/ P 1 > err, where err is a prescribed error bound (in the proposed examples an err = 0.001 has been considered), the outer iteration restarts considering a new P 0 = (1 − η) P 0 + η P 1 (0 < η < 1). The value P i such that ( P i − P i −1 )/ P i < err is the required elastoplastic buckling load.
53 54 55
4. Examples
56 57 58 59 60 61 62 63 64 65 66
In the numerical analyses an Aluminum alloy AL 14S-T6, characterized by a Young modulus E = 72 400 MPa, a Poisson ratio ν = 0.32, a yield stress σ0 = 430 MPa and a value for the Ramberg–Osgood coefficients n = 10.9, k = 0.3485, has been considered. The results are reported in terms of the buckling coefficient κ = P cr (1 − ν 2 )/ζ E, in which ζ = π 2 h2 /12a2 is a thickness parameter. Here the three different values ζ = 0.0001, ζ = 0.001 and ζ = 0.002 (corresponding to a value of the thickness/length ratio h/a = 0.011, h/a = 0.0349 and h/a = 0.0493, respectively) have been considered.
4.1. Rectangular plates
101 102
In order to validate the numerical procedure related to the proposed model, a first analysis concerning the elastoplastic buckling of a flat, rectangular, simply supported plate of length a, width b, thickness h, subjected to uniaxial compression has been investigated, and compared with the analytical solution:
N cr =
π 2 D b 2 m2 b2
a2
αxx + 2(αxy + 2αzz ) +
a2 b 2 m2
αyy
N cr =
2
π D b m b2
a2
2
+2+
a
2
b 2 m2
104 105 106 107 108
(27)
derived in [5] as an elastoplastic generalization of the corresponding elastic solution: 2
103
109 110 111 112 113 114
(28)
reported, among others, in [15]. In both Eqs. (27)–(28) m represents the longitudinal half-wave number characterizing the buckling mode, whereas αi j reported in Eq. (27) are the elastoplastic coefficients defined in Eq. (1). Fig. 6 reports the buckling coefficient κ obtained varying the b/a aspect ratio related to the first 5 modes (m = 1, 2, . . . , 5) by means of both Eq. (27) and the proposed formulation, and for both the FT and the DT models (together with the ones derived from a linearly elastic model, here reported for comparison purpose). The results show an excellent agreement between analytical and numerical analysis. It is worth noticing that, because both the analyses derive from the same kinematical hypothesis on the sinusoidal behavior of the out-of-plane displacement w, the comparison only measures the accuracy of the stiffness matrix evaluation and of the double iterations required for solving the eigenvalue problem numerically. The figure also shows the coincidence between
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.7 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
the minimum values are κmin = 3.32 when a/b = (0.72, 1.45, 2.16, 2.88, 3.6) in FT and κmin = 2.67 when a/b = (0.82, 1.62, 2.49, 3.3, 4.1) in DT. Fig. 7 reports the critical response of the same plate subjected to biaxial loads obtained varying the load ratio ξ/η , compared with analytical ones proposed by Durban and Zuckerman [8]. Again the results are in good agreement, including the growing discrepancy between the results computed by the two elasto-plastic models, which give rise the plastic buckling paradox for which “bifurcation loads calculated using the simplest flow theory of plasticity consistently overestimated buckling loads of plates and shells obtained in tests. Calculations based on the less respectable deformation theory of plasticity gave reasonably good agreement with test results” Bijlaard, [6]. For a better focus of such paradox, Figs. 8 and 9 report the limit values σcrmin /σ0 and σcrmax /σ0 for all of the examples considered and their position on the Ramberg–Osgood curve – Fig. 8 – and on the corresponding ( E s / E , E t / E ) curves – Fig. 9 – adopted in the analysis, showing that if the minimum values coincide for any analysis performed, the differences on the maximum values are equal to zero when the critical load is contained in the elastic range of the stress–strain curve, which otherwise increase, reaching in the worst case, obtained with the thickness parameter ζ = 0.002 and the load parameters ξ = −η = 1, a ratio close to three (i.e. the FT analysis returns a critical load triple of the corresponding value obtained through the DT analysis). For a better comprehension of the cause involving such paradox is essential to undertake a closer analysis of the instantaneous moduli αxx , 2(αxy + 2αzz ) and α y y which appear in Eq. (18), and point out the differences on such coefficients derived from the use of the two theories. Substituting Eq. (6) in (5), and taking into account the values of (s x , s y ) reported in Eq. (7) the following expressions for such coefficients can be obtained:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
1 Et ξ 2 ξ 2 − ξ η + η2 αxx = 4−3 1− ρ E¯ 1 Et αyy = 4−3 1− η2 ξ 2 − ξ η + η2 ρ E¯ 1 Et Et αxy = 2 − 2(1 − 2ν ) −3 1− ηξ ξ 2 − ξ η + η2 ρ E E¯
34 35 36 37 38 39 40 41 42 43
αzz =
44
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 110
E
elastic and elastoplastic results for the thinner plate ζ = 0.0001, for which in any performed analysis the coefficient κ assumes a minimum value κmin = 4 when a/b = (1, 2, 3, 4, 5), and a maximum value κmax = 4.5 when a/b = 1.41. Again, the curve related to m number of half-waves intersects the curve related to m + 1 half-wave when a/b = (1.41, 2.45, 3.46, 4.47). When the plate becomes thicker (ζ = 0.001), elastic and elastoplastic curves diverge, even though they maintain the same general aspects. Both the FT and DT curves are characterized by a minimum value, that is κmin = 3.92 when a/b = (0.97, 1.9, 2.91, 3.89, 4.84) for the FT case and κmin = 3.87 when a/b = (0.98, 1.95, 2.92, 3.9, 4.87) in DT analysis. The intersections between the curves also change: they are a/b = (1.32, 2.35, 3.34, 4.32) and a/b = (1.34, 2.37, 3.36, 4.34) in FT and DT case, respectively. The differences between the two models increase for the thickest plate (ζ = 0.002). In such a case
whose values depend, through ( E t , E s ), on the actual load coefficient P and, through the (ξ, η) coefficients, on the load ratio existing between the assigned nx and n y in-plane loads. The latter affect the slope of αi j as shown in Figs. 10–12. Coherently with the elastoplastic theories adopted, it is reasonable to expect that coefficients decline as the load parameter P increases. Figs. 10–12 show how it occurs for any coefficient, excluding both αzz ( F T ), that assume the constant value G¯ / E = 2(1 + ν ), and αxy which, since it contains a term which depend on the product ηξ , can increase if such coefficients are characterized by different signs. As a consequence, if subject to a compressive force in one direction and a tensile force in the other, the instantaneous modulus 2(αxy + 2αzz ) related to the FT model presents an unnatural evolution, contrary to what happens in the DT model, in which the increasing values of αxy are canceled by the decreasing values of αzz , the latter depending on the secant value G¯ = G s reported
48
67
109
Fig. 7. Example 1, Position of maximum/minimum σcr /σ0 on the Ramberg–Osgood curve for any load ratio and thickness parameter considered.
47
51
G¯
111
E
46
50
Et ρ = + (1 − 2ν ) 3 − (1 − 2ν ) ¯G E 2 Et 2 −3 1− ηξ ξ − ξ η + η E¯
45
49
7
112 113 114
(29)
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
8
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.8 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
Fig. 9. Example 1, Generalized constitutive variables
47
αi j for ξ = η = 1.
48 50 52 53
56 57 58 59 60 61 62 63 64 65 66
115 116
51
55
113 114
4.2. C-section columns
49
54
112
Fig. 8. Example 1, Position of maximum/minimum σcr /σ0 on the ( E t , E s )/ E curves for any load ratio and thickness parameter considered.
in Eq. (9). Such anomaly can be considered a concurrent cause of plastic buckling paradox. Fig. 13 shows how the differences between the two models decrease drastically by simply substituting G¯ = G s in the FT model. Finally, the results related to a single plate also show how the weight of the quadratic v-terms represented in Eq. (13) is negligible, consistently with von Kármán hypothesis for which, due to the buckling mode involving exclusively out-of-plane displacement, only the corresponding nonlinear strain related to the w-terms affect the buckling analysis of flat plates.
The next example analyzes the critical behavior of the thinwalled C-section columns represented in Fig. 14. As it is derived by a closed form solution, the proposed procedure does not depend on the discretization adopted in the analysis, and a model capable to determining the critical load related to any of the flexural, torsional or flexo-torsional mode, both local and global, can be obtained with the minimum number of nodes required to represent the section geometry that, for the proposed example, are the six nodes (A . . . F), again represented in Fig. 14. The related global stiffness matrix K considered in the eigenvector problem has then dimension (24 × 24). The reduced dimension of the stiffness matrix leads to a very fast numerical code, particularly suitable for parametric analysis. For instance, Fig. 15 shows the elastic buckling parameter κ obtained varying the length ratio b/a between 1 to 100, again considering the three thickness parameters ζ = 0.0001,
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.9 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32 33
Fig. 10. Example 1, Generalized constitutive variables
αi j for ξ = 0, η = 1.
98 99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 11. Example 1, Generalized constitutive variables
αi j when ξ = −1, η = 1.
132
JID:AESCTE
10
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.10 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35 36
101
Fig. 12. Example 1, Buckling factor
κ versus ξ/η load ratio considering G = G s in FT analyses.
37
103
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
102 104
ζ = 0.001 and ζ = 0.002. The numerical analysis has been performed considering both the classical von Kármán strain model, in the figures indicated with the apex c (classical) and the enriched model which takes into account the nonlinear contribute of the v-component, in the figure reported with the apex e (enriched). The results show the contrasting behavior for the three cases. When ζ = 0.0001, according with the length value, for short structures the buckling mode is web-triggered local (Fig. 16A) and involves only out-of-plane displacements (local-flexural mode), and a number of half-waves increasing – from 1 to 30 – with the length ratio. In such a phase the critical load remains substantially constant, equal to the minimum value related to the buckling curve related to the number of half-waves activated in buckling mode. For long structures the buckling mode becomes global, purely torsional (Fig. 16D) involving comparable out-ofplane and in-plane displacements. As a consequence in such phase the von Kármán simplification fails, and numerical analyses performed considering the non-linear strain term reported in Eq. (13) return a critical load sensibly inferior, and a lower value of the ratio b/alim characterizing the transition between local and global mode, that is equal to b/alim = 28 (with a critical mode characterized to 30 half-waves) under the von Kármán hypothesis, whereas it is equal to b/alim = 19 (and 22 half-waves) otherwise. For ζ = 0.001 and ζ = 0.002 the mode is flexural-torsional (Fig. 16B–C) also for short structures, and the differences between the elastic analyses performed with the two deformation models
are visible for any chosen ratio b/a. For long structures the mode becomes global, purely torsional in a first phase (Fig. 17D), purely flexural in a final phase (Fig. 17E), involving in any case comparable in-plane and out-of-plane displacement. When ζ = 0.0001 the elastoplastic results again coincide with those deriving from an elastic analysis. Fig. 17 shows the elastoplastic analyses performed for ζ = 0.001 and ζ = 0.002. Again the DT model, coupled with the enriched deformation model, returns the lowest critical load. However, the analyses performed considering the FT model follow the same critical path, and do not present the different paths observed on the single plate reported in Fig. 13. It is interesting to observe that for short structures the differences between the two elastoplastic models are more pronounced under the von Kármán hypothesis. Such differences tend to disappear for long structures, when the elastoplastic buckling load is closer to the elastic buckling load. Instead, in such a phase the differences generated by considering the contribute of v in Eq. (13) emerge.
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
4.3. Close-sections columns
124 125
Finally, a buckling analysis of five different thin-walled closed section beams, whose form and dimensions are reported in Fig. 18, has been performed. The sections are labeled by the number n of the sides characterizing the section, and the dimensions are such that all the figures are inscribed in a circle of radius a, so that the beams fill the same space, whereas the thickness variation is such that the sections maintain the same area. The analyses are then
126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.11 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
11
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44 45 46
110
Fig. 13. Example 1, Buckling load factor κ versus α /β aspect ratio for simply supported plates with three different thickness parameter ζ .
111 112
47
113
48
114
49
115
50
116
51
Fig. 15. Example 2, linear-elastic buckling load factor κ versus α /β aspect ratio considering both the Classical von Kármán (c) and the enriched (e) strain–displacement model.
52 53 54
118 119 120
55 56 57 58 59 60 61 62 63 64 65 66
117
Fig. 14. Example 2, geometrical aspect ratio.
performed, again in terms of κ vs. b/a, on equivalent beams, in terms of overall dimensions and weight. Fig. 19 reports the buckling load obtained through elastic analyses. The results, reported in Fig. 19 considering linearly elastic behavior, show a general behavior very close to what was observed in the previous example, i.e. short plates characterized by a critical load almost constant and a local, purely flexural, critical mode (Fig. 20A), with a number of half-waves growing with the length whereas, when the length ratio b/a reaches a limit value b/alim , the critical mode becomes global, and the critical load-factor κ
121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:3273 /FLA
12
[m5G; v1.149; Prn:9/03/2015; 15:29] P.12 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
Fig. 16. Example 2, local flexural (A), local flexo-torsional (B, C), global torsional (D) and global flexural (E) buckling mode.
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
82 83
decreases for long plates. Because close sections show a greater torsional stiffness, purely flexural (Fig. 20B) behavior also characterizes the critical mode of the thickest plates, that do not present the flexo-torsional buckling load observed in the previous example. The limit value b/alim depends on the rigidity of the plate, and decreases with increasing the number n of the sides of the section. Fig. 19 also shows how for short plates, when the mode is local and the corresponding critical load is affected by the ratio between the length and the thickness of the single plate forming the section, the critical load-factor κ is greater for the triangle-shape section and it decreases as n increases. Vice versa, for long plates, when the mode becomes global and the corresponding critical load depends on the global stiffness of the section, the triangle-shaped section presents the lower critical factor, whereas the greatest is obtained for the square-shaped section. The elastoplastic behavior is also affected by the superior torsional stiffness, as reported in Fig. 21 for the thickest plates (ζ = 0.002). It shows a general behavior, in terms of critical path and associated critical mode, not so different to the one observed in elastic analysis. As in the previous examples, in all the performed analyses the DT model furnishes the lower critical coefficient, albeit the differences between the two models do not diverge, remaining almost constant varying the b/a ratio. It is worth noting that, increasing the number n, the gap between elastic and elastoplastic results increases, becoming almost double for n = 8.
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
44 45
110
5. Conclusion
111
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
81
112
Based on the Kirchhoff plate theory and making use of the Kantorovich technique, the present work provides a nonlinear model capable of deriving the critical behavior of prismatic beams and stiffened plates made of nonlinear, elastoplastic materials. A numerical procedure based on a closed form solution of the equilibrium equations has been derived and, since in such methods the results are not affected by numerical instabilities and/or meshsize, the approach can be adopted with any required accuracy. For validation purposes, the proposed formulation was first used to determine the critical behavior of rectangular plates. The comparison with exact solutions available in literature showed the correctness of the proposed procedure and, through a more accurate analysis of the elastoplastic coefficients derived by both DT and FT theories, it has shed a new light on the causes of the plastic buckling paradox. The analyses performed on a C-section and five different closed sections columns have shown the inadequacy of von Kármán hypothesis on the determination of both elastic and elastoplastic buckling loads when the buckling mode involves comparable in-plane and out-of-plane displacements and highlighted the need, in such cases, for a refined model. Finally, the sensitivity
113 114 115 116 117 118
Fig. 17. Example 2, elasto-plastic buckling load factor κ versus α /β aspect ratio considering both the Classical Von Kármán (c) and the enriched (e) strain-displacement model.
119 120 121 122
analyses have highlighted the possibility of following the transition from global to local buckling modes, and consequently the possibility of determining the corresponding variations in critical loads, by modifying key geometrical parameters of the structures and therefore the method can be helpfully exploited in design optimization procedures.
123 124 125 126 127 128 129
Conflict of interest statement
130 131
None declared.
132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.13 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
13
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10 11
76
Fig. 18. Example 3, geometrical aspect ratio.
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 19. Example 3, linear-elastic buckling load factor
κ versus α /β aspect ratio.
132
JID:AESCTE
14
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.14 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15 16
81
Fig. 20. Example 3, local and global buckling mode.
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 21. Example 3, linear-elastic and elasto-plastic buckling load factor
κ versus α /β aspect ratio considering both FD and DT models.
132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.15 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
References
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
[1] Z.P. Bazant, L. Cedolin, Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories, Courier Dover Editions, 1991. [2] C.M. Wang, C.Y. Wang, J.N. Reddy, Exact Solution for Buckling of Structural of Members, CRC Series in Computational Mechanics and Applied Analysis, 2009. [3] P. Weißgraeber, C. Mittelstedt, W. Becker, Buckling of composite panels: a criterion for optimum stiffener design, Aerosp. Sci. Technol. 16 (2012) 10–18. [4] Matthias Heitmann, Peter Horst, A new analysis model for the effective stiffness of stiffened metallic panels under combined compression and shear stress, Aerosp. Sci. Technol. 10 (2006) 316–326. [5] G.H. Handelmann, W. Prager, Plastic buckling of a rectangular plate under edge thrusts, naca-report-946, 1949, pp. 479–507, http://naca.central.cranfield.ac.uk/ reports/1949/naca-report-946.pdf. [6] P.P. Bijlaard, On the plastic buckling of plates, J. Aeronaut. Sci. 17 (1950) 485–493. [7] F. Bleich, Buckling Strength of Metal Structures, McGraw-Hill, New York, 1952. [8] D. Durban, Z. Zuckerman, Elastoplastic buckling of rectangular plates in biaxial compression/tension, Int. J. Mech. Sci. 41 (1999) 751–765.
15
[9] P. Le Grognec, A. Le Van, Some new analytical results for plastic buckling and initial post-buckling of plates and cylinders under uniform compression, ThinWalled Struct. 47 (2009) 879–889. [10] J.W. Hutchinson, Plastic buckling, Adv. Appl. Mech. 14 (1970) 67–144. [11] X.M. Wang, J.C. Huang, Elastic/plastic buckling analyses of rectangular plates under biaxial loadings by the differential quadrature method, Thin-Walled Struct. 47 (2009) 879–889. [12] W. Zhang, X.W. Wang, Elastic/plastic buckling analysis of thick rectangular plates by using the differential quadrature method, Comput. Math. Appl. 61 (2011) 44–61. [13] M. Kadkhodayan, M. Maarefdoust, Elastic/plastic buckling of isotropic thin plates subjected to uniform and linearly varying in-plane loading using incremental and deformation theories, Aerosp. Sci. Technol. 32 (2014) 66–83. [14] E. Ruocco, M. Fraldi, Critical behavior of flat and stiffened shell structures through different kinematical models: a comparative investigation, Thin-Walled Struct. 60 (2012) 205–215. [15] S.P. Timoshenko, J.M. Gere, Theory of Elastic Stability, McGraw–Hill, New York, 1961.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132
[m5G; v1.149; Prn:9/03/2015; 15:29] P.1 (1-15)
Aerospace Science and Technology ••• (••••) •••–•••
1
Contents lists available at ScienceDirect
2
67 68
3
Aerospace Science and Technology
4 5
69 70 71
6
72
www.elsevier.com/locate/aescte
7
73
8
74
9
75
10
76
11 12 13 14
Elastoplastic buckling analysis of thin-walled structures
80 81
16
82
a r t i c l e
i n f o
83
a b s t r a c t
18 19 20 21
84
Article history: Received 21 July 2014 Accepted 4 March 2015 Available online xxxx
22 23 24 25 26 27
78 79
Eugenio Ruocco 1
15 17
77
Keywords: Elastic–plastic buckling Kantorovich technique Buckling of stiffened plates Ramberg–Osgood law Buckling paradox
28 29
The study is concerned with the elastoplastic buckling of thin-walled beams and stiffened plates, subjected to in-plane, uniformly distributed, uniaxial and biaxial load. The ruling differential equations have been solved analytically by using the Kantorovich technique and the obtained displacement field has been employed in a general procedure that, by using the framework derived by the finite element method, is able to analyze the elastoplastic buckling behavior of prismatic beams and stiffened plates with arbitrary cross-section. The inelastic effect is modeled through a stress–strain law of the Ramberg– Osgood type, and both the incremental deformation theory and the J 2 flow theory are here considered. The reliability of the numerical procedure is illustrated for rectangular plates, and the contradicting results obtained by using the two plastic theories are discussed in detail. Finally, the performance of the method is illustrated through the analysis of a C-section and five different closed section columns. © 2015 Published by Elsevier Masson SAS.
85 86 87 88 89 90 91 92 93 94 95
30
96
31
97
32
98
33 34
99
1. Introduction
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
Due to their relevance in the design of thin structures, the critical behavior of plates, both flat and stiffened, has been extensively studied in the past years, as documented by the number of text [1,2] and papers [3,4] on the issue. More recently, a number of authors have also investigated the effect of elastoplastic constitutive models on the prediction of the buckling load of thin structures, generally by employing models based on both the incremental J 2 deformation theory (DT) and the J 2 flow theory (FT). Starting from the pioneering works of – among others – Handelmann and Prager [5], Bijlaard [6] and Bleich [7] in the 1950s, different aspects of buckling behavior in elastoplastic range have been analyzed, and collected in papers dedicated to the analytical investigation of buckling [8] and postbuckling [9] of simple structures. Some of these works examine the causes of the “plastic buckling paradox”, as it was first revealed in a series of theoretical investigations on the buckling of a simply supported flat plate of infinite length, through both the FT [5] and the DT [6]. A comparison of the results shows that “the buckling load obtained from the “more rational” FT model always far exceeds the experimental observation, whereas the results obtained adopting the “simplest”, total strain DT, predict buckling loads lower than FT model and in any case in closest agreement with those obtained experimentally. Moreover, the differences between the two models increase as the level of plasticity
60 61 62 63 64 65 66
1
E-mail address: [email protected]. Tel.: +39 081 5010266.
http://dx.doi.org/10.1016/j.ast.2015.03.001 1270-9638/© 2015 Published by Elsevier Masson SAS.
increases, giving in some cases results with different order of magnitude” [10]. As reported in [8], such aspects drastically affect the behavior and the load carrying capacity of flat plates subject to uniaxial or biaxial in-plane load, regardless of the boundary conditions. For instance, Durban and Zuckerman found, for specified compression/tension ratios, an optimal loading path for the DT that has no correspondence in FT. The works of Wang et al., which analyze the elasto-plastic buckling of thin [11] and thick [12] flat plates by differential quadrature method, confirm the results obtained in [8] and show that the paradox also involves flat plates of different geometry, such as triangular or elliptical, with different boundary conditions. A comparison of incremental DT and flow rule has been recently proposed in [13], where a new algorithm, based on a Generalized Differential Quadrature (GDQ) discretization technique to solve the out-of-plane plate equation, has been presented. The aim of this work is to present a new analytical formulation to define the critical behavior of both flat and stiffened plates, as well as of thin walled columns with open or closed cross sections, and characterized by non-linear elastoplastic behavior. To this end a systematic procedure based on a finite element algorithm capable of accounting for both local and global phenomena has been developed. As they are analytical, the obtained results do not depend on the discretization adopted, and reliable solutions can be achieved with the minimum number of elements required to represent the structural geometry, with the double effect of stabilizing the numerical procedure and accelerating the computing time required for the eigenvalue and eigenvector analysis.
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
2
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.2 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
Fig. 1. Geometry and load condition of a rectangular plate.
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
78
For validation purposes, the proposed procedure has been applied to flat plates, and the results compared to the analytical and/or numerical solutions available in literature. In such a preliminary analysis particular attention is paid to the evaluation of the instantaneous elastic moduli, and to the differences between them deriving from the two models adopted. The analysis indicates that the (correct) use of the constant value of Lamè modulus G in FT shows, as consequence, an irregular behavior of the corresponding generalized elastoplastic modulus which, when the ratio between the in-plane loads is greater than a critical value, in plastic range increases its stiffness, in contrast to what happens to the other modulus. Having adopted the secant modulus G s , the corresponding modulus calculated in the DT does not present the same anomaly. Observing that just in correspondence of such critical ratio the buckling path corresponding to the two theories begins to diverge, it is possible to assert that the irregular behavior of the elastic modulus related to the G-coefficient is a cause, or at least a concomitant cause, of the plastic buckling paradox. Another aspect considered in this work is the adequacy of the von Kármán hypothesis in the negligible contribution of secondorder strain terms related to in-plane displacement. In a previous work devoted to the influence of such terms on the elastic buckling of shell and both flat and stiffened plates [14], the authors demonstrate how, when the critical mode involves comparable in-plane and out-of-plane displacement, the omission of the non-linear strain terms related to in-plane displacements can considerably overestimate the critical load. If such cases never happen for flat plates, for which the von Kármán simplification is thus correct, it is quite common for stiffened plates and prismatic beams, when buckling mode involves torsional, flexo-torsional or global flexural buckling displacement. In the present work is – to the author’s knowledge for the first time – the influence of such nonlinear terms on the elasto-plastic buckling of stiffened plates is presented and discussed. In order to capture the inelastic effect, the constitutive behavior is here modeled through a stress–strain law of the Ramberg– Osgood type that, with the basic equations ruling the elastoplastic buckling of Kirchhoff plates subjected to uniaxial or biaxial compression and under both the FT and DT hypothesis, are introduced in the following section. Section 3 reports an analytical solution based on the Kantorovich technique, and its generalization in a FEM-like procedure for the buckling analysis of structures ascribable to plates rigidly connected together along their edges. The same section includes the numerical procedure necessary to solve the eigenproblem associated to the corresponding elastoplastic stiffness matrix. In Section 4 parametric analyses related to a single plate with different boundary and load conditions, a simply supported open C-section and five different thin-walled closed section beams are reported and discussed. The last part of the paper contains some conclusive considerations.
2. Governing equations
79
Consider a rectangular isotropic plate of dimension (a, b) and uniform thickness h, schematically represented, together with the local reference system adopted in this work, in Fig. 1. The plate is subjected to biaxial in-plane compressive load nx = −ξ Ph and n y = −ηPh, with (ξ, η) fixed ratio parameters, such that (ξ = η = 1) describes equibiaxial compression, (ξ = 1, η = 0) uniaxial compression in x direction, and so on. In plane stress conditions the constitutive behavior can be defined as follows:
81
⎡
⎤
80
⎡
⎤⎡
⎤
σ˙ x αxx αxy αxz ε˙ x ⎣ σ˙ y ⎦ = E ⎣ αxy α y y α yz ⎦ ⎣ ε˙ y ⎦ τ˙xy αxz α yz αzz γ˙xy
σ
σe = P ξ − ξ η + η
2
(3)
σe2
=
ρ
⎣
91 93 94 95 96 97 98 100 102 103 105 107 108 109 110
(4)
112
αi j can
114
113 115 116 117 118
c y y c zz − c 2yz
c xz c yz − c xy c zz c xx c zz − c 2xz
symm
119
⎤
0 0
⎦
c xx c y y − c 2xy
120
(5)
121 122 123 124
where:
c xx = 1 − 3 1 −
c yy = 1 − 3 1 − c xy = −
89
111
12
αxx αxy αxz ⎣ αxy α y y α yz ⎦ αxz α yz αzz 1
88
106
⎤
⎡
87
104
S i j S kl ε˙ kl
For both the FT and DT models the instantaneous moduli be expressed in the form:
⎡
86
101
In Eqs. (2) and (3) (G , λ, G s , λs ) are the elastic and the secant Lamè coefficients, G t is the tangent shear modulus, S i j are the stress deviator components and σe is the equivalent stress that, for the load conditions considered, assumes the value: 2
85
99
(2)
2 e
and the deformation theory of plasticity with the Hencky constitutive relation:
84
92
S i j S kl ε˙ kl
σ˙ i j = 2G s ε˙ i j + λs δi j ε˙ kk − 3(G s − G t )
83
90
(1)
where E is the elastic Young modulus and αi j are instantaneous moduli depending on the plasticity theory considered to model material behavior. Here two plasticity theories are considered, namely the incremental, or flow, theory of plasticity, with the Prandtl–Reuss constitutive equation:
σ˙ i j = 2G ε˙ i j + λδi j ε˙ kk − 3(G − G t )
82
1 2
Et E¯ Et E¯
s2y
4 s2x 4
1 − (1 − 2ν )
125
,
126
c xz = 0,
127 128
, Et E
c yz = 0,
−3 1−
Et E¯
129
sx s y 2
130
131
,
132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.3 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
3
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
Fig. 2. Ramberg–Osgood stress–strain relation.
34 35 36
c zz =
37 38 39 40 41 42 43 44 45 46 47 48
ρ
101
Et
− 1 σx sx = = −ξ ξ 2 − ξ η + η2 2 σe − 1 σy sy = = −η ξ 2 − ξ η + η2 2 σe
55
E¯ = E ;
G¯ = G =
E
58
in flow theory, whereas assume the secant values
E¯ = E s ;
G¯ = G s =
61 62 63 64 65 66
E 2(1 + ν ) + 3( EE − 1)
(9)
ε=
σe E
+k
σ0 σe E σ0
Et
= 1 + nk
σe σ0
σe =1+k Es σ0 E
n−1
(10)
has been adopted to model the stress–strain elastoplastic curve. In Eq. (10) σ0 represents the nominal yield stress, σe is the equivalent stress as defined in Eq. (4) and (k, n) are constant describing the shape of the stress–strain elastoplastic curve, as shown in Fig. 2. In such a hypothesis the ratio of the elastic modulus E and the
106
n−1
107
(11)
2 1
ε˙ y = s˙ y, y +
2
111 112 113 114 115 116 117 118 119 120 122 123
(12)
s˙ 2y ,x + s˙ 2z,x
110
121
The non-linear expressions for the effective strain components contributing to the strain energy can be put in the following form:
ε˙ x = s˙ x,x +
108 109
whose dependent on the (k, n) coefficients are reported, in terms of σe /σ0 ratio, in Fig. 3, where it is possible to observe how, for low values of σe /σ0 , both E s and E t assume the same value of the elastic constant E, then both decrease when σe /σ0 achieve a limit value, with a velocity depending on the assigned values of (k, n). The critical behavior of the plate has been analyzed within the framework of the Kirchhoff–Love kinematic model, hence the dis˙ (x, y , z, t ) = [˙sx s˙ y s˙ z ] T can be expressed placement field velocity U in terms of the three components U˙ = U˙ (x, y , t ), V˙ = V˙ (x, y , t ) ˙ (x, y , t ) of the middle surface as: ˙ =W and W
1
103 105
˙ s˙ z = W
n
102 104
˙ ,y s˙ y = V˙ − z · W
in deformation theory. In this paper the Ramberg–Osgood elastoplastic equation
E
˙ ,x s˙ x = U˙ − z · W
s
59 60
(8)
2(1 + ν )
56 57
(7)
and the constitutive terms ( E¯ , G¯ ), that are the elastic values
53 54
(6)
are constant depending on both the stress state and the incremental theory considered, through the values assumed by
51 52
corresponding secant and tangent values ( E s , E t ) can be obtained in the forms:
G¯ E = c zz c xx c y y − c 2xy Et
49 50
100
s˙ 2y , y + s˙ 2z, y
124 125 126 127 128 129 130
131 132
JID:AESCTE
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.4 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
4
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
Fig. 3. E t / E, and E s / E curves varying the Ramberg–Osgood coefficients n, k.
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
99
γ˙xy = s˙ x, y + s˙ y,x + s˙ y,x s˙ y, y + s˙ z,x s˙ z, y γ˙xz = s˙ x,z + s˙ z,x + s˙ y,x s˙ y,z + s˙ z,x s˙ z,z γ˙ yz = s˙ y,z + s˙ z, y + s˙ y, y s˙ y,z + s˙ z, y s˙ z,z
56
˙ )= δΩ(U˙ , V˙ , W
59
ΩU =
64 65 66
(14)
1 2
+
h/2
ΩV =
2
σx h
106
(17)
˙ ,xx − h2 /12 W ˙ ,xxy y + ξ Ph W ˙ , y y − h2 /12 W ˙ ,yyyy = 0 +ηPh W
107 108 109 110
111 112 113 114
(18)
115
˙ can be obtained. for the out-of-plane displacement W It is worth noting that, because in Eq. (13) the non-linear contribution of the displacement v˙ has been considered, both the Eqs. (17) and (18) are slightly different to the classical equation reported in literature, with new terms usually neglected under the von Kármán hypothesis.
116
(15)
(16)
Eqs. (17) and (18) can be analytically solved by means of the Kantorovich technique with which, through a separation of (x, y ) variables, a system of partial differential equations can be transformed in a set of ordinary equations for which a closed form solution exists. If simply supported, the boundary conditions on the edges x = ±b/2 are automatically satisfied by the following displacement field:
117 118 119 120 121 122
s2y ,x
+
s2z,x
+ σyh
s2y , y
+
s2z, y
124 125
−h/2
105
for the in-plane eigenmode components (U˙ , V˙ ), and:
αxx W˙ ,xxxx + 2(αxy + 2αzz ) W˙ ,xxy y + α y y W˙ , y y y y
102 104
− h(ξ P V˙ ,xx + η P V˙ , y y ) = 0
Eh3 /12
101
123
σ˙ x ε˙ x + σ˙ x ε˙ x + τ˙xy γ˙xy + τ˙xz γ˙xz + τ˙ yz γ˙ yz dz
1
Eh αxx U˙ ,xx + (αxy + αzz ) V˙ ,xy + αzz U˙ , y y = 0 Eh α y y V˙ , y y + (αxy + αzz )U˙ ,xy + αzz V˙ ,xx
100
103
3. The Kantorovich method
represents the strain energy functional for the plate:
62 63
δ(ΩU + Ω V )dΓ = 0
where:
60 61
Γ
57 58
(13)
In Eq. (13) and in the following, the comma before subscripts indicates partial derivative, e.g. s˙ i , j is equivalent to ∂ s˙ i /∂ x j . Unlike the classical von Kármán strain–displacement hypothesis, which considers only the non-linear contribution deriving from the out-ofplane displacement s˙ z , in Eq. (13) the additional non-linear terms related to the in-plane component s˙ y are also considered. Even if negligible for the buckling analysis of flat plate, such a contribution can considerably affect the critical response of stiffened plates and prismatic beams, when buckling mode involves flexo-torsional or global flexural buckling displacement (Figs. 4A, 4B), remaining negligible in case of local flexural mode (Fig. 4C). The equilibrium state is thus sought by invoking the principle of minimum potential energy:
54 55
Kirchhoff’s kinematical model represented in Eq. (12), the minimization of the deformation energy with respect to arbitrary vari˙ ) returns the following local equations: ation of (U˙ , V˙ , W
is the potential energy for the plate subjected to uniform in-plane compressive stress and δ is the variational symbol. Substituting the
126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.5 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
5
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
Fig. 4. Pictorial scheme of global flexural (A), flexo-torsional (B) and local flexural (C) buckling mode.
14
80
15 16 17 18 19 20 21 22
81
U˙ (x, y ) = u˙ ( y ) cos
mπ x
82
b
83
mπ x V˙ (x, y ) = v˙ ( y ) sin
˙ (x, y ) = w ˙ ( y ) sin W
84
b mπ x
85 86
(19)
b
87
˙ ( y ) are unknown functions, and ˙ =w where u˙ = u˙ ( y ), v˙ = v˙ ( y ), w m is the assigned number of half-waves in x-direction related to the buckling mode. Indicating with λ = b/mπ the half-wave length in buckling mode, in virtue of the displacement field (19), the system (17) becomes the following linear, ordinary differential equation system related to the in-plane displacements (u˙ , v˙ ):
88
95
30
λ2 αzz u˙ , y y + λ(αxy + αzz ) v˙ , y − αxx u˙ = 0
31
λ2 (α y y − ηPh/ A ) v˙ , y y − λ(αxy + αzz )u˙ , y
23 24 25 26 27 28 29
32 33 34
− (αzz − α Ph/ A ) v˙ = 0
36 37 38 39
˙ ,yyyy + (α y y − ηκ P ) w
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
92 93 94
98
(20)
η h − ξ κ λ2 P
˙ =0 ˙ , y y + λ4 αxx D − ξ hλ2 P w − 2λ2 D (αxy + 2αzz ) w
(21)
99
Fig. 5. (a) Assembled thin-walled plates, geometry and local reference system. (b) Example 1, Buckling load factor κ versus α /β aspect ratio for simply supported plates with three different thicknesses parameter ζ .
which admits a closed form solution in the form:
u˙ = λ11 c 1 e ρ1 y + λ12 c 2 e ρ2 y + λ13 c 3 e ρ3 y + λ14 c 4 e ρ4 y
T˙ (x, y ) = t˙( y ) cos
v˙ = λ21 c 1 e ρ1 y + λ22 c 2 e ρ2 y + λ23 c 3 e ρ3 y + λ24 c 4 e ρ4 y
˙ = δ11 d1 e ω1 y + δ12 d2 e ω2 y + δ13 d3 e ω3 y + δ14 d4 e ω4 y w
(22)
depending on the known constant values (λi j , δ1 j ) and (ρ j , ω j ), which can be derived from the eigenvectors and eigenvalues matrices associated with the coefficients of (20) and (21), and by the unknown vector x = [c i di ] T , whose components can be determined by imposing the proper in-plane and out-of-plane boundary conditions on the edges y = ±a/2. Such boundary conditions, if collected in a single, square matrix in the form:
L·x=0
(23)
allow the analysis of single, flat, rectangular plates through the eigenvalue analysis of the matrix L with which it is possible to derive, for the simplest case, also a closed form solution [8]. The stability analysis of thin-walled structures with an arbitrary cross-section, as reported in Fig. 5, require the additional definition ˙ , y and of a local stiffness matrix k such that, indicating with ϑ˙ = w ˙ = [u˙ v˙ w ˙ ϑ˙ ] y=±a/2 the nodal displacecollecting in the vectors u ˙ ] y=±a/2 the corresponding dual forces, ment, and in p = [t˙ n˙ p˙ m it is possible to describe the nodal force–displacement relationship in the form
p=k·u
(24)
100 101 102 103
Under the displacement hypothesis (19), the corresponding Kirchhoff forces can be expressed as:
65 66
91
97
40 41
90
96
˙ whereas Eq. (18) becomes an ordinary differential equation in w:
35
89
˙ (x, y ) = n˙ ( y ) sin N P˙ (x, y ) = p˙ ( y ) sin
107
b mπ x
108 109
b mπ x
110 111 112
b mπ x
113
(25)
b
n˙ = A
αzz u˙ , y + αzz λ−1 v˙
p˙ = D
117 118
119
α y y w˙ , y y y − λ−2 αxy + 4λ−2 αzz − ηPh w˙ , y
116
α y y v˙ , y − αxy λ−1 u˙
114 115
˙ ( y ) are: ˙ =m where t˙ = t˙( y ), n˙ = n˙ ( y ), p˙ = p˙ ( y ) and m t˙ = A
105 106
mπ x
˙ (x, y ) = m ˙ ( y ) sin M
104
−2
120
121 122
(26)
123
In Eq. (26) A = E · h/(1 − ν 2 ) and D = E · h3 /12(1 − ν 2 ) representing respectively the extensional and the flexural plate stiffness. ˙ ) reported Substituting in Eq. (26) the displacement field (u˙ , v˙ , w in Eq. (22) it is possible to express the relationship between the nodal force and corresponding displacement, in the form reported in Eq. (24) [14]. Finally, it is possible to obtain the global stiffness matrix K in a FEM-like procedure, and the related plastic buckling parameter
125
˙ =D m
α y y w˙ , y y − αxy λ
˙ w
124 126 127 128 129 130 131 132
JID:AESCTE
AID:3273 /FLA
6
[m5G; v1.149; Prn:9/03/2015; 15:29] P.6 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
Fig. 6. Example 1, Buckling load factor
33
κ versus ξ/η load ratio for simply supported plates with three different thicknesses parameter ζ .
99
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
100
P as the lowest eigenvalue of K, and the corresponding buckling mode as the related eigenvector. However, for the non-linearity of the matrix K( P ), both because P is contained through (λi j , δ1 j ) in the general solution (22) and because it is contained, through σe , in the elastoplastic constitutive coefficients ( E t , E s ), the eigenvalue solution requires a double, nested iterative strategy based on a trial and error procedure. Assuming the Shanley concept of continued loading, i.e. no unloading occurs at the instant of the plastic buckling, the outer iteration starts by assuming a value of P = P 0 small enough to not cause plastic deformation and evaluating the related instantaneous moduli αi j ( P 0 ). With such values the inner iteration solves the corresponding eigenproblem, in which the critical load P 1 represents the lower eigenvalue. If ( P 1 − P 0 )/ P 1 > err, where err is a prescribed error bound (in the proposed examples an err = 0.001 has been considered), the outer iteration restarts considering a new P 0 = (1 − η) P 0 + η P 1 (0 < η < 1). The value P i such that ( P i − P i −1 )/ P i < err is the required elastoplastic buckling load.
53 54 55
4. Examples
56 57 58 59 60 61 62 63 64 65 66
In the numerical analyses an Aluminum alloy AL 14S-T6, characterized by a Young modulus E = 72 400 MPa, a Poisson ratio ν = 0.32, a yield stress σ0 = 430 MPa and a value for the Ramberg–Osgood coefficients n = 10.9, k = 0.3485, has been considered. The results are reported in terms of the buckling coefficient κ = P cr (1 − ν 2 )/ζ E, in which ζ = π 2 h2 /12a2 is a thickness parameter. Here the three different values ζ = 0.0001, ζ = 0.001 and ζ = 0.002 (corresponding to a value of the thickness/length ratio h/a = 0.011, h/a = 0.0349 and h/a = 0.0493, respectively) have been considered.
4.1. Rectangular plates
101 102
In order to validate the numerical procedure related to the proposed model, a first analysis concerning the elastoplastic buckling of a flat, rectangular, simply supported plate of length a, width b, thickness h, subjected to uniaxial compression has been investigated, and compared with the analytical solution:
N cr =
π 2 D b 2 m2 b2
a2
αxx + 2(αxy + 2αzz ) +
a2 b 2 m2
αyy
N cr =
2
π D b m b2
a2
2
+2+
a
2
b 2 m2
104 105 106 107 108
(27)
derived in [5] as an elastoplastic generalization of the corresponding elastic solution: 2
103
109 110 111 112 113 114
(28)
reported, among others, in [15]. In both Eqs. (27)–(28) m represents the longitudinal half-wave number characterizing the buckling mode, whereas αi j reported in Eq. (27) are the elastoplastic coefficients defined in Eq. (1). Fig. 6 reports the buckling coefficient κ obtained varying the b/a aspect ratio related to the first 5 modes (m = 1, 2, . . . , 5) by means of both Eq. (27) and the proposed formulation, and for both the FT and the DT models (together with the ones derived from a linearly elastic model, here reported for comparison purpose). The results show an excellent agreement between analytical and numerical analysis. It is worth noticing that, because both the analyses derive from the same kinematical hypothesis on the sinusoidal behavior of the out-of-plane displacement w, the comparison only measures the accuracy of the stiffness matrix evaluation and of the double iterations required for solving the eigenvalue problem numerically. The figure also shows the coincidence between
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.7 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
the minimum values are κmin = 3.32 when a/b = (0.72, 1.45, 2.16, 2.88, 3.6) in FT and κmin = 2.67 when a/b = (0.82, 1.62, 2.49, 3.3, 4.1) in DT. Fig. 7 reports the critical response of the same plate subjected to biaxial loads obtained varying the load ratio ξ/η , compared with analytical ones proposed by Durban and Zuckerman [8]. Again the results are in good agreement, including the growing discrepancy between the results computed by the two elasto-plastic models, which give rise the plastic buckling paradox for which “bifurcation loads calculated using the simplest flow theory of plasticity consistently overestimated buckling loads of plates and shells obtained in tests. Calculations based on the less respectable deformation theory of plasticity gave reasonably good agreement with test results” Bijlaard, [6]. For a better focus of such paradox, Figs. 8 and 9 report the limit values σcrmin /σ0 and σcrmax /σ0 for all of the examples considered and their position on the Ramberg–Osgood curve – Fig. 8 – and on the corresponding ( E s / E , E t / E ) curves – Fig. 9 – adopted in the analysis, showing that if the minimum values coincide for any analysis performed, the differences on the maximum values are equal to zero when the critical load is contained in the elastic range of the stress–strain curve, which otherwise increase, reaching in the worst case, obtained with the thickness parameter ζ = 0.002 and the load parameters ξ = −η = 1, a ratio close to three (i.e. the FT analysis returns a critical load triple of the corresponding value obtained through the DT analysis). For a better comprehension of the cause involving such paradox is essential to undertake a closer analysis of the instantaneous moduli αxx , 2(αxy + 2αzz ) and α y y which appear in Eq. (18), and point out the differences on such coefficients derived from the use of the two theories. Substituting Eq. (6) in (5), and taking into account the values of (s x , s y ) reported in Eq. (7) the following expressions for such coefficients can be obtained:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
1 Et ξ 2 ξ 2 − ξ η + η2 αxx = 4−3 1− ρ E¯ 1 Et αyy = 4−3 1− η2 ξ 2 − ξ η + η2 ρ E¯ 1 Et Et αxy = 2 − 2(1 − 2ν ) −3 1− ηξ ξ 2 − ξ η + η2 ρ E E¯
34 35 36 37 38 39 40 41 42 43
αzz =
44
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 110
E
elastic and elastoplastic results for the thinner plate ζ = 0.0001, for which in any performed analysis the coefficient κ assumes a minimum value κmin = 4 when a/b = (1, 2, 3, 4, 5), and a maximum value κmax = 4.5 when a/b = 1.41. Again, the curve related to m number of half-waves intersects the curve related to m + 1 half-wave when a/b = (1.41, 2.45, 3.46, 4.47). When the plate becomes thicker (ζ = 0.001), elastic and elastoplastic curves diverge, even though they maintain the same general aspects. Both the FT and DT curves are characterized by a minimum value, that is κmin = 3.92 when a/b = (0.97, 1.9, 2.91, 3.89, 4.84) for the FT case and κmin = 3.87 when a/b = (0.98, 1.95, 2.92, 3.9, 4.87) in DT analysis. The intersections between the curves also change: they are a/b = (1.32, 2.35, 3.34, 4.32) and a/b = (1.34, 2.37, 3.36, 4.34) in FT and DT case, respectively. The differences between the two models increase for the thickest plate (ζ = 0.002). In such a case
whose values depend, through ( E t , E s ), on the actual load coefficient P and, through the (ξ, η) coefficients, on the load ratio existing between the assigned nx and n y in-plane loads. The latter affect the slope of αi j as shown in Figs. 10–12. Coherently with the elastoplastic theories adopted, it is reasonable to expect that coefficients decline as the load parameter P increases. Figs. 10–12 show how it occurs for any coefficient, excluding both αzz ( F T ), that assume the constant value G¯ / E = 2(1 + ν ), and αxy which, since it contains a term which depend on the product ηξ , can increase if such coefficients are characterized by different signs. As a consequence, if subject to a compressive force in one direction and a tensile force in the other, the instantaneous modulus 2(αxy + 2αzz ) related to the FT model presents an unnatural evolution, contrary to what happens in the DT model, in which the increasing values of αxy are canceled by the decreasing values of αzz , the latter depending on the secant value G¯ = G s reported
48
67
109
Fig. 7. Example 1, Position of maximum/minimum σcr /σ0 on the Ramberg–Osgood curve for any load ratio and thickness parameter considered.
47
51
G¯
111
E
46
50
Et ρ = + (1 − 2ν ) 3 − (1 − 2ν ) ¯G E 2 Et 2 −3 1− ηξ ξ − ξ η + η E¯
45
49
7
112 113 114
(29)
115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
8
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.8 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
Fig. 9. Example 1, Generalized constitutive variables
47
αi j for ξ = η = 1.
48 50 52 53
56 57 58 59 60 61 62 63 64 65 66
115 116
51
55
113 114
4.2. C-section columns
49
54
112
Fig. 8. Example 1, Position of maximum/minimum σcr /σ0 on the ( E t , E s )/ E curves for any load ratio and thickness parameter considered.
in Eq. (9). Such anomaly can be considered a concurrent cause of plastic buckling paradox. Fig. 13 shows how the differences between the two models decrease drastically by simply substituting G¯ = G s in the FT model. Finally, the results related to a single plate also show how the weight of the quadratic v-terms represented in Eq. (13) is negligible, consistently with von Kármán hypothesis for which, due to the buckling mode involving exclusively out-of-plane displacement, only the corresponding nonlinear strain related to the w-terms affect the buckling analysis of flat plates.
The next example analyzes the critical behavior of the thinwalled C-section columns represented in Fig. 14. As it is derived by a closed form solution, the proposed procedure does not depend on the discretization adopted in the analysis, and a model capable to determining the critical load related to any of the flexural, torsional or flexo-torsional mode, both local and global, can be obtained with the minimum number of nodes required to represent the section geometry that, for the proposed example, are the six nodes (A . . . F), again represented in Fig. 14. The related global stiffness matrix K considered in the eigenvector problem has then dimension (24 × 24). The reduced dimension of the stiffness matrix leads to a very fast numerical code, particularly suitable for parametric analysis. For instance, Fig. 15 shows the elastic buckling parameter κ obtained varying the length ratio b/a between 1 to 100, again considering the three thickness parameters ζ = 0.0001,
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.9 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
9
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32 33
Fig. 10. Example 1, Generalized constitutive variables
αi j for ξ = 0, η = 1.
98 99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 11. Example 1, Generalized constitutive variables
αi j when ξ = −1, η = 1.
132
JID:AESCTE
10
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.10 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35 36
101
Fig. 12. Example 1, Buckling factor
κ versus ξ/η load ratio considering G = G s in FT analyses.
37
103
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
102 104
ζ = 0.001 and ζ = 0.002. The numerical analysis has been performed considering both the classical von Kármán strain model, in the figures indicated with the apex c (classical) and the enriched model which takes into account the nonlinear contribute of the v-component, in the figure reported with the apex e (enriched). The results show the contrasting behavior for the three cases. When ζ = 0.0001, according with the length value, for short structures the buckling mode is web-triggered local (Fig. 16A) and involves only out-of-plane displacements (local-flexural mode), and a number of half-waves increasing – from 1 to 30 – with the length ratio. In such a phase the critical load remains substantially constant, equal to the minimum value related to the buckling curve related to the number of half-waves activated in buckling mode. For long structures the buckling mode becomes global, purely torsional (Fig. 16D) involving comparable out-ofplane and in-plane displacements. As a consequence in such phase the von Kármán simplification fails, and numerical analyses performed considering the non-linear strain term reported in Eq. (13) return a critical load sensibly inferior, and a lower value of the ratio b/alim characterizing the transition between local and global mode, that is equal to b/alim = 28 (with a critical mode characterized to 30 half-waves) under the von Kármán hypothesis, whereas it is equal to b/alim = 19 (and 22 half-waves) otherwise. For ζ = 0.001 and ζ = 0.002 the mode is flexural-torsional (Fig. 16B–C) also for short structures, and the differences between the elastic analyses performed with the two deformation models
are visible for any chosen ratio b/a. For long structures the mode becomes global, purely torsional in a first phase (Fig. 17D), purely flexural in a final phase (Fig. 17E), involving in any case comparable in-plane and out-of-plane displacement. When ζ = 0.0001 the elastoplastic results again coincide with those deriving from an elastic analysis. Fig. 17 shows the elastoplastic analyses performed for ζ = 0.001 and ζ = 0.002. Again the DT model, coupled with the enriched deformation model, returns the lowest critical load. However, the analyses performed considering the FT model follow the same critical path, and do not present the different paths observed on the single plate reported in Fig. 13. It is interesting to observe that for short structures the differences between the two elastoplastic models are more pronounced under the von Kármán hypothesis. Such differences tend to disappear for long structures, when the elastoplastic buckling load is closer to the elastic buckling load. Instead, in such a phase the differences generated by considering the contribute of v in Eq. (13) emerge.
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123
4.3. Close-sections columns
124 125
Finally, a buckling analysis of five different thin-walled closed section beams, whose form and dimensions are reported in Fig. 18, has been performed. The sections are labeled by the number n of the sides characterizing the section, and the dimensions are such that all the figures are inscribed in a circle of radius a, so that the beams fill the same space, whereas the thickness variation is such that the sections maintain the same area. The analyses are then
126 127 128 129 130 131 132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.11 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
11
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44 45 46
110
Fig. 13. Example 1, Buckling load factor κ versus α /β aspect ratio for simply supported plates with three different thickness parameter ζ .
111 112
47
113
48
114
49
115
50
116
51
Fig. 15. Example 2, linear-elastic buckling load factor κ versus α /β aspect ratio considering both the Classical von Kármán (c) and the enriched (e) strain–displacement model.
52 53 54
118 119 120
55 56 57 58 59 60 61 62 63 64 65 66
117
Fig. 14. Example 2, geometrical aspect ratio.
performed, again in terms of κ vs. b/a, on equivalent beams, in terms of overall dimensions and weight. Fig. 19 reports the buckling load obtained through elastic analyses. The results, reported in Fig. 19 considering linearly elastic behavior, show a general behavior very close to what was observed in the previous example, i.e. short plates characterized by a critical load almost constant and a local, purely flexural, critical mode (Fig. 20A), with a number of half-waves growing with the length whereas, when the length ratio b/a reaches a limit value b/alim , the critical mode becomes global, and the critical load-factor κ
121 122 123 124 125 126 127 128 129 130 131 132
JID:AESCTE
AID:3273 /FLA
12
[m5G; v1.149; Prn:9/03/2015; 15:29] P.12 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
Fig. 16. Example 2, local flexural (A), local flexo-torsional (B, C), global torsional (D) and global flexural (E) buckling mode.
15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
82 83
decreases for long plates. Because close sections show a greater torsional stiffness, purely flexural (Fig. 20B) behavior also characterizes the critical mode of the thickest plates, that do not present the flexo-torsional buckling load observed in the previous example. The limit value b/alim depends on the rigidity of the plate, and decreases with increasing the number n of the sides of the section. Fig. 19 also shows how for short plates, when the mode is local and the corresponding critical load is affected by the ratio between the length and the thickness of the single plate forming the section, the critical load-factor κ is greater for the triangle-shape section and it decreases as n increases. Vice versa, for long plates, when the mode becomes global and the corresponding critical load depends on the global stiffness of the section, the triangle-shaped section presents the lower critical factor, whereas the greatest is obtained for the square-shaped section. The elastoplastic behavior is also affected by the superior torsional stiffness, as reported in Fig. 21 for the thickest plates (ζ = 0.002). It shows a general behavior, in terms of critical path and associated critical mode, not so different to the one observed in elastic analysis. As in the previous examples, in all the performed analyses the DT model furnishes the lower critical coefficient, albeit the differences between the two models do not diverge, remaining almost constant varying the b/a ratio. It is worth noting that, increasing the number n, the gap between elastic and elastoplastic results increases, becoming almost double for n = 8.
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
44 45
110
5. Conclusion
111
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
81
112
Based on the Kirchhoff plate theory and making use of the Kantorovich technique, the present work provides a nonlinear model capable of deriving the critical behavior of prismatic beams and stiffened plates made of nonlinear, elastoplastic materials. A numerical procedure based on a closed form solution of the equilibrium equations has been derived and, since in such methods the results are not affected by numerical instabilities and/or meshsize, the approach can be adopted with any required accuracy. For validation purposes, the proposed formulation was first used to determine the critical behavior of rectangular plates. The comparison with exact solutions available in literature showed the correctness of the proposed procedure and, through a more accurate analysis of the elastoplastic coefficients derived by both DT and FT theories, it has shed a new light on the causes of the plastic buckling paradox. The analyses performed on a C-section and five different closed sections columns have shown the inadequacy of von Kármán hypothesis on the determination of both elastic and elastoplastic buckling loads when the buckling mode involves comparable in-plane and out-of-plane displacements and highlighted the need, in such cases, for a refined model. Finally, the sensitivity
113 114 115 116 117 118
Fig. 17. Example 2, elasto-plastic buckling load factor κ versus α /β aspect ratio considering both the Classical Von Kármán (c) and the enriched (e) strain-displacement model.
119 120 121 122
analyses have highlighted the possibility of following the transition from global to local buckling modes, and consequently the possibility of determining the corresponding variations in critical loads, by modifying key geometrical parameters of the structures and therefore the method can be helpfully exploited in design optimization procedures.
123 124 125 126 127 128 129
Conflict of interest statement
130 131
None declared.
132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.13 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
13
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10 11
76
Fig. 18. Example 3, geometrical aspect ratio.
77
12
78
13
79
14
80
15
81
16
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 19. Example 3, linear-elastic buckling load factor
κ versus α /β aspect ratio.
132
JID:AESCTE
14
AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.14 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
67
2
68
3
69
4
70
5
71
6
72
7
73
8
74
9
75
10
76
11
77
12
78
13
79
14
80
15 16
81
Fig. 20. Example 3, local and global buckling mode.
82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65 66
131
Fig. 21. Example 3, linear-elastic and elasto-plastic buckling load factor
κ versus α /β aspect ratio considering both FD and DT models.
132
JID:AESCTE AID:3273 /FLA
[m5G; v1.149; Prn:9/03/2015; 15:29] P.15 (1-15)
E. Ruocco / Aerospace Science and Technology ••• (••••) •••–•••
1
References
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
[1] Z.P. Bazant, L. Cedolin, Stability of Structures: Elastic, Inelastic, Fracture and Damage Theories, Courier Dover Editions, 1991. [2] C.M. Wang, C.Y. Wang, J.N. Reddy, Exact Solution for Buckling of Structural of Members, CRC Series in Computational Mechanics and Applied Analysis, 2009. [3] P. Weißgraeber, C. Mittelstedt, W. Becker, Buckling of composite panels: a criterion for optimum stiffener design, Aerosp. Sci. Technol. 16 (2012) 10–18. [4] Matthias Heitmann, Peter Horst, A new analysis model for the effective stiffness of stiffened metallic panels under combined compression and shear stress, Aerosp. Sci. Technol. 10 (2006) 316–326. [5] G.H. Handelmann, W. Prager, Plastic buckling of a rectangular plate under edge thrusts, naca-report-946, 1949, pp. 479–507, http://naca.central.cranfield.ac.uk/ reports/1949/naca-report-946.pdf. [6] P.P. Bijlaard, On the plastic buckling of plates, J. Aeronaut. Sci. 17 (1950) 485–493. [7] F. Bleich, Buckling Strength of Metal Structures, McGraw-Hill, New York, 1952. [8] D. Durban, Z. Zuckerman, Elastoplastic buckling of rectangular plates in biaxial compression/tension, Int. J. Mech. Sci. 41 (1999) 751–765.
15
[9] P. Le Grognec, A. Le Van, Some new analytical results for plastic buckling and initial post-buckling of plates and cylinders under uniform compression, ThinWalled Struct. 47 (2009) 879–889. [10] J.W. Hutchinson, Plastic buckling, Adv. Appl. Mech. 14 (1970) 67–144. [11] X.M. Wang, J.C. Huang, Elastic/plastic buckling analyses of rectangular plates under biaxial loadings by the differential quadrature method, Thin-Walled Struct. 47 (2009) 879–889. [12] W. Zhang, X.W. Wang, Elastic/plastic buckling analysis of thick rectangular plates by using the differential quadrature method, Comput. Math. Appl. 61 (2011) 44–61. [13] M. Kadkhodayan, M. Maarefdoust, Elastic/plastic buckling of isotropic thin plates subjected to uniform and linearly varying in-plane loading using incremental and deformation theories, Aerosp. Sci. Technol. 32 (2014) 66–83. [14] E. Ruocco, M. Fraldi, Critical behavior of flat and stiffened shell structures through different kinematical models: a comparative investigation, Thin-Walled Struct. 60 (2012) 205–215. [15] S.P. Timoshenko, J.M. Gere, Theory of Elastic Stability, McGraw–Hill, New York, 1961.
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
17
83
18
84
19
85
20
86
21
87
22
88
23
89
24
90
25
91
26
92
27
93
28
94
29
95
30
96
31
97
32
98
33
99
34
100
35
101
36
102
37
103
38
104
39
105
40
106
41
107
42
108
43
109
44
110
45
111
46
112
47
113
48
114
49
115
50
116
51
117
52
118
53
119
54
120
55
121
56
122
57
123
58
124
59
125
60
126
61
127
62
128
63
129
64
130
65
131
66
132