Geometric and material nonlinear analyses of elastically restrained arches

Engineering Structures 29 (2007) 283–295 www.elsevier.com/locate/engstruct

Geometric and material nonlinear analyses of elastically restrained arches Y.-L. Pi ∗ , M.A. Bradford, F. Tin-Loi, R.I. Gilbert School of Civil and Environmental Engineering, The University of New South Wales, Sydney, Australia Received 9 July 2005; received in revised form 15 November 2005; accepted 13 January 2006 Available online 28 November 2006

Abstract An arch is often connected with other structural members and may be supported by elastic foundations that provide elastic-like restraints to the arch. These elastic restraints participate in the structural response of the arch, and may significantly influence its structural behaviour. A curved-beam element for the structural analysis of arches should therefore consider the effects of these elastic restraints. It is known that the elasto-plastic properties of materials also affect the behaviour of a structure significantly, and so the nonlinearities of the materials need to be included in the curved element for elasto-plastic analysis of an arch. An elastic curved-beam element for the in-plane nonlinear analysis of arches is extended here to an elasto-plastic curved-beam element by formulating the elastic restraints and the material nonlinearities into the element. Geometric nonlinearities in the element are based on an accurate rotation matrix of the two dimensional special orthogonal group that satisfies the desired orthogonality and unimodular conditions. Comparisons of results from the curved element with test results and with analytical solutions demonstrate that the curved-beam element can provide accurate predictions for elastic behaviour, and acceptablly accurate predictions for the elasto-plastic behaviour of elastically restrained arches. c 2007 Published by Elsevier Ltd  Keywords: Arches; Curved-beam; Elastic restraints; Elasto-plastic; Finite element; Geometric and material nonlinearity

1. Introduction Geometric nonlinearity is important for the nonlinear inplane analysis of arches. In the nonlinear range, particularly after buckling, the deformations of the arch increase rapidly and become very large, so that predicting the large deformation nonlinear behaviour correctly requires consideration of the effects of large deformations on the deformed curvature and on the axial deformations as pointed out by Pi and Trahair [1]. However, in conventional formulations for curvedbeam elements [2–6], the nonlinear strains under in-plane loading consist only of nonlinear membrane strains and linear bending strains. The higher-order bending strain components have been ignored in the conventional finite-element (FE) formulations of curved-beam elements [2–6]. Pi and Trahair [1] developed a curved-beam element for the in-plane nonlinear elastic analysis of arches by considering the effects of the higher-order curvature components due to bending and the axial deformations of the bending deformation. In addition ∗ Corresponding author.

E-mail address: [email protected] (Y.-L. Pi). c 2007 Published by Elsevier Ltd 0141-0296/$ - see front matter  doi:10.1016/j.engstruct.2006.01.016

to geometric nonlinearity, it is known that the material nonlinearity affects the elasto-plastic behaviour of structures significantly. In order to extend the curved-beam element for the nonlinear elasto-plastic analysis of arches, the material nonlinearities need to be included in the curved-beam element. Axial compression is the major primary action in an arch. In order to produce this axial compression, the supports of the arches are usually fixed or pin-ended. However, in practice, an arch may be supported by elastic foundations or by other structural members that provide elastic-like restraints to the arch. For example, the ends of the arch shown in Fig. 1 are elastically restrained. In many cases, by knowing the structural configuration that connects the arch, the stiffnesses of the elastic restraints can be estimated accurately. These elastic restraints participate in the structural response of the arch and may influence significantly its structural behaviour. The curved-beam element, therefore, also needs to consider the effects of these elastic restraints, but these elastic restraints are usually not included in formulations of curved elements. For example, the curved elements in commercial FE packages such as ABAQUS [7] and ANSYS [8] do not include these elastic

284

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

Fig. 1. Elastically restrained arches.

restraints. Because of this, when these FE packages are used, the elastic restraints need to be modeled using seperate spring elements. The purpose of this paper, then, is to extend quite significantly the elastic curved-beam element of Pi and Trahair [1] by formulating material nonlinearity and elastic restraints into the element for the geometric and material nonlinear large deformation analysis of elastically restrained arches, and to present some of its applications to engineering structures. 2. Finite element model 2.1. Nonlinear strains The Euler–Bernoulli hypothesis, that plane sections that are normal to a centroidal axis remain plane and normal to the deformed centroidal axis after deformation, is used here in the formulation of the curved-beam element. Three sets of axes are used to describe the geometry of the arch, as shown in Fig. 2. The first set is the fixed axes OY Z with basis vectors (PY , P Z ). The second set is a material axis system oys with basis vectors (p y , ps ). After the deformation, the origin o displaces radially v and axially w to o1 , and the oy and os axes move to o1 y1 and o1 s1 , respectively. The third set of axes o1 y1 s1 is also a material axis system with basis vectors (q y , qs ). The axes o1 y1 s1 and basis vectors q y , qs move with the arch during the deformation. The rotations from vectors p y , ps to vectors q y , qs can be described using an orthogonal rotation matrix R as [1]      qy py (v˜  )/(1 + ) (1 + w˜  )/(1 + ) = , (1) −(v˜  )/(1 + ) (1 + w˜  )/(1 + ) qs ps  where (1 + ) = (1 + w˜  )2 + v˜ 2 , v˜  = v  − wκ0 , w˜  = w + vκ0 , () ≡ d()/ds, and κ0 is the initial curvature of the centroidal axis of the arch and is given by κ0 = −1/R for a circular arch with an upward rise; and R is the radius of the initial curvature. The matrix R is a skew-symmetric and satisfies the orthogonal and unimodular conditions of the two dimensional

