Computing the structural buckling limit load by using dynamic relaxation method

International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 2 3 4 5 6 7 8 9 10 11 12 13 Q2 14 15 Q1 16 17 Q3 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

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Computing the structural buckling limit load by using dynamic relaxation method Mohammad Rezaiee-Pajand n, Hossein Estiri Civil Engineering Department, Ferdowsi University of Mashhad, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 11 August 2015 Received in revised form 22 January 2016 Accepted 31 January 2016

The numerical structural analysis schemes are extensively developed by progress of modern computer Q4 processing power. One of these approximate approaches is called "dynamic relaxation (DR) method." This technique explicitly solves the simultaneous system of equations. For analyzing the static structures, the DR strategy transfers the governing equations to the dynamic space. By adding the fictitious damping and mass to the static equilibrium equations, the corresponding artificial dynamic system is achieved. The static equilibrium path is required in order to investigate the structural stability behavior. This path shows the relationship between the loads and the displacements. In this way, the critical points and buckling loads of the non-linear structures can be obtained. The corresponding load to the first limit point is known as buckling limit load. For estimating the buckling load, the variable load factor is used in the DR process. A new procedure for finding the load factor is presented by imposing the work increment of the external forces to zero. The proposed formula only requires the fictitious parameters of the DR scheme. To prove the efficiency and robustness of the suggested algorithm, various geometric non-linear analyses are performed. The obtained results demonstrate that the new method can successfully estimate the buckling limit load of structures. & 2016 Published by Elsevier Ltd.

Keywords: Dynamic relaxation Buckling limit load Load factor Equilibrium path Fictitious parameters

1. Introduction The structural stability can be evaluated by estimating its critical load. In general, an eigenvalue problem should be solved to calculate this load [1]. In this approach, the aforesaid load makes the stiffness matrix singular, in which the related determinant is zero. It is worth emphasizing that this technique is unsuitable for large structures since it consumes a great deal of time [2]. Alternatively, it is possible to use the equilibrium path for computing the critical load. In this strategy, the critical load is the maximum load before reaching the first limit load or limit displacement points [3–6]. Based on this, three types of load–displacement curves are resulted from the structural responses. In the first one, no buckling occurs, whereas the second one, in which no bifurcation point emerges, has the limit buckling point. The third curve includes the bifurcation point in addition to the limit point. In the last case, the buckling may occur before or after the bifurcation points [5]. It can be emphasized that the bifurcation points generally appear in structures with initial defects [6]. So far, various solution techniques have been used for post-buckling analysis of structures. An improved arc-length method with high-efficiency n

Corresponding author. Tel./fax: þ 98 51 38412912. E-mail address: [email protected] (M. Rezaiee-Pajand).

was proposed by Zhu and Chu [7]. They utilized this procedure to analyze the composite structures. Lee et al. applied FSDT and the von-Kármán non-linear strain–displacement relationship for postbuckling analysis of the FGM plates under edge compression and thermal loads [8]. Asemi and Shariyat investigated effects of the materials auxetic on uniaxial and biaxial post-buckling behaviors of the FGM plates [9]. Uniaxial and biaxial post-buckling behaviors of the longitudinally graded plates were studied by Shariyat and Asemi [10]. The DR method is an approximate numerical procedure. For the first time, this scheme was introduced by Day in 1965 [11]. In this algorithm, a static problem is converted to a fictitious dynamic one. In other words, the mass and damping are added to the static relations for forming an artificial dynamic system [12]. This method is based on the second-order Richardson's approach which is developed by Frankel [13]. Then, Rushton applied this technique to non-linear problems. This researcher used the DR strategy to perform the geometric non-linear analysis of bending plates [14]. It should be noted that the mass and damping matrices play important roles in preserving the stability of the DR process. Furthermore, the time step has also influenced the convergence rate. Researchers presented various procedures for estimating the fictitious parameters of the DR scheme. In this paper, the DR process is employed in tracing the structural equilibrium path and computing the buckling load. At the

http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022 0020-7462/& 2016 Published by Elsevier Ltd.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 68 Q6 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

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

2

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 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

first stage, the formulas of the DR procedure are presented. Afterwards, the best schemes for calculating the fictitious parameters are reviewed. Then, a new approach is suggested for finding the load factor. To achieve this goal, the work increment of external forces is written in terms of the DR artificial parameters. By setting this increment to zero, a new formula is obtained for determining the load factor. It should be reminded that the work increment is dependent on the displacement and load factor increments. Moreover, the displacement increment is related to velocity. Furthermore, velocity depends on the fictitious mass and damping. As a result, the proposed relation is only associated with the parameters of the DR scheme. This factor is employed in estimating the structural buckling limit load. It should be noted that the main aim of this paper is the calculation of buckling limit load. This load is estimated by the equilibrium path of structures. In other words, at first, the equilibrium path of structures is plotted. Based on this path, the corresponding load to the first limit point is known as buckling limit load.