Fig. 2. Axes and position vectors.

special orthogonal group SO(2) [9] that RRT = RT R = I and det R = 1, where I is the unit matrix. Therefore, the invariant requirement needed for a rigid body rotation is satisfied and so the strains remain invariant during the rigid body rotations. The strains at an arbitrary point P at the cross-section (Figs. 1 and 2) can be obtained by position vector analysis, in conjunction with the rotation relation given by Eq. (1) as [1]  yy =  ys = sy = 0,

(2)

and the longitudinal normal strain ss is obtained as   1 2 1 2 v˜ (1 + w˜  ) − v˜  w˜   ss = w˜ + v˜ + w˜ − y 2 2 [(1 + w˜  )2 + v˜ 2 ]1/2   2 2 1/2 + κ0 [(1 + w˜ ) + v˜ ] − κ0 ,

(3)

which can be reduced to that of [10] and [11] for straight members by setting κ0 = 0. The strains ss include the secondorder membrane strains, as well as the second- and higherorder bending strains. No approximations are made to either the rotation matrix or the longitudinal normal strains. If a small rotation matrix is used [10,12], the higher-order bending strain will be lost.

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

285

Eq. (4) in the FE analysis for the nonlinear stress–strain behaviour of the concrete, the following assumptions and considerations can be used [17]. When the value of the mean cylinder compressive stress σc is not available, the cube strength σcu may be available and then the cylinder compressive stress σc is assumed to be given by [17,18] σc = 0.85σcu.

(5)

When the Young’s modulus is not given, its value is assumed to be [17,18] √ (6) E c = 0.043ρ 1.5 σc , where σc and E c are expressed in MPa and ρ is the density of the concrete and its unit needs to be kg/m3. When the density ρ is not available, ρ = 2400 kg/m3 is assumed. When the value of the strain c corresponding to the maximum stress of the concrete is not given, its value is assumed to be [17,18]

Fig. 3. Stress–strain curve for concrete.

2.2. Material nonlinearity

c = 0.002. Material nonlinearity affects the elasto-plastic behaviour of a structure significantly. Therefore, to perform a nonlinear elastoplastic analysis, material nonlinearies need to be included in the curved element. When material nonlinearity is considered, the constitutive equations for materials are usually used in association with their uniaxial stress–strain, respectively. Hence, stress–strain curves for some commonly used materials such as steel and concrete are formulated for use in the present curved-beam element. The experimental short term uniaxial stress–strain curve for concrete has essentially no linear range and the slope of the curve is continuous up to “failure”. The response of the concrete is nonlinear and, after the ultimate stress is reached, the material softens until it can no longer carry any stress. A large number of empirical stress–strain relationships represented by prescriptive equations have been proposed by many researchers over the years [17]. The following stress–strain equation proposed by Saenz [16] is used in the curved element for the nonlinear compressive behaviour of concrete: σ =

Ec  , 1 + (E c /E scnt − 2)/c + (/c )2

(4)

where σ is the stress,  is the strain, E c is the Young’s modulus of the concrete, E scnt is the secant modulus corresponding to the maximum stress σc and given by E scnt = σc /c as shown in Fig. 3, the maximum stress σc takes the value of the concrete mean cylinder compressive strength, and c is the strain corresponding to the maximum stress σc . The stress–strain curve given by this equation can represent the ascending and descending portions of the nonlinear relationship between the compressive stresses and strains for concrete as shown in Fig. 3, and is consistent with other similar empirical curves [17,18]. It can be seen from Eq. (4) that, to define a stress–strain curve, the Young’s modulus of the concrete E c , its compressive strength σc and the corresponding strain c are essential. However, a number of reported experimental studies do not necessarily give these values. For these cases, in order to use

(7)

When a uniaxial specimen is loaded in tension, its response is elastic until cracks form so quickly that it is very difficult to observe the actual behaviour. For the purpose of modelling, it is assumed that the concrete that is subjected to tension loses strength through a softening mechanism and that the open cracks can be respresented by a loss of elastic stiffness. It is also assumed that the cracks can close completely when the stress across them becomes compressive. The phenomenon of tension stiffening [19] is considered herein for the tensile behaviour of concrete, because tensile stresses are generated in the concrete beyond a crack due to the restraining action by the steel component and the transfer of stresses from the reinforcement and the adjacent uncracked concrete. The bilinear stress–strain curve in Fig. 3 is used to represent the tensile behaviour of the concrete where the Young’s modulus is the same as that for the compressive behaviour, σt is the tensile strength of the concrete, and t u is the maximum tensile strain. The (typical) value of the tensile strength σt can be obtained from splitting tests. In the absence of tested values, the value of σt is assumed to be given by [19] σt = 0.1σc .

(8)

The strain t corresponding to the tensile strength is t = σt /E c . The ratio t u /t of the maximum tensile strain to the strain t is chosen to be inversely proportional to the length of the finite element to avoid mesh dependence. For hot-rolled steel, a tri-linear elastic–plastic-strain hardening stress–strain curve shown in Fig. 4 can be used [13], where E s is the Young’s modulus of the steel, σ y0 is its yield stress, and  y is the strain at which yielding occurs. After the onset of yielding, the steel is assumed to be fully plastic until strain hardening starts at the strain h which is assumed to be n y times the yield strain  y . E tan is the tangent modulus during strain hardening, while the maximum strain u is assumed to be n s times the yield strain  y . The values of E,  y , E s , n y , n s need to be input from experimental results or using assumed values. The values n y = 11 and n s = 31 can usually be used [14,15].