successive iterative steps. In this way, new formulas were presented for calculating the artificial mass and damping [25]. Rezaiee-Pajand and Sarafrazi utilized the power iteration process for finding the minimum eigenvalue [26]. By minimizing the error between two successive iterations, Rezaiee-Pajand et al. proposed another approach for computing the damping matrix [27]. In another study, Rezaiee-Pajand and Sarafrazi set the damping parameter to zero, and proposed a formula for the time step ratio [28]. Alamatian recommended a new formula for estimating the fictitious mass of the kinetic DR technique [29]. Rezaiee-Pajand et al. found a new time step by minimizing the unbalanced energy function of the fictitious dynamic system, [30]. Finally, RezaieePajand and Rezaiee presented a new time step for the kinetic DR strategy [31]. In most of the DR approaches, the artificial damping is estimated by the Rayleigh's principle [20]. This relationship is demonstrated as follows: sffiffiffiffiffiffiffiffiffiffiffiffiffiffi XT F ð3Þ C ¼ 2M XT MX

2. The dynamic relaxation method

In this relation, the stiffness matrix and the internal force vector are shown by S and F, respectively. The absolute values of the stiffness matrix entries in a row are calculated and summed to find the mass of corresponding degree of freedom, as follows [12]:

The basic relations of the DR technique are achieved by using the central finite difference formulation. These equations are presented in the below form: : n þ 12

Xi

2mnii C nii t n : n  12 2t n n ¼ X þ r n ; r n ¼ pni  f i ; i ¼ 1; 2; :::; ndof 2mnii þC nii t n i 2mnii þ C nii t n i i : n þ 12

X ni þ 1 ¼ X ni þ t n þ 1 X i

; i ¼ 1; 2; :::; ndof

ð1Þ

t 2 X   S 4 j ¼ 1 ij ndof

mii ¼

ð4Þ

ð2Þ

In these relations, the i-th entry of the fictitious mass and damping matrices, the artificial time step, the i-th entry of the n internal and residual force vectors are mnii ,C nii ,t n ,f i and r ni , correspondingly. It should be noted that superscript n denotes the n-th iteration of the DR strategy. Additionally, the external force of the static system and the number of degrees of freedom are indicated by pni and ndof, respectively. Moreover, the displacement and _ correspondingly. Noticevelocity vectors are shown by X and X, ably, the DR process is explicit. In other words, only vector operators are required in the solution process. This characteristic is the result of the diagonal mass and damping matrices. It should be highlighted that Eqs. (1) and (2) are successively used to achieve an acceptable error. To estimate the fictitious parameters, various procedures are suggested. Brew and Brotton suggested a method in which the mass of each degree of freedom was assumed to be proportional to its diagonal element in the structural stiffness matrix [15]. Bunce computed the critical damping by utilizing the Rayleigh's principle [16]. By deploying Gershgorin's circle theory, Cassell and Hobbs calculated the fictitious mass [17]. Papadrakakis presumed that both mass and damping are proportional to the stiffness matrix diagonal entries [18]. Underwood suggested a well-known formulation for iterations of the DR method [12]. Besides, by using the Rayleigh's principle, the damping factor and the time step is calculated by Qiang [19]. Zhang and Yu proposed a new formula for estimating the artificial damping, which is based on the Rayleigh's principle [20]. Zhang et al. used nodal damping, in which the degrees of freedom for each node have similar damping [21]. Rezaiee-Pajand et al. and Kadkhodayan et al. minimized the residual force and proposed a new relation for calculating the time step [22,23]. Another well-known formulation is called the kinetic DR method. In this scheme, the damping is ignored. In other words, the damping matrix is set zero in Eq. (1) [24]. Rezaiee-Pajand and Alamatian minimized the displacement errors of the two

3. The proposed formulation Non-linear behavior of structures has various characteristics corresponding to load and displacement limit points in stable and unstable paths, buckling points and post-buckling behavior. Sometimes, the increase in the load and displacement simultaneously occurs in a stable path. Besides, the reduction in the load value and increase in the displacement may be observed in an unstable path. The existence of limit points in the structural behavior path leads to the analysis complexity. One application of the DR algorithm is tracing the structural equilibrium paths, and consequently estimating the buckling limit load. In common formulations of the DR method, the external load is constant in each loading increment. Therefore, this technique cannot compute the critical load of structures. Nevertheless, several automatic approaches for removing this difficulty have been proposed [3,32– 34]. In these strategies, variable external force is deployed in the DR process. The load factor is usually utilized for changing the load level. If the load factor is shown by λ, the following residual force can be obtained: R ¼ λP F

ð5Þ

In the last relation, the residual force and reference load vectors are shown by R and P, respectively. The reference load is also the external force. Fig. 1 illustrates the equilibrium path of a structure with a snap-through point. Firstly, the load and displacement increase simultaneously on the stable path. This increase proceeds until reaching the first limit point. In Fig. 1, a load limit point exists. Afterwards, the load decreases whereas the displacement increases. After limit load, the structure becomes unstable [6]. The part of this equilibrium path before the limit load is called the prebuckling region, and the section after this point is the postbuckling part. It is worth emphasizing that the load before the first load or displacement limit point denotes the buckling limit load of the structure [3–6]. This load is shown in Fig. 1 by Pcr [4,35].

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

In this study, a new formula for the load factor is achieved by setting the work increment of the external force to zero. The work of these forces can be calculated by employing Eq. (6). Note that Eq. (7) is the work increment. W ¼ XT ðλPÞ

wn ¼

)

ndof X

λ n pi X i

ð6Þ

i¼1

δW ¼ δXT ðδλPÞ

)

δ wn ¼

ndof X

δλn pi δX i

i¼1

)

δ wn ¼

ndof X

: n  12

δλn pi t n þ 1 X i

ð7Þ

i¼1

In these relations, λ, W, δw and δX are the load factor, the work of external forces, the work increment and the displacement increment, respectively. It is intended that the load factor increment, δλ, is calculated so that the work increment of external forces, Eq. (7), is zero. The residual force of the n-th step can be calculated by deploying Eq. (8). By inserting the velocity, obtained from Eq. (1), into relationship Eq. (7) and using equality Eq. (8), the next formula can be found: r ni ¼ r ni  1 þ δλ pi n

δwn ¼ t n þ 1

ndof X

ð8Þ "

2mnii C nii t n : n  12 2t n n pi X þ ðr n  1 þ 2mnii þC nii t n i 2mnii þ C nii t n i

δλ

i¼1



δλn ¼

ndof P i¼1

h i n1 pi 2t n r ni  1 þ ð2mnii  t n C nii ÞX : i 2 2mn þ t n C n ii

ii

2t n

ndof P i¼1

Step 2. - The tangent stiffness matrix and internal force vector are found. Step 3. - The mass and damping matrices are calculated.

ð9Þ

Step 5. - Eq. (5) is employed in the calculation of the residual force vector.  T    Step 6. - If RP T RP  o 10  4 , go to Step 8. Otherwise, the velocities

δλn pi Þ

ð10Þ

ð11Þ

Herein, the suggested load factor is multiplied by the reference force to calculate the load applied to the structure. By utilizing this load, the structure is analyzed, and the responses are found. The next increment will start using these results. It is worth

Pcr

Step 1. - The initial value of the velocity and displacement should be zero.

Step 4. - The load factor should be computed by using Eqs. (10) and (11).

ðpi Þ2 2mnii þ t n C nii

λn þ 1 ¼ λn þ δλn

emphasizing that the arc length methods can find both the static path and buckling limit load. However, in the arc length approach, both the number of iterations and convergence points depend on the magnitude of the arc length. In other words, the structural analyst should assume the magnitude of the initial arc length at the beginning of the analysis. As a result, a slight change in the arc length may result in divergence or incorrect buckling load. On the other hand, the proposed technique is only related to DR parameters. The suggested formula can be calculated by the mass, damping, external load and residual force; as Eq. (11) indicates it. Consequently, the analysis simply requires vector operators, whereas the analyzer role in the proposed method is much less than the common arc length technique. Based on the authors' formulation, the steps required for calculating the buckling loads of the structures are summarized in the following lines:

#

By setting Eq. (9) zero, two values for the load factor can be achieved. It can be proved that one of them is zero. The other one can be obtained from Eq. (10). Since the critical load is dependent on λ, the first value is not acceptable. In other words, setting δλ to zero leads to a constant load level, which is used in the common DR procedure. Finally, the load factor can be computed by the following equality:

should be updated by deploying Eq. (1). Step 7. - The nodal displacement should be updated by applying Eq. (2), and the analysis process continues from Step 2. Step 8. - The displacements and load factor of the current increment are printed. Step 9. - If the target load or the target displacement is found; the process is finished. Otherwise, λ ¼ λ þ Δλ and the analysis should continue from Step 2. It is important to note that each of the methods introduced in Section 2 can be deployed in Step 3 for estimating the fictitious mass and damping parameters. As it was mentioned, Δλ in Step 9 is the load factor increment, which is employed at the beginning of the next step. By using smaller values of this parameter, more accurate results will be obtained. Furthermore, the target load is

Upper limit point (buckling load)

Stiffening

Softening Lower limit point

Load

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 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

3

Pre-buckling region

Post-buckling region

Stable path

Unstable path

Stable path

Displacement

Fig. 1. Estimation of the buckling limit load by using DR scheme.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

Fig. 2. The arch truss geometry. Table 1 Nodal coordinates of the arch truss. Node

1 and 19

2 and 18

3 and 17

4 and 16

5 and 15

6 and 14

7 and 13

8 and 12

9 and 11

10

X (cm) Y (cm)

7 3429 0

7 3048 50.65

72667 34.75

7 2286 83.82

7 1905 65.3

7 1524 110.85

71143 87.99

7762 128.5

7 381 100.65

0 134.6

Table 2 The cross-section area of the members of the arch truss. Member

1–10 and 35

11 and 12

13–16

17 and 18

19–22

23 and 24

25 and 26

27 and 28

29–32

33 and 34

Cross-section areas (cm2)

51.61

64.52

83.87

96.77

103.23

161.29

193.55

258.06

290.32

309.68

300 250 200

Load(kN)

150 100 50 0 -50 -100 -150 0

25

50

75

100

125

150

175

200

225

250

275

300

2.5

2.75

Tip Deflection(cm)

Tip deflection

300 250 200 150 Load(kN)

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 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

100 50 0 -50 -100 -150

0

0.25

0.5

0.75

1

1.25

1.5

1.75

2

2.25

Horizontal Displacement- Node 12(cm)

Horizontal displacement of node 12 Fig. 3. The arch truss equilibrium path. (a) Tip deflection. (b) Horizontal displacement of node 12.

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

X-Y view

X-Z view

Fig. 4. The ceiling space truss geometry. (a) X–Y view, (b) X–Z view.