286

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

Fig. 6. Division of cross-sections. Fig. 4. Tri-linear stress–strain curve for steel.

Fig. 5. Ramberg–Osgood rounded stress–strain curve for structural steel.

For steel reinforcement, cold-formed steel section, and deck sheeting, the stress and strain curves display no typical yield characteristics. In these cases, the Ramberg–Osgood stress–strain curves (Fig. 5) are used, which can be expressed as [13]  n σ σ p = + , (9) Es 100 σ p where E s is the Young’s modulus, and σ p is a reference stress and usually takes the value of the 0.2% proof stress, i.e., σ p = σ0.2 . The elastic modulus E s and proof stress σ0.2 can be obtained from experiments, and the parameters p and n are chosen to match the experimental data. The longitudinal normal stress σss is a function of the longitudinal normal strain, and the relationship of the strain and stress can be expressed in an incremental form as δσss = E t δss ,

(10)

where E t is the tangent modulus of the relevant material. The curved element also includes piece-wise linear stress–strain curve input facilities. When the experimental

stress–strain curves are different to the above stress–strain curves, they can be input in a piece-wise linear fashion. Because of the path-dependent nature of plasticity, the crosssections are divided into elemental areas, or into grids of “monitoring points” or “fibres”, as shown in Fig. 6. The current stress, current yield stress, and the current level of equivalent plastic strain must be stored at each monitoring point and updated incrementally throughout the analysis. In this way, the spread of yielding throughout the cross-section can be captured, and the stress resultants and effective cross-sectional properties can be determined by numerical integration over the cross-section. The monitoring points can be obtained by firstly dividing a cross-section into uniformly straight or curved strips. The use of curved strips is for circular solid and circular hollow sections, and for the rounded corners of rectangular and square hollow sections. The monitoring point scheme in each strip depends on the position and orientation of the strip relative to the axis of bending of the cross-section, and the stress gradient through the section. The number and orientation of strips and the number and layout of the monitoring points can be determined according to the problem at hand. 2.3. Nonlinear equilibrium The nonlinear equilibrium equations for a curved element can be derived from the principle of virtual work, which requires that S dU = δss σss dV − δuT (q − p)ds 0 V

δuTk (Qk − Pk ) = 0 (11) − k=1,2

for all admissible sets of infinitesimal virtual displacements δv and δw, where δ() denotes the Lagrange operator of simultaneous variations, q is the vector of distributed loads, Qk is the vector of concentrated loads acting at ends of the element, p is the vector of elastic reactions due to distributed restraints, Pk is the vector of elastic reactions of the restraints at the ends of the element, and the vectors u and uk conjugate the external loads q and Qk and the reactions p and Pk due to

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

287

The variation of the strain δss can be written as δss = cBδϕ = cBNδr,

where c = {1, y}. Because the strains include the higher-order deformed curvature terms caused by the effects of axial deformations, the order of the bending strains is consistent with that of the membrane strains, so that the same low-order cubic polynomials can be used as the shape functions for both the v and w displacements. As a result, the membrane locking problem is avoided, and neither selective reduction integrations nor higher-order interpolation polynomials are needed. From Eq. (14),

Fig. 7. Loading on an element.

δp = sAδϕ = sANδr and



δPk = SAδϕ = SANδrk .

Fig. 8. Restraints on an element.

k=1,2

elastic restraints, respectively, and can be expressed as 

u = {v, v , w} , T

q = {q y , qs , m} , T



uk = {v, v Qk =



, w}Tk

(k = 1, 2);

{Q y , Q s M}Tk T

(k = 1, 2);

(12) (13)

Pk = {k V v, kΦ v  , k W w}Tk = STk uk

(k = 1, 2);

k=1,2

(20)

k=1,2

Substituting Eqs. (16) and (19) into Eq. (11) leads to   T δr [NT BT R − NT AT (q − p)]ds 0

p = {kv v, kφ v , kw w} = s u, T

(19)

−

(14)



NT AT (Qk − Pk ) = 0 ∀δr

(21)

k=1,2

where q y and Q y are the distributed and concentrated forces parallel to the oy axis, qs and Q s are the distributed and concentrated forces tangential to the os axis, and m and M are the distributed and concentrated moments, as shown in Fig. 7; s and Sk are vectors of the stiffness of the distributed and concentrated restraints, respectively, and are given by s = {kv , kφ , kw }T

and Sk = {k V , kΦ , k W }Tk

(k = 1, 2) (15)

in which kv , kφ , kw and k V , kΦ , k W are the elastic stiffnesses of the radial, rotational and axial distributed and concentrated restraints as shown in Fig. 8, respectively. It is noted that the virtual work due to reactions p and Pk of distributed and concentrated elastic restraints are generally not included in previous formulations [1–8]. The vectors δu and δuk of the virtual displacements conjugate to the external loads q and Qk , and the reactions Pk due to the elastic restraints can be expressed as δu = δuk = Aδϕ = ANδr,

(16)

where ϕ is the vector of general displacements of the centroids and is given by ϕ = {v, v  , v  , w, w , w }T = Nr,

(17)

where N is the shape function matrix whose elements are functions of s, and r are the nodal displacements r = {v1 , v1 , w1 , w1 , v2 , v2 , w2 , w2 }T ,

(18)

where v1 , v2 , w1 , w2 are the radial and axial displacements, v1 and v2 are the rotations, and w1 and w2 are the extensions of both ends of the element, which represent usual “engineering” freedoms.

where R is the vector of stress resultants given by R= cT σ d A.

(22)

A