8000 6000 4000

Load(kN)

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 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

5

2000 0 Horizontal displacement

-2000 -4000

Vertical deflection

-6000 -8000 -70

-60

-50

-40

-30

-20

-10

0

10

20

30

40

50

60

Displacement(mm)

Fig. 5. The ceiling space truss equilibrium path for node A.

X-Z view

X-Y view Fig. 6. The star truss. (a) X–Y view (b) X–Z view.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

6

the maximum force before the first limit point occurs. This load is the buckling limit load of the structure [3].

4.1. The arch truss In Fig. 2, the geometry of this 2D truss is demonstrated [36]. This structure includes 19 nodes and 35 elements. Node numbers are shown by the circles in Fig. 2. The nodal coordinates are inserted in Table 1. Moreover, the cross- section areas of the members are listed in Table 2. Furthermore, the elasticity modulus is 64.964 GPa. The references load P and maximum λ are 30 kN and 10, respectively. Therefore, the target load is achieved 300 kN. For finding the buckling limit load, it is presumed that Δλ ¼0.04. Fig. 3 portrays the arch truss equilibrium path. According to the structural equilibrium path, the load and deflection increase simultaneously. This process is continued until the deflection is equal to 38.44 cm. The corresponding load is 25011.894 kN, which shows the first snap-through point. Afterward, the structure enters the unstable region. At this state, the load reverses and deflection increases. The next limit point occurs in the load and deflection of  106.044 kN and 225.96 cm, respectively. Then, the load and deflection increase simultaneously again. It should be reminded; the horizontal displacement decrease after the second limit point. Based on the proposed technique, the buckling load is 25011.894 N.

4. Numerical samples In this section, the buckling limit loads of several frames, trusses and shells are calculated by deploying the proposed approach. A FORTRAN program is developed to achieve this aim. Recall that the aforesaid structures have geometric non-linear behavior. The suggested formulation is not related to the structural type. Hence, the analysis procedure can be conducted by utilizing any finite elements. For tracing the equilibrium path, Δλ is assumed to be one. Nevertheless, smaller values are applied in finding the buckling limit load. Table 3 The applied load cases of the star truss. Load case

Loading node number

Load value

Symmetric load [3,33]

1

P

4 3

Load(N)

2 1 0 -1 -2 -3

0

5

10

15

20

25

30

35

40

45

50

Tip Deflection(mm)

Tip deflection

4 3 2

Load(N)

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 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

1 0 -1 -2 -3 -0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Horizontal Displacement(mm)

Horizontal displacement of node 2 Fig. 7. The equilibrium path of the star truss. (a) Tip deflection (b) Horizontal displacement of node 2.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4.2. The ceiling space truss A three-dimensional truss is shown in Fig. 4. Previously, other researchers analyzed this structure [37]. The truss height is 9.17 m.

Moreover; the young's modulus and members' cross-section areas are 210 GPa and 0.01 m2, correspondingly. The references and target loads are 900 kN and 9000 kN, respectively. According to the suggested algorithm, the buckling load of this structure is

Fig. 8. The geodesic truss geometry.

10 8

Load(kN)

6 4 2 0 -2 -4

0

10

20

30

40

50

60

70

80

90

100

Tip Deflection(mm) Tip deflection

10 8 6 Load(kN)

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 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

7

4 2 0 -2 -4

-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

Horizontal Displacement- Node A(mm) Horizontal displacement of node A

Fig. 9. The geodesic truss equilibrium path. (a) Tip deflection (b) Horizontal displacement of node A.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

8

X-Yview

X-Zview

Fig. 10. The Schwedler's dome geometry. (a) X–Y view (b) X–Z view.

600 500 400 300 Load(N)

200 100 0 -100 -200 -300 -400

0

50

100

150

200

250

300

350

400

450

500

550

600

650

700

750

0.5

1

1.5

2

2.5

3

Deflection- Node 1(mm) Tip deflection

600 500 400 300 200 Load(N)

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 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

100 0 -100 -200 -300 -400 -4.5

-4

-3.5

-3

-2.5

-2

-1.5

-1

-0.5

0

Horizontal Displacement- Node 2(mm) Horizontal displacement of node 2

Fig. 11. The equilibrium path of the Schwedler's truss. (a) Tip deflection (b) Horizontal displacement of node 2.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

5050.306552 kN. It is worth emphasizing that the magnitude of the aforesaid load is compatible with other researchers' load [37]. For finding the buckling limit load, it is assumed that Δλ ¼ 0.04. The load-horizontal and vertical displacement curves are shown in Fig. 5. Note that these curves are similar to other researchers' curves [37]. 4.3. The star truss In this part, the star truss, shown in Fig. 6, is analyzed under symmetric load case. In previous research works, this truss was solved. Based on these studies, its buckling load is 3.1370844 N [3,33]. This 3D truss has 13 nodes, 24 members and 21 degrees of freedom. The moduli of elasticity and cross-section areas of members are 5 GPa and 2 mm2, correspondingly. In Table 3, the applied load patterns are introduced. The target load and maximum λ are 4 N and 10, respectively. For this truss, Δλ is 0.01 for finding the buckling limit load. The equilibrium path is illustrated in Fig. 7. Two limit load points exist in the path. The load factor at the first point is 3.155936 N, and the relevant displacement is 7.71 mm. At this limit point, the structure becomes unstable. The related buckling load is 3.155936 N for the symmetric load. Afterwards, the load increment decreases until another limit point is reached. The corresponding load factor and displacement at this point are  2.76064 N and 30.3 mm, correspondingly. At this state, the load reverses and begins to increase. The snap-back point occurs when the horizontal displacement is 0.4 mm.

members are 68950 N/mm2 and 6.5 cm2, correspondingly. The reference's load and maximum λ are 1000 N and 10, respectively. In this sample, the load factor increment is 0.04 for finding the buckling limit load. The boundary nodes of this truss are simply supported. To find the structural shape, Eq. (12) is utilized. In this relation, the units are in terms of the meter. The equilibrium path of the tip node and node A are demonstrated in Fig. 9. By employing the proposed method, it is concluded that the buckling load of this truss is 3170.887 N. x2 þ y2 þ ðz þ 7:2Þ2 ¼ 60:48

ð12Þ

4.5. The Schwedler's truss This structure is shown in Fig. 10. The Schwedler's dome is a truss with 240 members and 219 degrees of freedom [38]. This truss is a very large structure. The modulus of elasticity and the cross-section areas of members are 210 GPa and 4.5 cm2, correspondingly. The supports placed at the structural lowest level are simple. The references and target loads are 600 60 N and N, respectively. Fig. 11 shows the load–displacement curves of the tip node and node 2. The load and displacement of the first load limit points are 453.2 N and 107 mm, respectively. Moreover, the load and displacement of the second limit point are  356 N and 457 mm, correspondingly. The snap-back point related to the horizontal displacement occurs when the displacement reaches 4.336 mm. According to the suggested approach, the buckling load is 453.24924 N. It is worthwhile to mention that Δλ ¼ 0.04 is assumed in this sample for finding the buckling limit load.

4.4. The geodesic truss 4.6. The Toggle's frame This section is devoted to analyze the geodesic truss shown in Fig 8 [37]. The elasticity modulus and the cross-section area of the

Fig. 12. The Toggle's frame.

In what follows, the efficiency of the proposed algorithm in analyzing the frames is investigated. For this purpose, the developed program is verified by calculating the buckling load of the frame demonstrated in Fig. 12. This structure was analyzed previously [32,34,39,40]. In this paper, five elements are used to model each member of this frame. The young's modulus, the members' cross-section area and the moment of inertia are 71 GPa, 118.06 mm2 and 374.61 mm4. The maximum λ and load are 10 and 450 N, respectively. Furthermore, the load factor increment is 1 for finding the buckling limit load. By using the suggested strategy, it is concluded that the buckling load of this frame is 155.44414296 N. It is worth emphasizing that initial

450 400 350 300 Load(N)

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 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

9

250 200 150 100 50 0

0

2

4

6

8

10 12 Vertical Deflection(mm)

14

16

18

20

Fig. 13. The load-tip deflection curve for the Toggle's frame.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

10

4.7. The cylindrical roof

displacements are not required for analysis of this frame, in contrast to the columns. Fig. 13 demonstrates the equilibrium path corresponding to the tip node. Note that the obtained results are similar to the results of previous researches [32,34,39,40].

In what follows, the buckling loads of several shells are estimated. In these samples, the load factor increment is 1 for finding

Fig. 14. The cylindrical roof.

3 Node B 2.5 Node A

Load(kN)

2

1.5

1

0.5

0

0

3

6

9

12

15

18

21

24

27

30

Deflection(mm) Fig. 15. Load–deflection curves of the cylindrical roof (t¼ 12.7 mm).

1 Node B 0.8 Node A 0.6 Load(kN)

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 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

0.4 0.2 0 -0.2 -0.4

0

3

6

9

12

15

18

21

24

27

30

33

Deflection(mm) Fig. 16. Load–deflection curves of the cylindrical roof (t¼ 6.35 mm).

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

Fig. 17. The shallow spherical cap.

80 70 60 50 Load(kN)

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 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

11

40 30 20 10 0

0

25

50

75

100

125

150

175

200

225

250

275

300

325

Deflection(mm)

Fig. 18. The equilibrium path of the shallow spherical cap.

Fig. 19. The semi-cylindrical shell geometry.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

12

1500 Vertical Deflection 1250 Horizontal Displacement

Load

1000

UY

750 WZ 500

250

0 -10

0

10

20

30

40

50

60

70

80

90

100

110

120

130

140

150

Displacement Node A

1500

1250

1000 Load

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 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

UY

750 UX 500

250

0

0

2.5

5

7.5

10

12.5

15

17.5

20

22.5

25

Displacement Node B Fig. 20. The load–displacement curves for the semi-cylindrical shell. (a) Node A (b) Node B.