Eq. (21) has to hold for all admissible sets of virtual nodal displacements δr, and so the nonlinear equilibrium equations can be obtained as  NT BT Rds − f = 0, (23) 0

 where 0 NT BT Rds is the vector of the element internal forces conjugate to the vector of virtual nodal displacements δr, while 

NT AT (q − p)ds + NT AT (Qk − Pk ) (24) f= 0

k=1,2

is the conjugate vector of the external loads acting on the

 element, which includes the external forces 0 (−NT AT pds −  T T k=1,2 N A Pk produced by reactions p and Pk of the elastic restraints. 2.4. Consistent linearization The consistent linearization of the nonlinear equilibrium formulation plays a key role in numerical implementations employing an incremental-iterative solution procedure. A complete account of linearization procedures in the general context of an infinite dimensional manifold can be found in [20]. Linearization of the nonlinear equilibrium can be obtained by taking the variations of the nonlinear equilibrium equations (23) as

288





Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

  NT δBT R + BT δR ds

0





+ +



δA (p − q) + A (δp − δq) ds   NT δAT (Pk − Qk ) + AT (δPk − δQk ) = 0.

N 0





3. Implementation

T

T

T

(25)

k=1,2

j

kTi ri = λ1i fref ,

From Eqs. (19) and (22) δR = cT δσ d A = cT E t cBδθ d A = DBδθ = DBNδr, A

A

(26) where the matrix D is given by D= cT E t cd A.

(27)

A

The terms δBT R can be written as the identities δBT R = Mσ δθ = Mσ Nδr,

(28)

where Mσ is a symmetric matrix that accounts for the nonlinear effects of the stress resultants on the tangent stiffness matrix. Since the entries of the matrix A are constants, δA = 0. Substituting Eqs. (26) and (28) into Eq. (25) leads to the first term of Eq. (25) as   

NT δBT R + BT δR + AT δp ds + AT δPk 0



k=1,2

  = NT BT DB + Mσ + AT sA Nδrds 0

+ NT AT SAN. 

(29)

k=1,2

Since δA = 0, the consistent linearization of the second and third terms of Eq. (25) is then obtained as 

δp = NT AT δqds + NT AT δQk . (30) 0

The FE model can be implemented easily in an incrementaliterative procedure using a modified Newton–Raphson method [21]. The formulation of the incremental-iterative procedure can be established by combining Eqs. (23) and (31) as

k=1,2

Substituting Eqs. (29) and (30) into Eq. (25) leads to the tangent stiffness relationship kT δr = δp, where the tangent stiffness matrix kT is given by  kT = NT (BT DB + Mσ + AT sA)Nds 0

+ NT AT SAN.

(31)

(33)

where i and j denote the load step and the iteration within the load step, respectively, and the reference load fref can be specified as input data. Each load step i consists of the application of an increment of external load and subsequent iteration j to restore equilibrium. In the modified Newton–Raphson method, the tangent stiffness matrix kTi is formed at the commencement of the load step i and is then held constant throughout the equilibrium iteration. The load factor λ11 for the first iteration at the first load step is specified as input data. For the first iteration at the second and subsequent load steps, the load factor λ1i is determined by using the strategy  Nreal 1 1 , (34) λi = λi−1 Ndesired where Ndesired is the input desired number of iterations and Nreal is the actual number of iterations required for convergence in the load step (i − 1). At the end of the first iteration of each step, the modified Euclidean norm    N  1

rk 2   (35) ξ= N k=1  rk,ref  is used to check convergence, where N is the total number of degrees of freedom, rk is the increment of the displacement component k, and rk,ref is the largest displacement of the corresponding type. If ξ ≤ η, convergence has been deemed to be attained where η is a specified tolerance, typically within the range 10−1 –10−5 . If ξ > η, the iteration has not converged and needs to be continued. The convergence checking by Eq. (35) is repeated at the end of every iteration. When the iteration ( j − 1) does j not converge, the vector of unbalanced forces ψ i for the j th iteration is computed as     j NT BT R − AT (q − p) ds ψi = 0

(32)

−



i N A (Qk − Pk )

k=1,2

k=1,2

The tangent stiffness matrix kT is symmetric and includes

 the contributions of the elastic restraints 0 NT AT sANds +  T T k=1,2 N A SAN to the tangent stiffness matrix given by Eq. (32). It is worth pointing out that effects of elastic restraints on both the nonlinear equilibrium given by Eq. (23) and the tangent stiffness matrix need to be considered simultaneously.

T

,

T

(36)

j

where the matrices B, A and the stress resultants R are calculated based on the displacements ri−1 at the end of the last converged load step (i − 1). j The displacement increments rri due to the unbalanced j forces ψ i are obtained from j

j

kT i rri = ψ i .

(37)

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

289

During the iterative process, the change in the incremental j load factor λi is determined by the arc-length restrained method proposed by Crisfield [22], and the automatic incrementation of the arc length is adopted [22]. The sign of the load increment follows the sign of the determinant of the tangent stiffness matrix. j Once λi has been obtained, the increment of the displacement vector in the current iteration can be expressed as j

j

ri =

λi

λ1i

j

r1 i + rri .

(38)

Two methods can be used for updating the strains. The first method is to use the incremental strains when the iterations have converged, while the second method is to use the iterative strains in the each iteration. Since the present FE model includes the reactions due to elastic restraints in the external forces given in Eq. (23) and the unbalanced forces have to be calculated in every iteration, the second method is used here for more efficient unbalanced force corrections.

Fig. 9. Buckling load of rotationally restrained arches under uniform radial loads.

4. Applications

The analytical solution for antisymmetric bifurcation buckling is given by [24]

4.1. Nonlinear buckling of arches with rotational end restraints

A1 q¯ 2 + B1 q¯ + C1 = 0,

The FE model has been used to calculate the elastic buckling load of three groups of arches (0.5S/r x = 50, 75, 100, with S being the length of the arch and r x the radius of gyration of its cross-section) with end rotational restraints under a radial load q uniformly distributed around the arch axis. The included angles 2Θ of the three arch groups are assumed to increase from 2Θ = 0◦ to 2Θ = 180◦. In the FE analyses, the arches are assumed to have a uniform 260UB37 I-section [23], with the overall depth D = 0.256 m, the flange width B = 0.146 m, the flange thickness t f = 0.0109 m, web thickness tw =√0.0064 m, the radius of gyration of the cross-section r x = Ix /A = 0.108 m, and A and Ix are the area and second moment of area of the cross-section about its major principal axis. The lengths S = 10.8, 16.2, 21.6 m (Fig. 1) for each group of arches, respectively. Young’s modulus is assumed to be E = 200 000 MPa. For the investigation, the dimensionless ratio α of the bending stiffness E Ix of the arch to the stiffness k of the rotational restraint (defined by α = E Ix /k S) is assumed to be equal to 1. Eight elements were used in each FE model. Variations of the dimensionless buckling load q R/N Eηb with the included angle 2Θ obtained from FE results are compared in Fig. 9 with analytical solutions [24] for the in-plane buckling load of shallow arches, where R is the radius of the initial curvature of an arch, N Eηb is the second mode buckling load of a corresponding pin-ended column of the same length S and with the same end rotational restraints about its major principal axis, and N Eηb is given by N Eηb =

(ηb π)2 E Ix , (S/2)2

(39)

in which ηb is the solution of the transcendental equation ηb π tan(ηb π) = . (40) 1 + 2α(ηb π)2

(41)

where the dimensionless radial load q¯ is defined as q R − N¯ , (42) N¯ in which N¯ is the actual axial compressive force in the arch and the coefficients A1 , B1 and C1 are given by q¯ =

A1 = 5 + 40α + 60α 2 + 8α 2 ηb2 π 2 , B1 = C1 =

(43)

4(1 + 8α + 12α + 4α 2 ηb2 π 2 ), 12(1 + 2α + 4α 2 ηb 2 π 2 )(ηb2 π 2 ) , λ2 2

(44) (45)

in which λ is the modified slenderness of an arch defined by λ=

S2 RΘ 2 . = rx 4r x R

(46)

The solution from Eq. (41) is suitable for shallow arches with the modified slenderness λ ≥ λsb where λsb is given by λsb =

 (ηb π) 3(5 + 40α + 60α 2 + 8α 2 ηb2 π 2 )(1 + 2α + 4α 2 ηb 2 π 2 ) (1 + 8α + 12α 2 + 4α 2 ηb2 π 2 )

.

(47)

Arches with a modified slenderness λ < λsb buckle in a symmetric snap-through mode. To obtain the theoretical solution for the symmetric snap-through buckling load, one needs to solve the equation for nonlinear equilibrium and the equation for snap-through buckling equilibrium simultaneously [24]. Based on the solutions, the approximation for the symmetric snap-through buckling load of shallow arches with a modified slenderness λ < λsb was proposed [24] as q R = (0.15 + 0.006λ2)N Eηb .

(48)

It can be seen from Fig. 9 that, for shallow arches with included angle 2Θ smaller than 90◦ , the FE results are almost

290

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

Fig. 10. Snap-through buckling and postbuckling of rotationally restrained arches under central concentrated load.

Fig. 11. Bifurcation buckling and postbuckling of rotationally restrained arches under central concentrated load.

the same as the analytical solutions. It is worth pointing out that the FE results are the nonlinear buckling loads, but no theoretical solutions for the nonlinear buckling of deep arches with included angle larger than 90◦ are available in the open literature.

1 vc = 2 R μ

4.2. Nonlinear elastic analysis of arches with end rotational restraints The FE model has also been used to predict the nonlinear elastic behaviour of two arches with end rotational restraints and under a central concentrated load. The modified slenderness of the two arches is λ = 8 (included angle 2Θ = 18.4◦ and length S = 10.8 m) and λ = 20 (2Θ = 23◦ , S = 21.6 m). In the FE analysis, the 250UB37 I-section was used and the dimensionless ratio α of the bending stiffness E Ix of the arch to the stiffness k of the rotational restraint was assumed to equal 0.5. Four elements were used in the FE model and error tolerance η = 0.0001 was used. In the FE analysis, to trigger the antisymmetric bifurcation buckling, very small initial equal and opposite radial loads were applied at the quarter points from either support and then removed. Variations in the dimensionless central radial displacement vc / f with the dimensionless central load Q¯ obtained from the FE results for the nonlinear elastic responses of the two arches are compared with the analytical solutions in Figs. 10 and 11, where f is the initial rise of the arch, and the dimensionless load Q¯ is defined as Q R2 Θ Q(S/2)(R/2) = . Q¯ = E Ix 2E Ix

(49)

To show the effects of the elastic restraints, the FE results and analytical solutions for pin-ended arches (α = ∞) are also shown in Figs. 10 and 11. The analytical nonlinear relationship between the dimensionless load Q¯ and the displacement vc / f can be obtained from the expression for the nonlinear central radial displacement given by [24]

 1−

 1 (μΘ )2 Λ+ cos(μΘ ) 2

  Q¯ (sec(μΘ ) − 1)Λ − μΘ , + Ξ+ μΘ μΘ (1 + 2α)

(50)

where the dimensionless parameters μ, Λ and Ξ are defined by N¯ R 2 , E Ix 1 + 2α , Λ= 2α + tan(μΘ )/(μΘ ) μ2 =

(51) (52)

and Ξ =

tan(μΘ )[2α + tan(μΘ /2)/(μΘ )] . 2α + tan(μΘ )/(μΘ )

(53)

In Eq. (50), the parameter μ can be obtained from the nonlinear relationship between Q¯ and μ, given by A2 Q¯ 2 + B2 Q¯ + C2 = 0,

(54)

where the coefficients A2 , B2 and C2 are given by  1 μΘ − sin(μΘ ) cos(μΘ ) 2 Ξ A2 = 4 4 4μ Θ μΘ 2(1 − cos(μΘ ))2 Ξ +3 μΘ  sin(μΘ )[4 − cos(μΘ )] , − μΘ  1 (sin(μΘ ) cos(μΘ ) − μΘ ) ΛΞ B2 = 4 4 μ Θ 2 cos(μΘ )  (1 − cos(μΘ ))2 + Λ , 2 cos(μΘ )   μΘ 2 1 C2 = + 2 λ 4μ Θ 2   μΘ − sin(μΘ ) cos(μΘ ) 2 2(μΘ )2 Λ × − . μΘ cos2 (μΘ ) 3 −

(55)

(56)

(57)

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

291

Table 1 Dimensionless buckling loads Upper buckling load Q¯ Theory Elements 8

4

2

Lower buckling load Q¯ Theory Elements 8

Values Errors

3.37 –

3.36 0.3%

3.35 0.6%

3.45 2.4%

1.61 –

1.61 0%

Values Errors

5.87 –

5.86 0.17%

5.85 0.34%

– –

−1.415 –

−1.40 0.35%

Arch

λ=8 λ = 20

4

2 1.62 0.6%

−1.39 0.7%

1.82 13% – –

Fig. 12. Efficiency of the curved element.

Fig. 13. Efficiency of the curved element.

It can be seen from Figs. 10 and 11 that the FE results for the nonlinear behaviour of both arches are almost identical to the analytical solutions [24]. It can also be seen that the arch with the slenderness λ = 8 may buckle in a symmetric snap-through and snap-back mode, while the arch with the slenderness λ = 20 may buckle in an antisymmetric bifurcation mode. Comparisons of the results of the rotationally restrained arches with those of pin-ended arches demonstrate that the end rotational restraints have significant effects on the nonlinear behaviours of arches. To show the efficiency of the present curved beam-element, two-element and four-element models were also used to obtain FE results for the arch with the modified slenderness λ = 8 (included angle 2Θ = 18.4◦ and length S = 10.8 m) and their results are compared in Fig. 12 and Table 1. FE results for the antisymmetric buckling and postbuckling behaviour of the arch with λ = 20 obtained by the four-element model are shown in Fig. 13 and Table 1. It can be seen from Table 1 that the errors of four- and eight-element models range from 0% to 0.7%. This indicates that the predictions of the present curved-beam element are quite accurate. Because very small initial equal and opposite radial loads need to be applied at the quarter points from each support to trigger the antisymmetric bifurcation buckling, the two-element model cannot be used for antisymmetric bifurcation buckling analysis. It can be seen that the results of the four-element model are almost indentical with those of the eight element model. The result of the two-element model has some errors for the snap-through buckling loads. The errors of the two-element model increase with an increase

in the deformation as in Fig. 12. Therefore, to obtain accurate predictions, the two-element model should not be used. 4.3. Nonlinear elasto-plastic analyses of concrete-filled steel RHS and reinforced concrete rectangular section arches To investigate further the combined effects of the elastic restraints and material nonlinearities, the in-plane nonlinear elastic and inelastic responses of concrete-filled steel rectangular hollow section (RHS) pin-ended and rotationally restrained circular arches (α = 0.1) that are subjected to a central concentrated load have been investigated using the present FE model. An Australian RHS section (150 × 102.5) [25] was used for the bare RHS section arch and the concrete-filled steel RHS section arch. For comparison, the dimensions of the reinforced concrete section were chosen to be the same as the RHS. Their dimensions are shown in Fig. 14. The length of the arch S = 8 m. The material properties for the steel RHS section are: Young’s modulus E s = 200 000 MPa, yield stress σ y = 250 MPa, yield strain  y = 0.00125, a trilinear stress–strain curve (Fig. 4) was assumed, the strain at which hardening starts h = 0.01375, and the tangent modulus E tan = 6000 MPa. The material properties for the concrete were assumed to be: Young’s modulus E c = 31 950 MPa, compressive strength σc = 40 MPa, strain corresponding to the compressive strength c = 0.002, the compressive crush strain crush = 0.005, tensile strength σt = 4 MPa, and the ratio of the maximum tensile strain to the strain corresponding the tensile strength was assumed to be 11. Two arches with included angles

292

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

Fig. 14. Pin-ended arches with concrete-filled steel RHS, bare steel RHS, or reinforced concrete rectangular section.

Fig. 15. In-plane elastic behaviour of concrete-filled steel RHS, bare steel RHS, and reinforced concrete rectangular section arch.

2Θ = 20◦ and 2Θ = 120◦ were investigated. For comparison, the nonlinear elastic and inelastic behaviour of steel RHS arches with the same length, included angle and cross-section without the concrete infill and of reinforced concrete arches with the same dimensions and the same concrete properties but without the steel encasement were also investigated. The material properties of the reinforcement were assumed to be: Young’s modulus E s = 200 000 MPa, the proof stress σ p = σ0.2 = 450 MPa, and the Ramberg–Osgood rounded stress–strain curve (Fig. 5) given by Eq. (9) was assumed. The parameter p = 0.2, corresponding to the 0.2% proof stress σ0.2 , is used with n = 15. Eight elements were used for each arch and error tolerance η = 0.001 was used. Again, to trigger the bifurcation buckling, very small initial equal and opposite radial loads were applied at quarter points from either support

and then removed. The elastic nonlinear responses of these arches are shown in Fig. 15, while their inelastic counterparts are shown in Fig. 16. It can be seen that the concrete-filled steel RHS arches have the highest elastic and inelastic strengths, while the reinforced concrete arches have the lowest elastic and inelastic buckling strengths. The analytical elastic buckling loads for shallow pin-ended arches without restraints given by [26]    π4 Θ Q = 1.33 + 4.5 1 − 0.65 2 N p , with λs π 2 Np =

π 2 E Ix (S/2)2

are shown in Fig. 15.

(58)

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

293

Fig. 16. In-plane inelastic behaviour of concrete-filled steel RHS, bare steel RHS, and reinforced concrete rectangular section arch.

The approximations for the antisymmetric bifurcation buckling load of a non-shallow pin-ended arch without restraints that was proposed by Bradford et al. [26] are    π 2 Θ Q = 5.83 − 0.85 Θ − Np. (59) 2 π2 Approximations for the antisymmetric bifurcation buckling load for arches with included angle 2Θ = 120◦ given by Eq. (59) are shown in Fig. 15 as well. It can be seen that the FE results for the buckling loads are the same as the analytical and approximate results of Bradford et al. [26]. It can be seen from Figs. 15 and 16 that the rotational restraints affect both the buckling mode and the buckling load. While the pin-ended arches buckled in an antisymmetric bifurcation mode, the rotationally restrained arches buckled in symmetric snap-through mode. The elastic restraints increase the elastic and inelastic buckling load of the deep arches (2Θ = 120◦ ) significantly. Comparisons of Fig. 15 with Fig. 16 demonstrated that, for all arches, the inelastic load-carrying capacity is much lower than the corresponding elastic load-carrying capacity, because the material nonlinear inelasticity plays an important role. It can also be seen that the ranges of displacements vc for inelastic analysis are smaller than those for the elastic analysis, because the maximum strains for the materials are limited. When the maximum strains are reached, the inelastic analysis automatically stops. It is worth pointing out that the compressive crush strain for concrete crush = 0.005 was the dominant maximum strain. Because of this, strain hardening did not occur for steel in the reinforced concrete and concrete-filled steel RHS arches. For bare steel arches, strain hardening did not occur when inelastic

buckling occurred, but strain hardening did occur when the displacements became very large. 4.4. Comparison with experimental results Wang et al. [27] tested shallow parabolic tied reinforced concrete arches. The present FE model was used to analyze the nonlinear inelastic behaviour of a tested arch. The arch was simply supported at both ends and two horizontal steel ties of 125×75×12 angle section were welded to both ends of the arch (Fig. 17). The arch was subjected to a uniformly distributed vertical load along the entire span of the arch. The cross-section of the arch was rectangular, with a width B = 500 mm and a depth t = 50 mm. Two layers of 6 mm SL62@200 steel mesh were placed to provide nominal reinforcement with a top and bottom cover of 10 mm. The span and rise of the arch are L = 4000 mm and f = 200 mm, respectively. Because the stress–strain curve for concrete was not obtained in the test, the concrete stress–strain curve given in Fig. 5 was used in the FE computation. Young’s modulus and the cylinder compressive strength of the concrete from the standard cylinder test were E c = 25 300 MPa and σc = 34.3 MPa, respectively. The strain c corresponding to the strength σc was assumed to be c = 0.0035 in the FE runs. The Ramberg–Osgood stress–strain curve was used for the steel mesh. The proof stress σ p = σ0.2 = 400 MPa and Young’s modulus E s = 205 000 MPa were obtained from tests. The parameters p = 0.2 and n = 15 were assumed. The tied arch can be considered to be restrained elastically at its ends by the horizontal tie. The stiffness can be calculated from the dimensions and Young’s modulus of the angle tie as E A/L = 113 000 kN/m, where Young’s modulus of elasticity for steel was assumed to be E s = 200 000 MPa. Eight elements

294

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295

Fig. 17. Parabolic arch tested by Wang et al. [27].

Fig. 18. Inelastic buckling and postbuckling of a tied parabolic arch: comparison with experimental results.

were used and the error tolerance η = 0.001 was used. The FE results are compared with the experimental results in Fig. 18, and it can be seen that the predicted maximum load carrying capacity q = 36.9 kN/m is slightly higher than the experimental load carrying capacity q = 34.2 KN/m (with an error 7.9%), which is acceptably close to the experimental result. One of the possible error sources may be that the stress–strain curve given in Fig. 5 was somewhat different to the real curve. In addition, the angle section tie was under eccentric tension, and so the effective area of the angle section as a tie may be smaller than its nominal area and so the restraining stiffness is somewhat smaller than the above nominal value. 5. Conclusions This paper has compresively extended an elastic curvedbeam element of Pi and Trahair [1] by formulating material nonlinearity and elastical restraints into the element for in-plane elasto-plastic large deformation analysis of elastically restrained arches. The stress–strain curves of some commonly used materials have been formulated in the element. An input

facility for a piece-wise linear stress–strain curve is also available in the present element. To consider the path-dependence of material plasticity, the cross-sections were divided into elemental areas to capture the current stress, current yield stress, and the current level of equivalent plastic strain throughout the analysis. The FE model was applied to the nonlinear elastic analysis of arches with end rotational restraints, to the nonlinear elastoplastic analyses of concrete-filled steel RHS and reinforced concrete rectangular section arches, and to the nonlinear buckling calculations of arches with rotational end restraints. The FE model was also used to predict the nonlinear behaviour of a tested shallow parabolic tied reinforced concrete arch. The comparisons of the FE results with the results of the analytical solutions for the nonlinear analysis of elastically restrained arches and with the test results demonstrated that the curvedbeam element provides accurate predictions for the nonlinear elastic and acceptably accurate predictions for the elasto-plastic behaviour of elastically restrained arches. The present curvedbeam element can be further developed to include the ability to handle initial geometric imperfections and residual stresses. Acknowledgments This work has been supported by the Australian Research Council through Discovery Projects awarded to the authors and through an Australian Research Council Federation Fellowship awarded to the second author and an Australian Professorial Fellowship awarded to the fouth author. References [1] Pi Y-L, Trahair NS. Non-linear buckling and postbuckling of elastic arches. Engineering Structures 1998;20(7):571–9. [2] Wen RK, Suhendro B. Nonlinear curved-beam element for arch structures. Journal of Structural Engineering, ASCE 1991;117(11):3496–515. [3] Noor AK, Peters JM. Mixed model and reduced/selective integration displacement models for nonlinear analysis of curved beams. International Journal for Numerical Methods in Engineering 1981;17(4):615–31. [4] Stolarski H, Belytschko T. Membrane locking and reduced integration for curved elements. Journal of Applied Mechanics, Transaction ASME 1982;49(1):172–6.

Y.-L. Pi et al. / Engineering Structures 29 (2007) 283–295 [5] Calhoun PR, DaDeppo DA. Nonlinear finite element analysis of clamped arches. Journal of Structural Engineering ASCE 1983;109(3):599–612. [6] Elias ZM, Chen K-L. Nonlinear shallow curved-beam finite element. Journal of Engineering Mechanics, ASCE 1988;114(6):1076–87. [7] ABAQUS. ABAQUS Theory Manual, Version 6.4. Pawtucket (RI, USA): Hibbit, Karlsson and Sorensen Inc.; 2003. [8] ANSYS. ANSYS Multiphysics 8.0. Canonsberg (PA, USA): Ansys Inc.; 2003. [9] Burn RP. Groups. A path to geometry. Cambridge (UK): Cambridge University Press; 2001. [10] Timoshenko SP, Gere JM. Theory of elastic stability. 2nd ed. New York (NY): McGraw-Hill; 1961. [11] Thompson JMT, Hunt GW. A general theory of elastic stability. New York (NY): John Wiley & Sons; 1973. [12] Valsov VZ. Thin-walled elastic beams. 2nd ed. Jerusalem: Israel Program for Scientific Translation; 1961. [13] Lemaitre L, Chaboche J-L. Mechanics of solid materials. Cambridge (UK): Cambridge University Press; 1994. [14] Bild S, Chen G, Trahair NS. Out-of-plane strength of steel beams. Journal of Structural Engineering, ACSE 1992;118(8):1987–2003. [15] Pi Y-L, Trahair NS. Nonlinear inelastic analysis of steel beam-columns. II: Applications. Journal of Structural Engineering, ACSE 1994;120(7): 2062–85. [16] Saenz LP. Discussion of equation for the stress–strain curve of concrete. Journal of The American Concrete Institute 1964;61(9):1227–39. [17] Warner RF, Rangan BV, Hall AS, Faulkes KA. Concrete structures.

295

Melbourne (Australia): Longman; 1999. [18] Oehlers DJ, Bradford MA. Composite steel and concrete structural members. Oxford (UK): Pergamon; 1995. [19] Gilbert RI, Warner RF. Tension stiffening in reinforced concrete slabs. Journal of the Structural Division, ASCE 1978;104(ST12):1885–900. [20] Marsden J, Hughes TJR. Topics in the mathematical foundation of elasticity. In: Knops RJ, editor. Nonlinear analysis and mechanics: HeriotWatt Symposium, vol II. Research note in Mathematics, vol. 27. London (UK): Pitman; 1978. [21] Zienkiewicz OC, Taylor RL. The finite element method. 4th ed. New York (NY, USA): McGraw-Hill; 1989. [22] Crisfield MA. A fast incremental/iterative solution procedure that handles snap-through. Computers and Structures 1981;13(1–3):55–62. [23] BHP. Hot rolled and structural steel products. 2000 ed. Melbourne (Australia): BHP Co. Pty Ltd; 2000. [24] Pi Y-L, Bradford MA, Tin-Loi F. In-plane stability of elastically restrained arches. UNICIV report. Sydney (Australia): School of Civil and Environmental Engineering, The University of New South Wales; 2005. [25] Tubemakers of Australia. Tubeline safe load tables. 2nd ed. Sydney (Australia): Tubemakers of Australia Limited; 1984. [26] Bradford MA, Uy B, Pi Y-L. In-plane stability of arches under a central concentrated load. Journal of Engineering Mechanics, ASCE 2002;128(7):710–9. [27] Wang T, Bradford MA, Gilbert RI. Short-term experimental study of a shallow reinforced concrete parabolic arch. Australian Journal of Structural Engineering 2005;6(1):53–60.