the buckling limit load. At first, the cylindrical roof shown in Fig. 14 is investigated. The edges in Y direction are simply supported, and other boundaries are free. Due to symmetry, only onequarter of the roof is modeled with a 10  10 mesh. Furthermore, the elasticity modulus and the Poisson's ratio are 310275 GPa and 0.3, respectively. Two different thickness are considered for this structure, 12.7 mm and 6.35 mm [33,41]. For these thicknesses, when the suggested scheme is used, the critical loads are 2.228628 kN and 0.589068 kN, correspondingly. It should be noted that the target load is 3 kN and 1 kN for the thin and thick thicknesses, respectively. Moreover, maximum λ is 10 for both. The static equilibrium paths of these structures are shown in Fig. 15 and Fig. 16. Based on Fig. 16, a local jump occurs in the load–displacement curve. It is worthwhile to mention that the critical loads are correctly calculated even with the aforesaid jump. 4.8. The shallow spherical cap In Fig. 17, the geometry of this shell is illustrated. Due to the symmetry, only one-quarter of this structure is analyzed by using 200 elements. The four edges of this shell are simply supported. Moreover, the thickness, the young's modulus and the Poisson's ratio are 99.45 mm, 68.944 GPa and 0.3 [42]. The maximum load and lambda are 80 kN and 10, respectively. For finding the buckling limit load, the load factor increment is 1. Based on proposed

Fig. 21. The deformed shape of the semi-cylindrical shell.

method, the critical load of this structure is 49912.60928 N. Furthermore, the static equilibrium path is shown in Fig. 18. 4.9. The semi-cylindrical shell At this stage, the benchmark structure shown in Fig. 19 is analyzed. In edges parallel with Y direction, all the nodes for the Z translations and also rotations around axis Y are restrained. Additionally, node A is placed in the middle of the free edge. The

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

modulus of elasticity, the thickness and the Poisson's ratio are 2.0685  107, 0.03 and 0.3, respectively [43,44]. These values are dimensionless. Due to the symmetry, only half of the structure is modeled by a 20  12 mesh. The reference load and maximum λ

are 150 and 10, correspondingly. In these samples, the load factor increment is 1 for finding the buckling limit load. In Fig. 20, the equilibrium paths, corresponding to the vertical and horizontal displacements of nodes A and B, are shown. Accordingly, it is

Fig. 22. Two-bay cylindrical roof.

100 Node A

90

Node B 80

Load(kN)

70 60 50 40 30 20 10 0

0

3

6

9

12

15

18

21

24

27

30

Deflection(mm) Tip deflection

100 90 80 70 Load(kN)

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 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

13

60 50 40 30 20 10 0 -0.4

-0.35

-0.3

-0.25

-0.2

-0.15

-0.1 -0.05 0 0.05 Horizontal Displacement(mm)

0.1

0.15

0.2

0.25

0.3

Horizontal displacement of node B Fig. 23. The equilibrium path for the two-bay roof. (a) Tip deflection (b) Horizontal displacement of node B.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

14

obvious that this structure has a snap-back point in the equilibrium path of the node B. The load corresponding to this point is the buckling limit load. This load is 4410.856. Fig. 21 demonstrates the shape of the structure under the maximum load. 4.10. The two-bay cylindrical roof In this sample, the shell demonstrated in Fig. 22 is analyzed. This structure has not been solved previously. Due to symmetry, only one-quarter of this shell is modeled by using a 15  10 mesh. The elasticity modulus, the Poisson's ratio and the thickness are 20 GPa, 0.3 and 3 cm, correspondingly. Four corner nodes, C, D, E and F, are simply supported. Nodes A and B are placed in the

middle of the edges. The target load and maximum λ are 100 kN and 10, correspondingly. The load factor increment is 1 for finding the buckling limit load. In the static equilibrium path of this structure, a displacement limit point emerges at first. This curve is portrayed in Fig. 23. The load corresponding to this point is 57.5822 kN. Then, the load limit point appears. The load related to this point is 70.9 kN. The final deformed shape of this shell is illustrated in Fig. 24. 4.11. The dome shell The dome shell shown in Fig. 25 is analyzed in this example. The young's modulus, the Poisson's ratio and the thickness are 2  1011 N/cm2, 0.3 and 0.3 mm, respectively. This structure is analyzed with 360 elements. The boundaries of this shell are simply supported. The shell widths in X and Y directions are 55 cm and 35 cm, correspondingly. Moreover, the structure's height is 10 cm [45]. The load is applied to node A. The reference and maximum loads are 30 kN and 300 kN, correspondingly. In this example, the load factor increment is 1 for finding the buckling limit load. The load–displacement curve of node A is depicted in Fig. 26. In these equilibrium paths, two load limit points are existed. The corresponding loads are 159.2899 kN and 45.5826 kN, respectively. Additionally, the deflections are 0.851 mm and 2.34 mm for point A, correspondingly. According to the suggested procedure, the buckling limit load is 159.2899 kN. 4.12. The second-order equation

Fig. 24. The deformed shape of the two-bay roof.

The final example is the simplest mathematical model. Here, the following equation is solved [46]. This system has one degree of freedom.  X 2 þ 2X  Y ¼ 0

ð13Þ

The tangent stiffness matrix is the first-order derivative of function with respect to X. This matrix is shown in Eq. (14). Moreover, the following substitutions must be introduced. These substitutions are used to in Eqs. (2) and (5). According to Eq. (14), the displacement vector and applied load are the variable and function of Eq. (13), i.e. X and Y, respectively. The variable X range is chosen from 0 to 3. The load factor increment is 0.001 for finding the summit point. The obtained curve of author's method is illustrated in Fig. 27. The results are compatible with the exact solution.   ½K t  ¼ 2  2X; λP ¼ Y; fX g ¼ X ð14Þ

Fig. 25. The dome shell.

300 250 200 Load(kN)

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 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

150 100 50 0 -50

0

0.5

1

1.5

2

2.5

3

3.5

Vertical Deflection(mm)

Fig. 26. The equilibrium path of the dome shell.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎ 1.5 1 0.5 0 -0.5 Y

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 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

15

-1 -1.5 -2 -2.5 -3 -3.5

0

0.5

1

1.5 X

2

2.5

3

Fig. 27. The second-order equation curve.

5. Conclusions In this paper, the DR algorithm is employed to compute the buckling limit load of structures. This load is estimated by the equilibrium path of structures. By setting the work increment of external loads to zero, a new method is proposed for calculating the load factor increment. The suggested technique can find the structural buckling load without complex matrix operations. Prior to reaching the load or displacement limits points; the maximum load is considered as the buckling limit load of the structure. In other words, at first, the equilibrium path of structures is plotted. Based on this path, the corresponding load to the first limit point is known as buckling limit load. Several trusses, frames and shells are analyzed in order to evaluate the capability of the proposed strategy. In this study, both finite element formulation and large deformations are considered. The suggested scheme automatically computes the buckling limit load and only requires the fictitious parameters of the DR iterations. It is noted that the eigenvalue problem is usually solved for finding the buckling load in the other approaches. The conventional process is very time-consuming for large structures. In contrast, the authors' method does not solve the corresponding eigenvalue problem.

References [1] S. Novoselac, T. Ergić, P. Baličević, Linear and non-linear buckling and post buckling analysis of a bar with the influence of imperfections, Tehnički Vjesnik 19 (2012) 695–701. [2] M. Rezaiee-Pajand, H.R. Vejdani Noghreiyan, A.R. Naghavi, Four new methods for finding structural critical points, Mech. Based Design Struct. Mach. 41 (2013) 399–420. [3] J. Alamatian, Displacement-based methods for calculating the buckling load and tracing the post-buckling regions with dynamic relaxation method, Comput. Struct. 114–115 (2013) 84–97. [4] Y.-L. Pi, M.A. Bradford, Multiple unstable equilibrium branches and non-linear dynamic buckling of shallow arches, Int. J. Non-Linear Mech. 60 (2014) 33–45. [5] C.-N. Chen, A finite element study on bifurcation and limit point buckling of elastic-plastic arches, Comput. Struct. 60 (1996) 189–196. [6] R. Szilard, Critical load and post-buckling analysis by fem using energy balancing technique, Comput. Struct. 20 (1985) 277–286. [7] Z. Ju-fen, C. Xiao-ting, An improved arc-length method and application in the post-buckling analysis for composite structures, Appl. Math. Mech. 23 (2002) 1081–1088. [8] Y.Y. Lee, X. Zhao, J.N. Reddy, Postbuckling analysis of functionally graded plates subject to compressive and thermal loads, Comput. Methods Appl. Mech. Eng. 199 (2010) 1645–1653. [9] K. Asemi, M. Shariyat, Three-dimensional biaxial post-buckling analysis of heterogeneous auxetic rectangular plates on elastic foundations by new criteria, Comput. Methods Appl. Mech. Eng. (2016), http://dx.doi.org/10.1016/j. cma.2015.12.026.

[10] M. Shariyat, K. Asemi, Uniaxial and biaxial post-buckling behaviors of longitudinally graded rectangular plates on elastic foundations according to the 3D theory of elasticity, Composite Struct. (2016), http://dx.doi.org/10.1016/j. compstruct.2016.01.065. [11] A.S. Day, An introduction to dynamic relaxation, Engineer 219 (1965) 218–221. [12] P. Underwood, Dynamic relaxation(in structural transient analysis), Computational Methods for Transient Analysis(A 84  29160 12 64), North-Holland, Amsterdam (1983), p. 245–265. [13] S.P. Frankel, Convergence rates of iterative treatments of partial differential equations, Math. Tables Other Aids Comput. 4 (1950) 65–75. [14] K.R. Rushton, Large deflecion of variable-thickness plates, Int. J. Mech. Sci. 10 (1968) 723–735. [15] J.S. Brew, D.M. Brotton, Non-linear structural analysis by dynamic relaxation, Int. J. Numer. Methods Eng. 3 (1971) 463–483. [16] J.W. Bunce, A note on the estimation of critical damping in dynamic relaxation, Int. J. Numer. Methods Eng. 4 (1972) 301–303. [17] A.C. Cassell, R.E. Hobbs, Numerical stability of dynamic relaxation analysis of non-linear structures, Int. J. Numer. Methods Eng. 10 (1976) 1407–1410. [18] M. Papadrakakis, A method for the automatic evaluation of the dynamic relaxation parameters, Comput. Methods Appl. Mech. Eng. 25 (1981) 35–48. [19] S. Qiang, An adaptive dynamic relaxation method for non-linear problems, Comput. Struct. 30 (1988) 855–859. [20] L.G. Zhang, T.X. Yu, Modified adaptive dynamic relaxation method and its application to elastic-plastic bending and wrinkling of circular plates, Comput. Struct. 33 (1989) 609–614. [21] L.C. Zhang, M. Kadkhodayan, Y.W. Mai, Development of the maDR method, Comput. Struct. 52 (1994) 1–8. [22] M. Kadkhodayan, J. Alamatian, G.J. Turvey, A new fictitious time for the dynamic relaxation (DXDR) method, Int. J. Numer. Methods Eng. 74 (2008) 996–1018. [23] M. Rezaiee-Pajand, J. Alamatian, The better time step for geometric non-linear analysis by using of dynamic relaxation method, Modarres Tech. Eng. 19 (2005) 61–74. [24] B.H.V. Topping, P. Ivanyi, Computer Aided Design of Cable Membrane Structures, Saxe-Coburg Publications, United Kingdom, 2008. [25] M. Rezaiee-Pajand, J. Alamatian, The dynamic relaxation method using new formulation for fictitious mass and damping, Struct. Eng. Mech. 34 (2010) 109–133. [26] M. Rezaiee-Pajand, S.R. Sarafrazi, Non-linear structural analysis using dynamic relaxation method with improved convergence rate, Int. J. Comput. Methods 7 (2010) 627–654. [27] M. Rezaiee-Pajand, M. Kadkhodayan, J. Alamatian, L.C. Zhang, A new method of fictitious viscous damping determination for the dynamic relaxation method, Comput. Struct. 89 (2011) 783–794. [28] M. Rezaiee-Pajand, S.R. Sarafrazi, Non-linear dynamic structural analysis using dynamic relaxation with zero damping, Comput. Struct. 89 (2011) 1274–1285. [29] J. Alamatian, A new formulation for fictitious mass of the Dynamic Relaxation method with kinetic damping, Comput. Struct. 90–91 (2012) 42–54. [30] M. Rezaiee-Pajand, M. Kadkhodayan, J. Alamatian, Timestep selection for dynamic relaxation method, Mech. Based Design Struct. Mach. 40 (2012) 42–72. [31] M. Rezaiee-Pajand, H. Rezaee, Fictitious time step for the kinetic dynamic relaxation method, Mech. Adv. Mater. Struct. 21 (2012) 631–644. [32] M. Rezaiee-Pajand, J. Alamatian, Automatic DR structural analysis of snapthrough and snap-back using optimized load increments, J. Struct. Eng. 137 (2011) 109–116. [33] K.S. Lee, S.E. Han, T. Park, A simple explicit arc-length method using the dynamic relaxation method with kinetic damping, Comput. Struct. 89 (2011) 216–233.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 Q7112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

16

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 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

M. Rezaiee-Pajand, H. Estiri / International Journal of Non-Linear Mechanics ∎ (∎∎∎∎) ∎∎∎–∎∎∎

[34] K.-S. Lee, S.-E. Han, J.-W. Hong, Post-buckling analysis of space frames using concept of hybrid arc-length methods, Int. J. Non-Linear Mech. 58 (2014) 76–88. [35] Y.-B. Yang, M.-S. Shieh, Solution method for non-linear problems with multiple critical points, AIAA J. 28 (1990) 2110–2116. [36] M.A.M. Torkamani, J.-H. Shieh, Higher-order stiffness matrices in non-linear finite element analysis of plane truss structures, Eng. Struct. 33 (2011) 3516–3526. [37] G. Ramesh, C.S. Krishnamoorthy, Inelastic post-buckling analysis of truss structures by dynamic relaxation method, Int. J. Numer. Methods Eng. 37 (1994) 3633–3657. [38] G. Ramesh, C.S. Krishnamoorthy, Post-buckling analysis of structures by dynamic relaxation, Int. J. Numer. Methods Eng. 36 (1993) 1339–1364. [39] F. Williams, An approach to the non-linear behaviour of the members of a rigid jointed plane framework with finite deflections, Q. J. Mech. Appl. Math. 17 (1964) 451–469. [40] R.D. Wood, O. Zienkiewicz, Geometrically non-linear finite element analysis of beams, frames, arches and axisymmetric shells, Comput. Struct. 7 (1977) 725–735.

[41] R. Levy, W.R. Spillers, Analysis of Geometrically Non-linear Structures, Chapman & Hall, United States, 1995. [42] D. Boutagouga, A. Gouasmia, K. Djeghaba, Geometrically non-linear analysis of thin shell by a quadrilateral finite element with in-plane rotational degrees of freedom, Eur. J. Comput. Mech. 19 (2010) 707–724. [43] H.-M. Jeon, Y. Lee, P.-S. Lee, K.-J. Bathe, The MITC3 þ shell element in geometric non-linear analysis, Comput. Struct. 146 (2015) 91–104. [44] K.Y. Sze, X.H. Liu, S.H. Lo, Popular benchmark problems for geometric nonlinear analysis of shells, Finite Elements Anal. Design 40 (2004) 1551–1569. [45] M. Rezaiee-Pajand, M. Ghalishooyan, M. salehi, Comprehensive evaluation of structural geometrical non-linear solution techniques Part II: Comparing efficiencies of the methods, Struct. Eng. Mech. 48 (2013) 879–914. [46] E. Carrera, A study on arc-length-type methods and their operation failures illustrated by a simple model, Comput. Struct. 50 (1994) 217–229.

Please cite this article as: M. Rezaiee-Pajand, H. Estiri, Computing the structural buckling limit load by using dynamic relaxation method, International Journal of Non-Linear Mechanics (2016), http://dx.doi.org/10.1016/j.ijnonlinmec.2016.01.022i

67 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 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132