- Email: [email protected]
Slip locking in finite elements for composite beams with deformable shear connection
Finite Elements in Analysis and Design 40 (2004) 1907 – 1930 www.elsevier.com/locate/!nel
Slip locking in !nite elements for composite beams with deformable shear connection Andrea Dall’Astaa;∗ , Alessandro Zonab a
Dipartimento di Progettazione e Costruzione dell’Ambiente, Universita di Camerino, Viale della Rimembranza, 63100 Ascoli Piceno, Italy b Istituto di Scienza e Tecnica delle Costruzioni, Universita Politecnica delle Marche, Via Brecce Bianche, 60131 Ancona, Italy Received 15 September 2002; received in revised form 21 October 2003; accepted 12 January 2004
Abstract In this paper, the authors discuss problems that may occur in the !nite element analysis of composite beams with deformable shear connection, where slip arises between beam components while contact is preserved. E2ects of locking problems, comparable to the membrane and shear locking that may develop in !nite elements for common beams, are illustrated. These locking problems mostly a2ect the description of the curvature and interface slip. A strategy to avoid locking by means of a calibrated choice of the displacement shape functions is then proposed and a comparison between di2erent displacement-based elements, in terms of local error, global error and convergence rate, is presented. The authors also illustrate the e2ects of slip locking on mixed !nite elements and comparisons between displacement-based and mixed elements are shown. ? 2004 Elsevier B.V. All rights reserved. Keywords: Composite beams; Shear connection; Locking; Slip oscillations; Mixed formulation; Beam !nite elements
1. Introduction Composite beams, made by two components joined by shear connectors to form an interacting unit, are frequently used in structural engineering. Common examples are steel-concrete composite decks for bridges construction. The behaviour of composite beams is in=uenced by the type of connection between the two parts. Flexible shear connectors permit the development of partial composite action [1] and, for a realistic analysis, the structural model must account for the interlayer slip between the components. Therefore, a simple analytical formulation that captures the salient features of the global ∗
Corresponding author. Tel.: +39-071-220-4551; fax: +39-071-220-4576. E-mail address: [email protected] (A. Dall’Asta).
0168-874X/$ - see front matter ? 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2004.01.007
1908
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
x G1
z
y1
w1 k y2
h
w2 k
G2
y
vj v'
y vj
v' δ
Fig. 1. Newmark kinematical model.
behaviour (bending behaviour and interface slip) is very useful both in the analysis and design of this structural system. If the study is limited to the in-plane bending behaviour, i.e. symmetric cross-sections are considered and no torsion and out-of-symmetry-plane bending occur, the analysis of this type of structure can be based on the Newmark kinematical model [2]: the Euler–Bernoulli linear beam theory (small displacements, rotations and strains) is used to model the two parts of the composite beam; the e2ects of the deformable shear connection are accounted for by using an interface model with distributed bond, preserving the contact between the components (Fig. 1). The interface slip, computed from the di2erence of small displacements of the two parts at the interface, is accordingly small. Once a kinematical model is de!ned for the problem, numerical solutions can be obtained by means of the !nite element method. The !nite element analysis of composite beams makes it possible to obtain useful information about the structural behaviour but requires an eHcient and reliable element. Quite a number of works on this topic can be found in the literature. Di2erent displacement-based !nite element models have been adopted (see [3] and the references quoted there). Finite elements based on the displacement approach have a simple formulation but their behaviour is not always satisfactory. More recently, models that attempt to overcome the limitations of the displacement-based formulation have been proposed: Salari et al. [4], Salari et Spacone [5], Ayoub et Filippou [6], Ayoub [7], Cas et al. [8] and Dall’Asta et Zona [9]. All of these !nite element models indicate the great interest in solving the convergence problems that a2ect the analysis of composite beams with deformable shear connection. This paper aims at clarifying some of the aspects involved in the !nite element solution. For typical structural problems, the best approximation property of the !nite element method [10] is a signi!cant guarantee of reasonable behaviour even when very crude meshes are employed. Generally, when two or more !elds are coupled in a structural model, the constraint between !elds may cause troubles when applied to !nite dimensional spaces, as those used in !nite element approximations. The accuracy of the solution depends on some characteristic parameters involved in the coupled terms. For limit values, such as zero or in!nity, further relations between unknowns reduce the dimension of the space in which the solution is sought [11,12]. In some cases, the dimension goes to zero and the model completely “locks”. In general, a sti2ener response and spurious strains
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1909
can be observed when the phenomenon occurs. Examples of these problems for beam elements include membrane locking [13,14] that may a2ect curved beam elements by varying the curvature; shear locking [14,15] that may develop in Timoshenko beam elements by varying the shear sti2ness; the eccentricity issue [16–18] that may involve the ordinary Euler–Bernoulli beam by varying the origin of the reference system. In all of these cases, the strain parameters involved are contributed by di2erent displacement !elds. A similar problem may a2ect composite beam elements. In fact a coupling arises between the =exural displacement !eld and the axial displacements !elds of the two components due to the deformable shear connection. As shown in this paper, if the approximations of the axial and =exural displacement !elds are not consistent, then the error of the !nite element solution depends strongly on the connection sti2ness, slip oscillation and poor approximation of the curvature occur for high values of the connection sti2ness. This is a problem of particular interest in the non-linear analysis of composite beams [3] because the non-linear behaviour of the shear connection is characterised by a constitutive law with a sti2ness varying from very a high value when the slip is close to zero, to a very low value when the connection is close to failure. Since the distribution of interface slip varies greatly along the beam, the solution obtained, which depends strongly on the connection sti2ness, is not accurate. In this paper, the authors analyse the e2ectiveness of !nite elements based on the Newmark kinematical model [2] in linear analyses. The inconveniences that arise in the displacement-based element, with shape functions of the lowest degree required by the problem (8DOF element), are illustrated. A strategy for avoiding slip locking by means of a calibrated choice of axial and transverse displacement shape functions is proposed. The authors also illustrate the e2ects of slip locking in mixed elements based on the 8DOF displacement !eld approximation. The proposed comparisons between displacement and mixed elements are intended to clarify whether a mixed formulation is a2ected by slip locking or if it can be used to reduce or eliminate slip locking problems.
2. Kinematical model and main equations 2.1. Kinematical model In the reference con!guration, the composite beam occupies the region V = A × [0; L] generated by translating its symmetrical cross-section A along a rectilinear axis, orthogonal to the cross-section. The cross-section is divided in two parts A = A1 ∪ A2 . We introduce an orthonormal reference frame {O; X; Y; Z}, where Z is orthogonal to A. The unit vectors parallel to the axis X , Y , Z are denoted by i, j, k, respectively. The connection is assumed to lie on a plane that is orthogonal to the Y -axis and permits only discontinuities parallel to the beam axis. It is also assumed that YZ is a symmetry plane for the problem and shear deformation cannot occur; in other words, every component deforms like plane Kirchho2 rod, and consequently the displacement !eld has the following representation: u (y; z) = v(z)j + [w (z) + (y − y)v (z)]k
in A × [0; L] ( = 1; 2);
(1)
where the prime denotes the derivatives with respect to the z co-ordinate and u ( = 1; 2) can be deduced from the three scalar functions v and w describing the displacement of the centroid of A
1910
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
in Y and Z directions (y is the ordinate of the centroid of A , see Fig. 1). The functions furnishing the displacement !eld can be collected in the vector uT = [w1
w2
v]:
(2)
2.2. Strain and stress parameters The kinematical assumptions produce a reduced set of deformations in the beam. More precisely, only the axial strain z (y; z) = w + (y − y)v = (z) + (y − y )(z)
in A × [0; L] ( = 1; 2)
(3)
exists ( = w is the axial strain evaluated at the centroid of A and = −v is the curvature), while the slip at the interface is expressed by the following relation: (z) = w2 (z) − w1 (z) + hv (z)
in [0; L];
(4)
where h=y2 −y1 is the distance of the two components centroids (Fig. 1). The parameters describing the strain measures can be collected in the vector ” which is related to the displacements u by means of the di2erential operator D (compatibility equations): 9 0 0 1 w1 0 2 0 9 w = Du: ”= (5) = 2 2 0 0 −9 v −1 1 h9 The Virtual Work Principle permits identifying the dynamic entities that are duals of the kinematic entities introduced by the kinematical assumptions. Such dynamic entities are resultants of the stresses which make work (active stress). The Virtual Work Principle speci!ed for the problem under examination becomes: L L r · Duˆ d z = g · Huˆ d z ∀uˆ in [0; L] (6) 0
0
(uˆ is the admissible displacement !eld) where the resultant of active stress and external forces are collected in the two vectors rT = [N1
N2
g T = [gz1
gz2
M12
fs ];
(7)
gy
mx ]
(8)
with (see [9] for more details) N ( = 1; 2) the axial forces on A , M12 the summation of the bending moments M on A , fs the interface shear force, gz1 , gz2 , gy , mx the two axial distributed loads, the transverse distributed load and the distributed couple, respectively (for the sake of brevity
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
no forces are prescribed at the end sections), and the di2erential operator 1 0 0 0 1 0 H= 0 0 1 0
0
1911
(9)
9
is de!ned. 2.3. Constitutive law The linear generalised constitutive relation between generalised strain and stress can be written in the form: r = D”
(10)
after introducing the symmetric operator: 0 (ES)1 0 (EA)1 (EA)2 (ES)2 0 0 ; D= (ES) (EJ ) 0 (ES) 1 2 12 0
0
0
(11)
k
where the following quantities, containing the inertial properties of the cross-section and constitutive parameters, are introduced: E dA ; (12a) (EA) = A (ES) = E (y − y ) dA ; (12b) A (EJ ) = E (y − y)2 dA : (12c) A
E is the elastic modulus of the material of the -th part, (EJ )12 = (EJ )1 + (EJ )2 , k is the connection sti2ness. Since the reference points with ordinates y are assumed to coincide with the centroids of the relevant components, the static modules (12b) are zero and the matrix D is diagonal. 3. Eectiveness of the displacement-based nite element approach 3.1. Displacement-based 8nite element formulation The problem can be approached by assuming the displacements as unknown and obtaining the equilibrium conditions in terms of displacements. The variational condition (6) coincides with the stationary condition of the total potential energy functional. It is natural to seek the solution in the
1912
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
subspace of U = W1 × W2 × V = H 1 × H 1 × H 2 consisting of functions satisfying the kinematical boundary conditions. An approximated solution can be calculated using the displacement approach of the !nite element method. The displacement formulation of the !nite element method introduces a polynomial approximation of the displacement !eld at the interior of each element: u˜ = Nd ;
(13)
where N is the matrix of the shape functions and d is the vector of the nodal displacements parameters (hereinafter •˜ denotes the approximation of !eld •). The approximated strain vector ”˜u and stress vector are obtained from the approximated displacement !eld: ”˜u = DNd = Bd;
(14)
r˜u = D”˜u = DBd:
(15)
The projection of the variational condition (6) in the !nite dimensional space U˜ de!ned by d leads to the following algebraic problem: KPE d − f = 0; where
KPE = f=
0
Le
B T DB d z;
(17a)
(HN )T g d z
(17b)
0 Le
(16)
are the element sti2ness matrix and the element loads vector, respectively. 3.2. Slip locking in the 8DOF displacement-based element Di2erent choices can be made to approximate the displacement !eld in the interior of each element [3,9]. The simplest model providing the lowest regularity required involves an 8DOF !nite element where the nodal parameters are the following: d8TDOF = [w1L w2L vL ’L w1R w2R vR ’R ]
(18)
(the apex L denotes the left node and R denotes the right node of the !nite element). In this case, de=ections are described by a third-order polynomial ensuring the continuity of the !rst derivative and the axial displacements of the components are described by linear functions. The matrix N of the shape functions becomes $1 0 0 0 $2 0 0 0 ; 0 0 0 $ 0 0 0 $ N = (19) 1 2 0 %3 %4 0 0 %1 %2 0
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930 2000 300
1400 mm 2
250 1400 mm 2
22
Slab concrete Ec = 30500 MPa Slab rebars steel Es = 200000 MPa
14
1600
1913
50
Beam steel Es = 200000 MPa
500
Fig. 2. Cross-section and materials properties of tested beams.
where: 1 (1 − &); (20a) 2 1 $2 (z) = (1 + &); (20b) 2 1 %1 (z) = (2 − 3& + &3 ); (20c) 4 Le (20d) %2 (z) = (1 − & − &2 + &3 ); 8 1 (20e) %3 (z) = (2 + 3& − &3 ); 4 Le (−1 − & + &2 − &3 ) %4 (z) = (20f) 8 (& = 2z=Le − 1, & ∈ [ − 1; 1], z ∈ [0; Le ]) are the shape functions and Le is the element length. The sti2ness matrix is reported in the appendix. The numerical results show that the solution accuracy obtained by the 8DOF element strongly depends on the connection sti2ness and the error grows when the connection sti2ness becomes high. As an example, a simply supported beam is considered (L = 30 m, see Fig. 2 for details of the cross-section) and the diagrams of the curvature and slip are reported in Fig. 3 for two values of the shear connection sti2ness, measured by the non-dimensional parameter L [19], where 1 h2 1 2 =k (21) + + (EA)1 (EA)2 (EJ )12 $1 (z) =
and L is equal to the beam length in the case of a simply supported static scheme. From the numerical results, it is evident that for low-sti2ness shear connection ( L = 1), the 8DOF element is able to approximate both the exact curvature and slip, even when just a few elements are used. However, for high connection sti2ness ( L = 10), the curvature computed by the 8DOF element is forced to be constant along the element and an oscillatory trend arises in the slip description. In this regards, it can be observed that the variational Eq. (6) deriving from the displacement approach may be posed in the following form: L L L D0 D0 u · D0 uˆ d z + k D1 u · D1 uˆ d z = g · Huˆ d z ∀u; ˆ (22) 0
0
0
1914
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L 0.0
0.0
0.2
0.2
0.4
exact
0.8
2+2 8DOF
αL = 1
1+1 8DOF
0.6 0.8
2+2 8DOF
1.0 1.2
4+4 8DOF
0.4
0.6
χ/χmax
χ/χmax
L
1.0
1+1 8DOF
1.4
1.2
1.5
1.5
αL = 10
exact
2+2 8DOF exact 1+1 8DOF 2+2 8DOF
0.5
δ/δ max
1.0
0.0 -0.5
1+1 8DOF
0.5
δ/δ max
1.0
4+4 8DOF
0.0
-0.5
-1.0
-1.0
αL = 1
-1.5
αL = 10
exact
-1.5
Fig. 3. Locking problems of the 8DOF element.
where the term depending on the sti2ness parameter k has been separated, the di2erential operator has been split in 9 0 0 0 (23a) D0 = ; 0 9 0 0 −92 D1 = [ − 1
1
h9]
(23b)
and D0 is constituted by the !rst three rows and columns of D. In the limit case where k → ∞, the solution must lie in the subspace U0 of the functions satisfying the condition [11]: L
0
D1 u · D1 u d z = 0:
(24)
This constraint induces a linear relation between the displacements !elds ensuring that the slip is null, i.e. lim = lim (w2 − w1 + hv ) = 0: (25) k →∞
k →∞
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1915
So that one displacement !eld, as an example w1 , can be deduced from the other displacement !elds and U0 can be assumed as U0 = W2 × V . In this case, as in other multi!eld problems [12], troubles arise in applying the constraint in the !nite dimensional approximating space. In fact, the !nite dimensional spaces U˜ = W˜ 1 × W˜ 2 × V˜ generated by the polynomial sequences reduce, in the limit case k → ∞, to U˜ 0 = W˜ 2 × V˜ 0 where V˜ 0 ⊂ V˜ . In other words, the constraint does not simply induce a linear relation between the unknowns but also reduces the dimension of the space V˜ in which the transversal displacements are sought. For the sake of brevity, proof of what asserted is given for the particular case of the 8DOF element; the demonstration, however, can be generalized in a simple way in order to consider other !nite elements. In the 8DOF element, the slip locking problem derives from the di2erent dimensions of the spaces involved in the expression of the coupling term , where the term provided by the axial displacements lies in the two-dimensional space of !rst-order polynomials, while the term furnished by the transverse displacement lies in the three-dimensional space of the second-order polynomials. In fact, the substitution of the 8DOF displacement !eld, whose shape functions are grouped in matrix (19), in the slip expression (4) leads to (z) = a
1+& 1−& 3h c(1 − &2 ) +b + 2 2 2Le
(26)
or the equivalent form, obtained from the previous equation by collecting the terms of the same degree: h 3h 2 1 1 a + b + 3 c + (a − b)& − (z) = c& ; (27) 2 Le 2 2Le where a = −w1R + w2R + h’R ;
(28a)
b = −w1L + w2L + h’L ;
(28b)
Le (’L + ’R ): (28c) 2 Condition (25) leads to three equations ensuring that the second-order polynomial describing the slip is null. The !rst two (i.e. a = 0 and b = 0) involve both parameters describing axial and transverse displacement and express the kinematical restraint provided by a rigid shear connection: c = v R − vL −
w1L − w2L = h’L ;
(29a)
w1R − w2R = h’R ;
(29b)
while the third equation (i.e. condition on the quadratic term c = 0) involves only parameters related to transverse displacement: Le (’R + ’L ) = vR − vL ; 2
(29c)
1916
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
such condition does not re=ect the physical relation arising from axial and transverse displacement but introduces a further constraint on the de=ection DOF. The substitution of this last condition (29c) in the expression of the 8DOF de=ection v(z) leads to the following reduced form: v(z) = vL + ’L Le 18 (3 + 2& − &2 ) + ’R Le 18 (1 + 2& + &2 );
(30)
where the third-order terms, proportional to the parameter c (28c), disappear and v(z) is described by a second-order polynomial only; in other words, the e2ective DOF are reduced. Consequently, the curvature is constrained to be a constant along the element. In the limit condition k → ∞, this a2ects the slip description producing a spurious oscillatory trend in the solution [12]. The situation described occurs every time the contribution in slip at the interface deriving from axial displacements and from de=ection are described by means of polynomials of di2erent order (i.e. w and v are polynomials of di2erent order). This leads to one or more conditions involving only axial or transverse displacement !elds that cannot express the kinematical restraint provided by a rigid shear connection but introduces further constraint on the axial or on the de=ection !elds, respectively. Obviously, the poorer the element, the more evident the consequent locking will be. 3.3. Re8ned displacement elements By starting from the 8DOF element, the problem can be overcome by introducing additional internal nodes [3,9]. If an internal with a DOF is introduced for each of the two axial displacements, the 10DOF element is obtained [3,9]. In this element, both contributions to slip description, deriving from axial and transverse approximations, are second-order polynomials. Other elements based on the choice of calibrated internal nodes can be de!ned. If a further node, with one transverse and one rotational DOF, is used in the transverse !eld description (v is now a !fth-order polynomial), then three internal nodes are needed for each of the two axial displacement !elds. The resulting element (16DOF) describes the slip with a fourth-order polynomial [3,9]. 3.4. Comparisons between displacement elements This section reports some numerical results in order to compare the behaviour of the approximate solution with respect to the connection sti2ness and the convergence rate performed by the three previous !nite elements, 8DOF, 10DOF and 16DOF. Firstly, a direct comparison between the solution obtained with the 8DOF, 10DOF and 16DOF !nite elements is reported in Fig. 4, where the same problem previously analysed in Fig. 2 is considered and a two-element discretization is adopted. It is evident that the accuracy drastically improves for high connection sti2ness by using the 10DOF or the 16DOF !nite elements. In fact, elements with internal nodes furnish almost the same accuracy both for low and high sti2ness levels. For a more detailed analysis of the behaviour of di2erent elements, a cantilever with length L=2 is considered. In this case, the exact (exponential) solution di2ers more remarkably from the solutions provided by the shape functions and the approximation errors are larger and more evident. The same cross-section of the previous applications is considered and the connection sti2ness is described by the parameter L earlier introduced.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L
1917
L
0.0
0.0
0.2
0.2
10DOF 0.6
10DOF
0.8
exact
16DOF exact
1.0
αL= 1
α L = 10 1.2
1.5
1.5
8DOF 10DOF 16DOF exact
0.5 0.0
1.0 0.5
δ/δ max
1.0
δ/δ max
8DOF 0.6 0.8
16DOF
1.0 1.2
0.4
8DOF
χ/χ max
χ/χ max
0.4
10DOF
16DOF
0.0
exact
-0.5
-0.5
-1.0
-1.0
8DOF
-1.5
αL = 1
-1.5
α L = 10
Fig. 4. Comparisons between displacement elements: curvature and slip in the simply supported beam.
Fig. 5 describes the curvature obtained for four di2erent values of the connection sti2ness parameter ( L = 1, 7, 20, 50). This comparison intends to analyse the sensitivity with respect to L and to present the distribution of the local error along the whole beam. Discretisation adopts only one element; it is evident that the 16DOF element has more DOF describing the curvature (third-order polynomial) and potentially can better !t the exact solution with respect to both the 8DOF and 10DOF elements furnishing a linear description of the curvature. It can be observed that 10DOF and 16DOF elements present solutions with the same order of accuracy for all the connection sti2ness considered. It is also interesting to observe the di2erences between the 8DOF and the 10DOF elements, for which the curvature is described by the same polynomial order. Fig. 6 reports results concerning the interface slip. In this case, the shape of the solution change very considerably by varying from low to high connection sti2ness values. The 8DOF and 10DOF elements can provide a parabolic description while the 16DOF element can describe the interface slip by means of a fourth-order polynomial. The 16DOF element furnishes satisfactory results at all the connection sti2ness levels, while the 8DOF element local error is larger than the local error deriving from the 10DOF element.
1918
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2
L/2 1.2
1.2
16DOF
1.0
0.8
10DOF
0.6
8DOF
0.4 0.2
χ/χ max
χ/χ max
0.8
-0.2
αL = 1
-0.4
-0.4
1.0
1.0
16DOF
αL = 7
16DOF
0.8
10DOF
0.6
8DOF
0.4
χ/χ max
0.8
χ/χ max
exact
1.2
1.2
10DOF
0.6
8DOF
0.4 0.2
0.2
exact
0.0
-0.4
8DOF
0.4
0.0
0.0
-0.2
10DOF 0.6
0.2
exact
-0.2
16DOF
1.0
αL = 20
0.0 -0.2 -0.4
exact
αL = 50
Fig. 5. Comparisons between displacement elements: curvature in the cantilever.
The global error is analysed in Figs. 7–9 and in Tables 1–4. The error is evaluated both for the simply supported static scheme and the cantilever. Fig. 7 and Tables 1 and 2 report the energy error e=
u − u˜h ;
u
(31)
where
u =
0
L
DDu · Du d z:
(32)
Fig. 8 and Tables 3 and 4 describe the results obtained for the maximum displacement (mid-span for the simply supported beam and end displacement for the cantilever). The results concerning the 16DOF element are not showed for the simply supported beam because they are out of scale (the error is too small to be depicted). The error drastically reduces by adopting elements with internal nodes (diagrams are traced in a logarithmic scale). It can be observed that the error always increases when the connection sti2ness increases only in the case of the 8DOF element.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2
L/2 1.2
1.2
αL = 2
0.8
δ/δ max
0.6 0.4
8DOF
0.8
exact 16DOF
exact 0.6
16DOF 0.4
8DOF
0.2
αL = 7
1.0
1.0
δ/δ max
1919
0.2
10DOF
10DOF
0.0
0.0
1.5
1.5
8DOF
8DOF
1.0
1.0
exact 0.5
10DOF
0.0
δ/δmax
δ/δmax
0.5
16DOF
10DOF
-0.5
-0.5
16DOF
-1.0
-1.0 -1.5
exact 0.0
α L = 20
α L = 50
-1.5
Fig. 6. Comparisons between displacement elements: slip in the cantilever. 100.000%
100.000% 1 8DOF
1+1 8DOF
2 8DOF
10.000%
2+2 8DOF
energy error
energy error
10.000% 1.000%
1+1 10DOF 4+4 8DOF
0.100% 2+2 10DOF
0.010%
1 10DOF
1.000%
4 8DOF 2 10DOF
0.100%
1 16DOF
0.010%
4 10DOF 2 16DOF
4+4 10DOF
0.001%
1
10
αL
100
0.001%
1
10
100
αL
Fig. 7. Comparisons between displacement elements: energy error versus connection sti2ness.
Finally, Fig. 9 describes the convergence rate by reporting the energy error versus the total DOF (the internal nodes DOF are computed in the total DOF). As in the previous case, the results for the 16DOF element are not showed for the simply supported beam because they are too small. The di2erent performances of the three elements are, once again, evident.
1920
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930 100.000%
100.000%
1 8DOF
1+1 8DOF
v
10.000%
2+2 8DOF
nodal error
nodal error
10.000%
4+4 8DOF
1.000% 0.100%
4 8DOF
1.000%
1 10DOF
0.100%
2 10DOF
1+1 10DOF
0.010%
2 8DOF
v
4 10DOF
0.010%
1 16DOF
2+2 10DOF
0.001%
2 16DOF
0.001% 1
10
100
1
10
αL
100
αL
Fig. 8. Comparisons between displacement elements: nodal error versus connection sti2ness.
100.000%
100.000% 8DOF αL = 20
8DOF αL = 20
10.000%
8DOF αL = 7
energy error
energy error
10.000% 1.000%
10DOF αL = 7
0.100%
8DOF αL = 7
1.000%
10DOF αL = 7 10DOF αL = 20
0.100%
16DOF αL = 20
0.010%
0.010% 10DOF αL = 20
total DOF
0.001% 0
5
10
15
20
25
0
30
16DOF αL = 7
total DOF
0.001%
5
10
15
20
25
30
Fig. 9. Comparisons between displacement elements: energy error versus total DOF.
Table 1 Energy error in the simply supported beam under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1+1 el. 2+2 el. 4+4 el.
1.06 0.07 0.01
2.24 0.43 0.10
4.68 1.09 0.27
7.43 1.79 0.45
10.28 2.50 0.62
13.30 3.26 0.79
14.81 3.71 0.88
16.42 4.52 1.04
16.81 4.93 1.21
10DOF
1+1 el. 2+2 el. 4+4 el.
1.03 0.07 0.00
0.88 0.06 0.00
0.86 0.06 0.00
0.89 0.06 0.01
0.93 0.06 0.01
0.99 0.07 0.01
1.03 0.07 0.01
1.06 0.07 0.01
1.07 0.07 0.01
16DOF
1+1 el.
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1921
Table 2 Energy error in the cantilever under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1 el. 2 el. 4 el.
3.13 0.21 0.02
6.24 0.65 0.11
10.25 1.36 0.27
14.46 2.23 0.47
19.52 3.45 0.74
27.38 5.70 1.20
33.26 7.68 1.59
41.61 11.44 2.48
43.99 17.31 8.75
10DOF
1 el. 2 el. 4 el.
3.10 0.20 0.02
4.85 0.33 0.03
6.07 0.45 0.03
6.55 0.55 0.04
6.31 0.63 0.06
5.31 0.69 0.07
4.51 0.65 0.09
3.36 0.41 0.08
2.96 0.25 0.04
16DOF
1 el. 2 el. 4 el.
0.00 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.00
0.01 0.00 0.00
0.04 0.00 0.00
0.09 0.00 0.00
0.13 0.01 0.00
0.11 0.02 0.00
0.04 0.01 0.00
Table 3 Nodal de=ection error in the simply supported beam under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1+1 el. 2+2 el. 4+4 el.
0.02 0.01 0.00
1.19 0.35 0.09
3.80 1.00 0.25
6.90 1.69 0.42
10.35 2.38 0.59
14.39 3.10 0.74
16.63 3.53 0.82
19.18 4.37 0.95
19.83 4.83 1.12
10DOF
1+1 el. 2+2 el. 4+4 el.
0.05 0.00 0.00
0.22 0.01 0.00
0.25 0.01 0.00
0.21 0.01 0.00
0.13 0.00 0.00
0.05 0.00 0.00
0.02 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.00
16DOF
1+1 el.
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
Table 4 Nodal de=ection error in the cantilever under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1 el. 2 el. 4 el.
0.11 0.01 0.00
2.19 0.40 0.09
5.63 1.07 0.25
9.61 1.84 0.43
14.38 2.74 0.63
21.07 4.06 0.89
25.52 5.09 1.08
31.34 7.10 1.49
32.91 8.02 1.85
10DOF
1 el. 2 el. 4 el.
0.08 0.01 0.00
0.72 0.05 0.00
1.25 0.10 0.01
1.52 0.15 0.01
1.49 0.19 0.02
1.09 0.21 0.02
0.75 0.19 0.03
0.24 0.09 0.02
0.06 0.03 0.01
16DOF
1 el. 2 el.
0.00 0.00
0.00 0.00
0.00 0.00
0.01 0.00
0.02 0.00
0.03 0.00
0.05 0.00
0.04 0.01
0.02 0.00
1922
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
4. Eectiveness of nite element mixed approaches 4.1. Three-8eld 8nite element mixed formulation The problem can be approached by assuming the displacement, strain and stress !elds as unknowns and solved using the compatibility condition, balance condition and the generalised constitutive law. An approximated solution can be obtained using a three-!eld formulation of the !nite element method. The method introduces polynomial approximations of the displacement !eld, of the strain !eld and of the stress !eld at the interior of each element: u˜ = Nd ;
”˜ = Ee; r˜ = Ss;
(33a,b,c)
where N , E and S are the matrices of shape functions and d , e and s are the vectors of displacement, strain and stress nodal parameters. The projection of the weak form of the balance, compatibility and constitutive equations (details in [9]) in the !nite dimensional space de!ned by d , e, s leads to the following solving system: BV T s − f = 0;
(34a)
V − EV T s = 0; De
(34b)
V − Bd V = 0; Ee
(34c)
where BV = EV = DV =
Le
0
Le
0
0
Le
S T B d z;
(35)
S T E d z;
(36)
E T DE d z:
(37)
Since the stress and strain !elds can be inter elements discontinuous, the system of equations (34) can be simpli!ed performing a static condensation of the stress and strain nodal parameters [9]. The resulting system has the same shape of the solving system derived in displacement approach (16) and can be easily used in a displacement-based !nite element program. 4.2. Reduction of slip locking in mixed elements The authors have already presented and tested three-!eld mixed elements for composite beams with deformable shear connection [9] based on the locking-free displacement !eld of the 10DOF
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1923
Table 5 Shape functions polynomials degree Element
8DOF HW011 HW010
Displacements
Strains
Stresses
w
v
N
M
fs
1 1 1
3 3 3
0 0 0
1 1 1
2 1 0
0 0 0
1 1 1
2 1 0
displacement element. In this paper, the e2ects of slip locking in mixed elements will be illustrated by introducing mixed elements based on the 8DOF displacement !eld. Once the displacement !eld is selected, the Limitation Principles [20,21] must be considered for a proper de!nition of the strain and stress approximations. Roughly speaking, the Limitations Principles state that the displacement approach and the three-!eld mixed approach become identical if the strain approximation ”˜ and the stress approximation r˜ are equal or richer than the strain ”˜u and the stress r˜u calculated from the approximation of the displacement !eld. Di2erent mixed elements were tested by the authors using the same shape functions for conjugated strains and stresses (i.e. the matrices E and S are equal). The tests showed that it is necessary to reduce approximation of the slip (and the interface force) in order to reduce locking. As a result, four choices are possible: constant or linear slip (interface force fs ) combined with constant axial strains (axial forces N ) and constant or linear curvature (moment M12 ). However, elements with constant curvature have a too poor behaviour if compared to the element with linear curvature, consequently only the results for elements with linear curvature are discussed in this paper. In Table 5, a comparisons between the shape functions of the 8DOF element and of the three-!eld mixed element with linear curvature and reduced slip (HW010 and HW011) is illustrated. The strain and stress polynomial degree of the displacement-based element are indicated in italics because they are not independent from the displacement shape functions but derived from them (using compatibility and linear constitutive law). It is worth saying that if the displacement–stress two-!eld mixed formulation is adopted and if the same representations is used for the stress both in the three-!eld and two-!eld formulations, then the same results are obtained as a consequence of the Limitation Theorems (two-!eld and three-!eld formulations coincide in the case considered). Similar considerations apply to the displacement– strain two-!eld mixed formulation, previously introduced by the authors [22] before the more general three-!eld formulation. 4.3. Comparisons between mixed and displacement elements This section reports some numerical results in order to compare the behaviour of the approximate solution, furnished by the !nite elements 8DOF, HW010 and HW011, with respect to the variation of the connection sti2ness. As in the previous case, a cantilever with length L=2 is considered. Fig. 10 describes the curvature obtained for four di2erent values of the connection sti2ness parameter ( L = 1, 7, 20, 50). For low sti2ness values ( L = 1 and L = 7), the three elements give
1924
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2
L/2 1.2
1.2 1.0
exact
1.0
HW011 HW010
0.4
χ/χmax
χ/χmax
8DOF
0.6
0.2
8DOF
0.0
0.0 -0.2
1.0
αL= 1
-0.4
0.8
exact
8DOF HW011
χ/χ max
0.6
0.4
8DOF
0.4 0.2
HW011
0.0
0.0
-0.4
αL = 7
1.0
exact
0.6
-0.2
HW010
1.2
0.8
0.2
HW011
0.4 0.2
1.2
χ/χmax
0.6
-0.2 -0.4
exact
0.8
0.8
-0.2
αL = 20
HW010
-0.4 -0.6
HW010 αL = 50
Fig. 10. Comparisons between mixed and displacement elements: curvature in the cantilever.
the same results (both displacement element and mixed elements describe the curvature with a linear function). When sti2ness increases, it is evident that the curvatures computed by the 8DOF and by the HW011 are forced to be constant along the beam axis, while the HW010 linear curvature does not degenerate into a constant. Fig. 11 describes the interface slip obtained for the same four di2erent values of the connection sti2ness parameter ( L = 1, 7, 20, 50). For low sti2ness values, the 8DOF element is locally more accurate than the mixed elements, but when sti2ness increases the HW010 element furnishes lower local errors. Notice that the mixed elements cannot satisfy the boundary conditions. Finally, in Figs. 12 and 13, the same problem is analysed using a two-element discretisation: it can be observed that a more re!ned mesh reduces the local errors but do not eliminate the described problems. In fact, the curvature trends for the 8DOF and HW011 elements are forced to be a constant along each element, while slip oscillations once again arise. The convergence curves for the mixed elements, for the sake of brevity not reported here, are very close to the curves of the 8DOF element.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2 1.4
L/2 1.4
αL = 1
1.2
1.0
HW010
δ/δ max
δ/δmax
0.8
exact
αL = 7
8DOF
1.2
8DOF
1.0
0.6
1925
HW010
0.8
HW011
0.6
exact
HW011 0.4
0.4
0.2
0.2
0.0
0.0
2.5
1.5
HW011
HW011
2.0
HW010
0.5
1.5
δ/δmax
δ/δmax
1.0
exact 0.0
HW010
1.0 0.5 0.0
exact
-0.5
-0.5 -1.0
8DOF αL = 20
8DOF
-1.0 -1.5
αL = 50
Fig. 11. Comparisons between mixed and displacement elements: slip in the cantilever.
5. Conclusions In this paper, the authors discuss problems that may take place in the !nite element analysis of composite beams with deformable shear connection, where slip occurs between beam components while contact is preserved. A coupling between the two axial displacement !elds and the =exural displacement !eld arise in this beam model due to the deformable interface. If the axial and =exural displacement shape functions are not consistent, then locking problems occur for high values of the connection sti2ness (slip locking). Slip locking is comparable to membrane and shear locking that may develop in common beam elements, and it mostly a2ect the description of the curvature and interface slip. This is a problem of particular interest in the non-linear analysis of composite beams, because the non-linear behaviour of the shear connection is characterised by a constitutive law with a sti2ness varying from very a high value when the slip is close to zero, to a very low value when the connection is close to failure. Since the distribution of interface slip varies greatly along the beam, the solution obtained, which depends strongly on the connection sti2ness, is not accurate. A strategy to avoid locking by means of a calibrated choice of the
1926
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
αL/2
L/2
1.2
1.2
exact
1.0
1.0
0.8
exact
0.8
8DOF
HW011 HW010
0.4
χ/χmax
χ/χmax
8DOF 0.6
0.2
HW010
0.0
αL = 1
-0.2
1.2 1.0
HW011
0.4 0.2
0.0 -0.2
0.6
αL = 7
1.2
exact
1.0
0.8
exact
0.8
8DOF
HW011
0.4
HW010
0.2
0.6
HW011 HW010
0.4 0.2
0.0 -0.2
χ/χmax
χ/χmax
8DOF 0.6
0.0
αL = 20
-0.2
αL = 50
Fig. 12. Comparisons between mixed and displacement elements: curvature in the cantilever (two-element mesh).
displacement shape functions is then proposed and a comparison between di2erent displacement-based elements, in terms of local error, global error and convergence rate, is illustrated. The authors also analyse the e2ect of locking on mixed elements. Starting from the displacement !eld of the element with shape functions of the lowest degree required by the problem (8DOF element), three-!eld mixed elements are introduced. Since the 8DOF displacement element is a2ected by locking, the applications intend to clarify if the mixed elements are still a2ected by slip locking. The numerical tests performed show that a mixed element based on a locking-su2ering displacement element can still be a2ected by locking. It is possible to reduce locking by lowering the degree of the shape functions for the slip and interface force, but this leads to poor description of these entities. Locking-free displacement elements, formulated accordingly to the indications of consistence between axial and =exural !elds, and mixed elements based on locking-free displacement shape functions, allow to obtain more accurate and reliable results in both linear and non-linear analyses.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2 1.2
L/2 1.4
αL = 1
8DOF 1.0
exact
δ/δmax
δ/δmax
αL = 7
HW011
1.2
1.0 0.8 0.6
8DOF
0.8
HW010
0.6
exact
HW010
0.4
1927
0.4
HW011
0.2
0.2 0.0
0.0
3.0
2.0
HW011 2.5
1.5
HW010 exact
δ/δmax
δ/δmax
1.0 0.5
1.5
HW010
1.0 0.5
0.0
0.0 -0.5 -1.0
HW011
2.0
8DOF
-0.5
αL = 20
-1.0
exact 8DOF
αL = 50
Fig. 13. Comparisons between mixed and displacement elements: slip in the cantilever (two-element mesh).
Appendix A A.1. The 8DOF element sti>ness matrix
K8DOF = where
KLL
KLR
T KLR
KRR
EA1 L 0 KLL = 0 0
;
(A.1)
0
0
EA2 L
0
0 0
EJ12 L3 EJ12 6 2 L
12
L 0 3 L − 0 3 +k h EJ12 6 2 2 L hL EJ12 4 − L 12
L − 3 L 3 h − 2 hL 12
h 2 h − 2 6h2 5L h2 10
hL − 12 hL 12 ; 2 h 10 2h2 L 15
(A.2)
1928
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L 0 0 0 3 L EA2 − 0 0 3 L +k h EJ12 EJ12 − 0 12 3 −6 2 2 L L hL EJ12 EJ12 0 −6 2 4 − L L 12 L EA1 − 0 0 0 L 6 L EA − 2 0 0 0 − 6 L +k KLR = h EJ12 EJ12 0 0 −12 6 3 2 2 L L hL EJ12 EJ12 0 0 −6 2 2 L L 12 EA1 L 0 KRR = 0 0
L − 3 L 3 h 2 hL 12
L 6 L 6 h − 2 hL − 12 −
hL − 12 hL 12 ; 2 h − 10 2h2 L 15 h hL − 2 12 h −hL 2 12 : 6h2 h2 − 5L 10 h2 h2 L − − 10 30
h − 2 h 2 6h2 5L h2 − 10
(A.3)
(A.4)
A.2. The HW010 mixed element sti>ness matrix
KV HW010 = where
KLL
KLR
T KLR
KRR
EA1 L 0 KLL = 0 0 EA1 L 0 KRR = 0 0
;
(A.5)
0
0
EA2 L
0
0 0
EJ12 L3 EJ12 6 2 L
12
0
0
EA2 L
0
0 0
EJ12 12 3 L EJ12 −6 2 L
L h L 0 − 4 4 2 L h L 0 − − 0 +k 4 4 2 ; EJ12 2 h h h 6 2 − 0 L 2 2 L EJ12 4 0 0 0 0 L L L h 0 − − 0 4 4 2 L L h 0 − 0 +k 4 4 2 ; EJ12 2 h h h −6 2 0 − L 2 2 L EJ12 4 0 0 0 0 L 0
(A.6)
(A.7)
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
KLR =
−
EA1 L 0
−
0
0
0
EA2 L
0
0
0
0
0
0
EJ12 L3 EJ12 −6 2 L
−12
EJ12 L2 EJ12 2 L 6
L 4 L − +k 4 h 2 0
L 4 L 4 h − 2 0
−
h 2 h 2 h2 − L 0 −
0 0 : 0 0
1929
(A.8)
References [1] D.J. Oehlers, M.A. Bradford, Elementary Behaviour of Composite Steel and Concrete Structural Members, Butterworth-Heinemann, London, 2000. [2] N.M. Newmark, C.P. Siess, I.M. Viest, Tests and analysis of composite beams with incomplete interaction, Proc. Soc. Exp. Stress Anal. 9 (1) (1951) 75–92. [3] A. Dall’Asta, A. Zona, Non-linear analysis of composite beams by a displacement approach, Comput. Struct. 80 (27–30) (2002) 2217–2228. [4] M.R. Salari, E. Spacone, P.B. Shing, D.M. Frangopol, Non-linear analysis of composite beams with deformable shear connectors, J. Struct. Eng. ASCE 124 (10) (1998) 1148–1158. [5] M.R. Salari, E. Spacone, Finite element formulations of one-dimensional elements with bond-slip, Eng. Struct. 23 (7) (2001) 815–826. [6] A. Ayoub, F.C. Filippou, Mixed formulation of nonlinear steel–concrete composite beam element, J. Struct. Eng. ASCE 126 (3) (2000) 371–381. [7] A. Ayoub, A two-!eld mixed variational principle for partially connected composite beams, Finite Elem. Anal. Des. 37 (2001) 929–959. [8] B. Cas, F. Saje, M. Saje, I. Planinc, Nonlinear !nite element analysis of composite planar frames with interlayer slip, in: B.H.V. Topping, Z. Bittnar (Eds.), Proceedings of the Sixth International Conference on Computational Structures Technology, Civil-Comp Press, Stirling, Scotland, 2000, Paper 135. [9] A. Dall’Asta, A. Zona, Three-!eld mixed formulation for the non-linear analysis of composite beams with deformable shear connection, Finite Elem. Anal. Des. 40 (4) (2004) 425– 448. [10] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, New York, 2000. [11] D. Chapelle, Reliability of !nite element methods for thin shells, in: B.H.V. Topping (Ed.), Computational Mechanics for the Twenty-First Century, Saxe-Coburg Publications, Edinburgh, 2000, pp. 99–108. [12] G. Prathap, B.P. Naganarayana, Stress oscillations and spurious load mechanisms in variationally inconsistent assumed strain formulations, Int. J. Numer. Methods Eng. 33 (1992) 2181–2197. [13] M.A. Cris!eld, Non-linear Finite Element Analysis of Solids and Structures, Vol. 1: Essentials, Wiley, New York, 1991. [14] L. Yunhua, Explanation and elimination of shear locking and membrane locking with !eld consistency approach, Comput. Methods Appl. Mech. Eng. 162 (1998) 249–269. [15] J.N. Reddy, On locking-free shear deformable beam !nite elements, Comput. Methods Appl. Mech. Eng. 149 (1997) 113–132. [16] I.J. Blaauwendraad, Realistic analysis of reinforced concrete framed structures, Heron 18 (4) (1972) 1–31. [17] A.K. Gupta, P.S. Ma, Error in eccentric beam formulation, Int. J. Numer. Methods Eng. 11 (1977) 1473–1483. [18] M.A. Cris!eld, The eccentricity issue in the design of plate and shell elements, Commun. Appl. Numer. Methods 7 (1991) 47–56. [19] U.A. Girhammar, D. Pan, Dynamic analysis of composite members with interlayer slip, Int. J. Solids Struct. 30 (6) (1993) 797–823.
1930
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
[20] B.M. Fraeijs De Veubeke, Displacement and equilibrium models in the !nite element method, in: O.C. Zienkiewicz, G.S. Hollister (Eds.), Stress Analysis, Wiley & Sons, London, 1965, pp. 145–196 Reprinted in: Int. J. Numer. Methods Eng. 52 (2001) 287–342. [21] H. Stolarski, T. Belytschko, Limitation Principles for mixed !nite elements based on the Hu-Washizu variational formulation, Comput. Methods Appl. Mech. Eng. 60 (1987) 195–216. [22] A. Dall’Asta, A. Zona, Finite elements for the analysis of composite beams with interlayer slip, Proceedings of the XVIII Congresso C.T.A., Vol. 2, Venezia, Italy, 26–28 September 2001, pp. 365–374.
Slip locking in !nite elements for composite beams with deformable shear connection Andrea Dall’Astaa;∗ , Alessandro Zonab a
Dipartimento di Progettazione e Costruzione dell’Ambiente, Universita di Camerino, Viale della Rimembranza, 63100 Ascoli Piceno, Italy b Istituto di Scienza e Tecnica delle Costruzioni, Universita Politecnica delle Marche, Via Brecce Bianche, 60131 Ancona, Italy Received 15 September 2002; received in revised form 21 October 2003; accepted 12 January 2004
Abstract In this paper, the authors discuss problems that may occur in the !nite element analysis of composite beams with deformable shear connection, where slip arises between beam components while contact is preserved. E2ects of locking problems, comparable to the membrane and shear locking that may develop in !nite elements for common beams, are illustrated. These locking problems mostly a2ect the description of the curvature and interface slip. A strategy to avoid locking by means of a calibrated choice of the displacement shape functions is then proposed and a comparison between di2erent displacement-based elements, in terms of local error, global error and convergence rate, is presented. The authors also illustrate the e2ects of slip locking on mixed !nite elements and comparisons between displacement-based and mixed elements are shown. ? 2004 Elsevier B.V. All rights reserved. Keywords: Composite beams; Shear connection; Locking; Slip oscillations; Mixed formulation; Beam !nite elements
1. Introduction Composite beams, made by two components joined by shear connectors to form an interacting unit, are frequently used in structural engineering. Common examples are steel-concrete composite decks for bridges construction. The behaviour of composite beams is in=uenced by the type of connection between the two parts. Flexible shear connectors permit the development of partial composite action [1] and, for a realistic analysis, the structural model must account for the interlayer slip between the components. Therefore, a simple analytical formulation that captures the salient features of the global ∗
Corresponding author. Tel.: +39-071-220-4551; fax: +39-071-220-4576. E-mail address: [email protected] (A. Dall’Asta).
0168-874X/$ - see front matter ? 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.!nel.2004.01.007
1908
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
x G1
z
y1
w1 k y2
h
w2 k
G2
y
vj v'
y vj
v' δ
Fig. 1. Newmark kinematical model.
behaviour (bending behaviour and interface slip) is very useful both in the analysis and design of this structural system. If the study is limited to the in-plane bending behaviour, i.e. symmetric cross-sections are considered and no torsion and out-of-symmetry-plane bending occur, the analysis of this type of structure can be based on the Newmark kinematical model [2]: the Euler–Bernoulli linear beam theory (small displacements, rotations and strains) is used to model the two parts of the composite beam; the e2ects of the deformable shear connection are accounted for by using an interface model with distributed bond, preserving the contact between the components (Fig. 1). The interface slip, computed from the di2erence of small displacements of the two parts at the interface, is accordingly small. Once a kinematical model is de!ned for the problem, numerical solutions can be obtained by means of the !nite element method. The !nite element analysis of composite beams makes it possible to obtain useful information about the structural behaviour but requires an eHcient and reliable element. Quite a number of works on this topic can be found in the literature. Di2erent displacement-based !nite element models have been adopted (see [3] and the references quoted there). Finite elements based on the displacement approach have a simple formulation but their behaviour is not always satisfactory. More recently, models that attempt to overcome the limitations of the displacement-based formulation have been proposed: Salari et al. [4], Salari et Spacone [5], Ayoub et Filippou [6], Ayoub [7], Cas et al. [8] and Dall’Asta et Zona [9]. All of these !nite element models indicate the great interest in solving the convergence problems that a2ect the analysis of composite beams with deformable shear connection. This paper aims at clarifying some of the aspects involved in the !nite element solution. For typical structural problems, the best approximation property of the !nite element method [10] is a signi!cant guarantee of reasonable behaviour even when very crude meshes are employed. Generally, when two or more !elds are coupled in a structural model, the constraint between !elds may cause troubles when applied to !nite dimensional spaces, as those used in !nite element approximations. The accuracy of the solution depends on some characteristic parameters involved in the coupled terms. For limit values, such as zero or in!nity, further relations between unknowns reduce the dimension of the space in which the solution is sought [11,12]. In some cases, the dimension goes to zero and the model completely “locks”. In general, a sti2ener response and spurious strains
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1909
can be observed when the phenomenon occurs. Examples of these problems for beam elements include membrane locking [13,14] that may a2ect curved beam elements by varying the curvature; shear locking [14,15] that may develop in Timoshenko beam elements by varying the shear sti2ness; the eccentricity issue [16–18] that may involve the ordinary Euler–Bernoulli beam by varying the origin of the reference system. In all of these cases, the strain parameters involved are contributed by di2erent displacement !elds. A similar problem may a2ect composite beam elements. In fact a coupling arises between the =exural displacement !eld and the axial displacements !elds of the two components due to the deformable shear connection. As shown in this paper, if the approximations of the axial and =exural displacement !elds are not consistent, then the error of the !nite element solution depends strongly on the connection sti2ness, slip oscillation and poor approximation of the curvature occur for high values of the connection sti2ness. This is a problem of particular interest in the non-linear analysis of composite beams [3] because the non-linear behaviour of the shear connection is characterised by a constitutive law with a sti2ness varying from very a high value when the slip is close to zero, to a very low value when the connection is close to failure. Since the distribution of interface slip varies greatly along the beam, the solution obtained, which depends strongly on the connection sti2ness, is not accurate. In this paper, the authors analyse the e2ectiveness of !nite elements based on the Newmark kinematical model [2] in linear analyses. The inconveniences that arise in the displacement-based element, with shape functions of the lowest degree required by the problem (8DOF element), are illustrated. A strategy for avoiding slip locking by means of a calibrated choice of axial and transverse displacement shape functions is proposed. The authors also illustrate the e2ects of slip locking in mixed elements based on the 8DOF displacement !eld approximation. The proposed comparisons between displacement and mixed elements are intended to clarify whether a mixed formulation is a2ected by slip locking or if it can be used to reduce or eliminate slip locking problems.
2. Kinematical model and main equations 2.1. Kinematical model In the reference con!guration, the composite beam occupies the region V = A × [0; L] generated by translating its symmetrical cross-section A along a rectilinear axis, orthogonal to the cross-section. The cross-section is divided in two parts A = A1 ∪ A2 . We introduce an orthonormal reference frame {O; X; Y; Z}, where Z is orthogonal to A. The unit vectors parallel to the axis X , Y , Z are denoted by i, j, k, respectively. The connection is assumed to lie on a plane that is orthogonal to the Y -axis and permits only discontinuities parallel to the beam axis. It is also assumed that YZ is a symmetry plane for the problem and shear deformation cannot occur; in other words, every component deforms like plane Kirchho2 rod, and consequently the displacement !eld has the following representation: u (y; z) = v(z)j + [w (z) + (y − y)v (z)]k
in A × [0; L] ( = 1; 2);
(1)
where the prime denotes the derivatives with respect to the z co-ordinate and u ( = 1; 2) can be deduced from the three scalar functions v and w describing the displacement of the centroid of A
1910
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
in Y and Z directions (y is the ordinate of the centroid of A , see Fig. 1). The functions furnishing the displacement !eld can be collected in the vector uT = [w1
w2
v]:
(2)
2.2. Strain and stress parameters The kinematical assumptions produce a reduced set of deformations in the beam. More precisely, only the axial strain z (y; z) = w + (y − y)v = (z) + (y − y )(z)
in A × [0; L] ( = 1; 2)
(3)
exists ( = w is the axial strain evaluated at the centroid of A and = −v is the curvature), while the slip at the interface is expressed by the following relation: (z) = w2 (z) − w1 (z) + hv (z)
in [0; L];
(4)
where h=y2 −y1 is the distance of the two components centroids (Fig. 1). The parameters describing the strain measures can be collected in the vector ” which is related to the displacements u by means of the di2erential operator D (compatibility equations): 9 0 0 1 w1 0 2 0 9 w = Du: ”= (5) = 2 2 0 0 −9 v −1 1 h9 The Virtual Work Principle permits identifying the dynamic entities that are duals of the kinematic entities introduced by the kinematical assumptions. Such dynamic entities are resultants of the stresses which make work (active stress). The Virtual Work Principle speci!ed for the problem under examination becomes: L L r · Duˆ d z = g · Huˆ d z ∀uˆ in [0; L] (6) 0
0
(uˆ is the admissible displacement !eld) where the resultant of active stress and external forces are collected in the two vectors rT = [N1
N2
g T = [gz1
gz2
M12
fs ];
(7)
gy
mx ]
(8)
with (see [9] for more details) N ( = 1; 2) the axial forces on A , M12 the summation of the bending moments M on A , fs the interface shear force, gz1 , gz2 , gy , mx the two axial distributed loads, the transverse distributed load and the distributed couple, respectively (for the sake of brevity
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
no forces are prescribed at the end sections), and the di2erential operator 1 0 0 0 1 0 H= 0 0 1 0
0
1911
(9)
9
is de!ned. 2.3. Constitutive law The linear generalised constitutive relation between generalised strain and stress can be written in the form: r = D”
(10)
after introducing the symmetric operator: 0 (ES)1 0 (EA)1 (EA)2 (ES)2 0 0 ; D= (ES) (EJ ) 0 (ES) 1 2 12 0
0
0
(11)
k
where the following quantities, containing the inertial properties of the cross-section and constitutive parameters, are introduced: E dA ; (12a) (EA) = A (ES) = E (y − y ) dA ; (12b) A (EJ ) = E (y − y)2 dA : (12c) A
E is the elastic modulus of the material of the -th part, (EJ )12 = (EJ )1 + (EJ )2 , k is the connection sti2ness. Since the reference points with ordinates y are assumed to coincide with the centroids of the relevant components, the static modules (12b) are zero and the matrix D is diagonal. 3. Eectiveness of the displacement-based nite element approach 3.1. Displacement-based 8nite element formulation The problem can be approached by assuming the displacements as unknown and obtaining the equilibrium conditions in terms of displacements. The variational condition (6) coincides with the stationary condition of the total potential energy functional. It is natural to seek the solution in the
1912
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
subspace of U = W1 × W2 × V = H 1 × H 1 × H 2 consisting of functions satisfying the kinematical boundary conditions. An approximated solution can be calculated using the displacement approach of the !nite element method. The displacement formulation of the !nite element method introduces a polynomial approximation of the displacement !eld at the interior of each element: u˜ = Nd ;
(13)
where N is the matrix of the shape functions and d is the vector of the nodal displacements parameters (hereinafter •˜ denotes the approximation of !eld •). The approximated strain vector ”˜u and stress vector are obtained from the approximated displacement !eld: ”˜u = DNd = Bd;
(14)
r˜u = D”˜u = DBd:
(15)
The projection of the variational condition (6) in the !nite dimensional space U˜ de!ned by d leads to the following algebraic problem: KPE d − f = 0; where
KPE = f=
0
Le
B T DB d z;
(17a)
(HN )T g d z
(17b)
0 Le
(16)
are the element sti2ness matrix and the element loads vector, respectively. 3.2. Slip locking in the 8DOF displacement-based element Di2erent choices can be made to approximate the displacement !eld in the interior of each element [3,9]. The simplest model providing the lowest regularity required involves an 8DOF !nite element where the nodal parameters are the following: d8TDOF = [w1L w2L vL ’L w1R w2R vR ’R ]
(18)
(the apex L denotes the left node and R denotes the right node of the !nite element). In this case, de=ections are described by a third-order polynomial ensuring the continuity of the !rst derivative and the axial displacements of the components are described by linear functions. The matrix N of the shape functions becomes $1 0 0 0 $2 0 0 0 ; 0 0 0 $ 0 0 0 $ N = (19) 1 2 0 %3 %4 0 0 %1 %2 0
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930 2000 300
1400 mm 2
250 1400 mm 2
22
Slab concrete Ec = 30500 MPa Slab rebars steel Es = 200000 MPa
14
1600
1913
50
Beam steel Es = 200000 MPa
500
Fig. 2. Cross-section and materials properties of tested beams.
where: 1 (1 − &); (20a) 2 1 $2 (z) = (1 + &); (20b) 2 1 %1 (z) = (2 − 3& + &3 ); (20c) 4 Le (20d) %2 (z) = (1 − & − &2 + &3 ); 8 1 (20e) %3 (z) = (2 + 3& − &3 ); 4 Le (−1 − & + &2 − &3 ) %4 (z) = (20f) 8 (& = 2z=Le − 1, & ∈ [ − 1; 1], z ∈ [0; Le ]) are the shape functions and Le is the element length. The sti2ness matrix is reported in the appendix. The numerical results show that the solution accuracy obtained by the 8DOF element strongly depends on the connection sti2ness and the error grows when the connection sti2ness becomes high. As an example, a simply supported beam is considered (L = 30 m, see Fig. 2 for details of the cross-section) and the diagrams of the curvature and slip are reported in Fig. 3 for two values of the shear connection sti2ness, measured by the non-dimensional parameter L [19], where 1 h2 1 2 =k (21) + + (EA)1 (EA)2 (EJ )12 $1 (z) =
and L is equal to the beam length in the case of a simply supported static scheme. From the numerical results, it is evident that for low-sti2ness shear connection ( L = 1), the 8DOF element is able to approximate both the exact curvature and slip, even when just a few elements are used. However, for high connection sti2ness ( L = 10), the curvature computed by the 8DOF element is forced to be constant along the element and an oscillatory trend arises in the slip description. In this regards, it can be observed that the variational Eq. (6) deriving from the displacement approach may be posed in the following form: L L L D0 D0 u · D0 uˆ d z + k D1 u · D1 uˆ d z = g · Huˆ d z ∀u; ˆ (22) 0
0
0
1914
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L 0.0
0.0
0.2
0.2
0.4
exact
0.8
2+2 8DOF
αL = 1
1+1 8DOF
0.6 0.8
2+2 8DOF
1.0 1.2
4+4 8DOF
0.4
0.6
χ/χmax
χ/χmax
L
1.0
1+1 8DOF
1.4
1.2
1.5
1.5
αL = 10
exact
2+2 8DOF exact 1+1 8DOF 2+2 8DOF
0.5
δ/δ max
1.0
0.0 -0.5
1+1 8DOF
0.5
δ/δ max
1.0
4+4 8DOF
0.0
-0.5
-1.0
-1.0
αL = 1
-1.5
αL = 10
exact
-1.5
Fig. 3. Locking problems of the 8DOF element.
where the term depending on the sti2ness parameter k has been separated, the di2erential operator has been split in 9 0 0 0 (23a) D0 = ; 0 9 0 0 −92 D1 = [ − 1
1
h9]
(23b)
and D0 is constituted by the !rst three rows and columns of D. In the limit case where k → ∞, the solution must lie in the subspace U0 of the functions satisfying the condition [11]: L
0
D1 u · D1 u d z = 0:
(24)
This constraint induces a linear relation between the displacements !elds ensuring that the slip is null, i.e. lim = lim (w2 − w1 + hv ) = 0: (25) k →∞
k →∞
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1915
So that one displacement !eld, as an example w1 , can be deduced from the other displacement !elds and U0 can be assumed as U0 = W2 × V . In this case, as in other multi!eld problems [12], troubles arise in applying the constraint in the !nite dimensional approximating space. In fact, the !nite dimensional spaces U˜ = W˜ 1 × W˜ 2 × V˜ generated by the polynomial sequences reduce, in the limit case k → ∞, to U˜ 0 = W˜ 2 × V˜ 0 where V˜ 0 ⊂ V˜ . In other words, the constraint does not simply induce a linear relation between the unknowns but also reduces the dimension of the space V˜ in which the transversal displacements are sought. For the sake of brevity, proof of what asserted is given for the particular case of the 8DOF element; the demonstration, however, can be generalized in a simple way in order to consider other !nite elements. In the 8DOF element, the slip locking problem derives from the di2erent dimensions of the spaces involved in the expression of the coupling term , where the term provided by the axial displacements lies in the two-dimensional space of !rst-order polynomials, while the term furnished by the transverse displacement lies in the three-dimensional space of the second-order polynomials. In fact, the substitution of the 8DOF displacement !eld, whose shape functions are grouped in matrix (19), in the slip expression (4) leads to (z) = a
1+& 1−& 3h c(1 − &2 ) +b + 2 2 2Le
(26)
or the equivalent form, obtained from the previous equation by collecting the terms of the same degree: h 3h 2 1 1 a + b + 3 c + (a − b)& − (z) = c& ; (27) 2 Le 2 2Le where a = −w1R + w2R + h’R ;
(28a)
b = −w1L + w2L + h’L ;
(28b)
Le (’L + ’R ): (28c) 2 Condition (25) leads to three equations ensuring that the second-order polynomial describing the slip is null. The !rst two (i.e. a = 0 and b = 0) involve both parameters describing axial and transverse displacement and express the kinematical restraint provided by a rigid shear connection: c = v R − vL −
w1L − w2L = h’L ;
(29a)
w1R − w2R = h’R ;
(29b)
while the third equation (i.e. condition on the quadratic term c = 0) involves only parameters related to transverse displacement: Le (’R + ’L ) = vR − vL ; 2
(29c)
1916
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
such condition does not re=ect the physical relation arising from axial and transverse displacement but introduces a further constraint on the de=ection DOF. The substitution of this last condition (29c) in the expression of the 8DOF de=ection v(z) leads to the following reduced form: v(z) = vL + ’L Le 18 (3 + 2& − &2 ) + ’R Le 18 (1 + 2& + &2 );
(30)
where the third-order terms, proportional to the parameter c (28c), disappear and v(z) is described by a second-order polynomial only; in other words, the e2ective DOF are reduced. Consequently, the curvature is constrained to be a constant along the element. In the limit condition k → ∞, this a2ects the slip description producing a spurious oscillatory trend in the solution [12]. The situation described occurs every time the contribution in slip at the interface deriving from axial displacements and from de=ection are described by means of polynomials of di2erent order (i.e. w and v are polynomials of di2erent order). This leads to one or more conditions involving only axial or transverse displacement !elds that cannot express the kinematical restraint provided by a rigid shear connection but introduces further constraint on the axial or on the de=ection !elds, respectively. Obviously, the poorer the element, the more evident the consequent locking will be. 3.3. Re8ned displacement elements By starting from the 8DOF element, the problem can be overcome by introducing additional internal nodes [3,9]. If an internal with a DOF is introduced for each of the two axial displacements, the 10DOF element is obtained [3,9]. In this element, both contributions to slip description, deriving from axial and transverse approximations, are second-order polynomials. Other elements based on the choice of calibrated internal nodes can be de!ned. If a further node, with one transverse and one rotational DOF, is used in the transverse !eld description (v is now a !fth-order polynomial), then three internal nodes are needed for each of the two axial displacement !elds. The resulting element (16DOF) describes the slip with a fourth-order polynomial [3,9]. 3.4. Comparisons between displacement elements This section reports some numerical results in order to compare the behaviour of the approximate solution with respect to the connection sti2ness and the convergence rate performed by the three previous !nite elements, 8DOF, 10DOF and 16DOF. Firstly, a direct comparison between the solution obtained with the 8DOF, 10DOF and 16DOF !nite elements is reported in Fig. 4, where the same problem previously analysed in Fig. 2 is considered and a two-element discretization is adopted. It is evident that the accuracy drastically improves for high connection sti2ness by using the 10DOF or the 16DOF !nite elements. In fact, elements with internal nodes furnish almost the same accuracy both for low and high sti2ness levels. For a more detailed analysis of the behaviour of di2erent elements, a cantilever with length L=2 is considered. In this case, the exact (exponential) solution di2ers more remarkably from the solutions provided by the shape functions and the approximation errors are larger and more evident. The same cross-section of the previous applications is considered and the connection sti2ness is described by the parameter L earlier introduced.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L
1917
L
0.0
0.0
0.2
0.2
10DOF 0.6
10DOF
0.8
exact
16DOF exact
1.0
αL= 1
α L = 10 1.2
1.5
1.5
8DOF 10DOF 16DOF exact
0.5 0.0
1.0 0.5
δ/δ max
1.0
δ/δ max
8DOF 0.6 0.8
16DOF
1.0 1.2
0.4
8DOF
χ/χ max
χ/χ max
0.4
10DOF
16DOF
0.0
exact
-0.5
-0.5
-1.0
-1.0
8DOF
-1.5
αL = 1
-1.5
α L = 10
Fig. 4. Comparisons between displacement elements: curvature and slip in the simply supported beam.
Fig. 5 describes the curvature obtained for four di2erent values of the connection sti2ness parameter ( L = 1, 7, 20, 50). This comparison intends to analyse the sensitivity with respect to L and to present the distribution of the local error along the whole beam. Discretisation adopts only one element; it is evident that the 16DOF element has more DOF describing the curvature (third-order polynomial) and potentially can better !t the exact solution with respect to both the 8DOF and 10DOF elements furnishing a linear description of the curvature. It can be observed that 10DOF and 16DOF elements present solutions with the same order of accuracy for all the connection sti2ness considered. It is also interesting to observe the di2erences between the 8DOF and the 10DOF elements, for which the curvature is described by the same polynomial order. Fig. 6 reports results concerning the interface slip. In this case, the shape of the solution change very considerably by varying from low to high connection sti2ness values. The 8DOF and 10DOF elements can provide a parabolic description while the 16DOF element can describe the interface slip by means of a fourth-order polynomial. The 16DOF element furnishes satisfactory results at all the connection sti2ness levels, while the 8DOF element local error is larger than the local error deriving from the 10DOF element.
1918
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2
L/2 1.2
1.2
16DOF
1.0
0.8
10DOF
0.6
8DOF
0.4 0.2
χ/χ max
χ/χ max
0.8
-0.2
αL = 1
-0.4
-0.4
1.0
1.0
16DOF
αL = 7
16DOF
0.8
10DOF
0.6
8DOF
0.4
χ/χ max
0.8
χ/χ max
exact
1.2
1.2
10DOF
0.6
8DOF
0.4 0.2
0.2
exact
0.0
-0.4
8DOF
0.4
0.0
0.0
-0.2
10DOF 0.6
0.2
exact
-0.2
16DOF
1.0
αL = 20
0.0 -0.2 -0.4
exact
αL = 50
Fig. 5. Comparisons between displacement elements: curvature in the cantilever.
The global error is analysed in Figs. 7–9 and in Tables 1–4. The error is evaluated both for the simply supported static scheme and the cantilever. Fig. 7 and Tables 1 and 2 report the energy error e=
u − u˜h ;
u
(31)
where
u =
0
L
DDu · Du d z:
(32)
Fig. 8 and Tables 3 and 4 describe the results obtained for the maximum displacement (mid-span for the simply supported beam and end displacement for the cantilever). The results concerning the 16DOF element are not showed for the simply supported beam because they are out of scale (the error is too small to be depicted). The error drastically reduces by adopting elements with internal nodes (diagrams are traced in a logarithmic scale). It can be observed that the error always increases when the connection sti2ness increases only in the case of the 8DOF element.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2
L/2 1.2
1.2
αL = 2
0.8
δ/δ max
0.6 0.4
8DOF
0.8
exact 16DOF
exact 0.6
16DOF 0.4
8DOF
0.2
αL = 7
1.0
1.0
δ/δ max
1919
0.2
10DOF
10DOF
0.0
0.0
1.5
1.5
8DOF
8DOF
1.0
1.0
exact 0.5
10DOF
0.0
δ/δmax
δ/δmax
0.5
16DOF
10DOF
-0.5
-0.5
16DOF
-1.0
-1.0 -1.5
exact 0.0
α L = 20
α L = 50
-1.5
Fig. 6. Comparisons between displacement elements: slip in the cantilever. 100.000%
100.000% 1 8DOF
1+1 8DOF
2 8DOF
10.000%
2+2 8DOF
energy error
energy error
10.000% 1.000%
1+1 10DOF 4+4 8DOF
0.100% 2+2 10DOF
0.010%
1 10DOF
1.000%
4 8DOF 2 10DOF
0.100%
1 16DOF
0.010%
4 10DOF 2 16DOF
4+4 10DOF
0.001%
1
10
αL
100
0.001%
1
10
100
αL
Fig. 7. Comparisons between displacement elements: energy error versus connection sti2ness.
Finally, Fig. 9 describes the convergence rate by reporting the energy error versus the total DOF (the internal nodes DOF are computed in the total DOF). As in the previous case, the results for the 16DOF element are not showed for the simply supported beam because they are too small. The di2erent performances of the three elements are, once again, evident.
1920
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930 100.000%
100.000%
1 8DOF
1+1 8DOF
v
10.000%
2+2 8DOF
nodal error
nodal error
10.000%
4+4 8DOF
1.000% 0.100%
4 8DOF
1.000%
1 10DOF
0.100%
2 10DOF
1+1 10DOF
0.010%
2 8DOF
v
4 10DOF
0.010%
1 16DOF
2+2 10DOF
0.001%
2 16DOF
0.001% 1
10
100
1
10
αL
100
αL
Fig. 8. Comparisons between displacement elements: nodal error versus connection sti2ness.
100.000%
100.000% 8DOF αL = 20
8DOF αL = 20
10.000%
8DOF αL = 7
energy error
energy error
10.000% 1.000%
10DOF αL = 7
0.100%
8DOF αL = 7
1.000%
10DOF αL = 7 10DOF αL = 20
0.100%
16DOF αL = 20
0.010%
0.010% 10DOF αL = 20
total DOF
0.001% 0
5
10
15
20
25
0
30
16DOF αL = 7
total DOF
0.001%
5
10
15
20
25
30
Fig. 9. Comparisons between displacement elements: energy error versus total DOF.
Table 1 Energy error in the simply supported beam under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1+1 el. 2+2 el. 4+4 el.
1.06 0.07 0.01
2.24 0.43 0.10
4.68 1.09 0.27
7.43 1.79 0.45
10.28 2.50 0.62
13.30 3.26 0.79
14.81 3.71 0.88
16.42 4.52 1.04
16.81 4.93 1.21
10DOF
1+1 el. 2+2 el. 4+4 el.
1.03 0.07 0.00
0.88 0.06 0.00
0.86 0.06 0.00
0.89 0.06 0.01
0.93 0.06 0.01
0.99 0.07 0.01
1.03 0.07 0.01
1.06 0.07 0.01
1.07 0.07 0.01
16DOF
1+1 el.
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1921
Table 2 Energy error in the cantilever under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1 el. 2 el. 4 el.
3.13 0.21 0.02
6.24 0.65 0.11
10.25 1.36 0.27
14.46 2.23 0.47
19.52 3.45 0.74
27.38 5.70 1.20
33.26 7.68 1.59
41.61 11.44 2.48
43.99 17.31 8.75
10DOF
1 el. 2 el. 4 el.
3.10 0.20 0.02
4.85 0.33 0.03
6.07 0.45 0.03
6.55 0.55 0.04
6.31 0.63 0.06
5.31 0.69 0.07
4.51 0.65 0.09
3.36 0.41 0.08
2.96 0.25 0.04
16DOF
1 el. 2 el. 4 el.
0.00 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.00
0.01 0.00 0.00
0.04 0.00 0.00
0.09 0.00 0.00
0.13 0.01 0.00
0.11 0.02 0.00
0.04 0.01 0.00
Table 3 Nodal de=ection error in the simply supported beam under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1+1 el. 2+2 el. 4+4 el.
0.02 0.01 0.00
1.19 0.35 0.09
3.80 1.00 0.25
6.90 1.69 0.42
10.35 2.38 0.59
14.39 3.10 0.74
16.63 3.53 0.82
19.18 4.37 0.95
19.83 4.83 1.12
10DOF
1+1 el. 2+2 el. 4+4 el.
0.05 0.00 0.00
0.22 0.01 0.00
0.25 0.01 0.00
0.21 0.01 0.00
0.13 0.00 0.00
0.05 0.00 0.00
0.02 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.00
16DOF
1+1 el.
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
Table 4 Nodal de=ection error in the cantilever under uniform load L (%) 1
3
5
7
10
15
20
50
100
8DOF
1 el. 2 el. 4 el.
0.11 0.01 0.00
2.19 0.40 0.09
5.63 1.07 0.25
9.61 1.84 0.43
14.38 2.74 0.63
21.07 4.06 0.89
25.52 5.09 1.08
31.34 7.10 1.49
32.91 8.02 1.85
10DOF
1 el. 2 el. 4 el.
0.08 0.01 0.00
0.72 0.05 0.00
1.25 0.10 0.01
1.52 0.15 0.01
1.49 0.19 0.02
1.09 0.21 0.02
0.75 0.19 0.03
0.24 0.09 0.02
0.06 0.03 0.01
16DOF
1 el. 2 el.
0.00 0.00
0.00 0.00
0.00 0.00
0.01 0.00
0.02 0.00
0.03 0.00
0.05 0.00
0.04 0.01
0.02 0.00
1922
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
4. Eectiveness of nite element mixed approaches 4.1. Three-8eld 8nite element mixed formulation The problem can be approached by assuming the displacement, strain and stress !elds as unknowns and solved using the compatibility condition, balance condition and the generalised constitutive law. An approximated solution can be obtained using a three-!eld formulation of the !nite element method. The method introduces polynomial approximations of the displacement !eld, of the strain !eld and of the stress !eld at the interior of each element: u˜ = Nd ;
”˜ = Ee; r˜ = Ss;
(33a,b,c)
where N , E and S are the matrices of shape functions and d , e and s are the vectors of displacement, strain and stress nodal parameters. The projection of the weak form of the balance, compatibility and constitutive equations (details in [9]) in the !nite dimensional space de!ned by d , e, s leads to the following solving system: BV T s − f = 0;
(34a)
V − EV T s = 0; De
(34b)
V − Bd V = 0; Ee
(34c)
where BV = EV = DV =
Le
0
Le
0
0
Le
S T B d z;
(35)
S T E d z;
(36)
E T DE d z:
(37)
Since the stress and strain !elds can be inter elements discontinuous, the system of equations (34) can be simpli!ed performing a static condensation of the stress and strain nodal parameters [9]. The resulting system has the same shape of the solving system derived in displacement approach (16) and can be easily used in a displacement-based !nite element program. 4.2. Reduction of slip locking in mixed elements The authors have already presented and tested three-!eld mixed elements for composite beams with deformable shear connection [9] based on the locking-free displacement !eld of the 10DOF
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
1923
Table 5 Shape functions polynomials degree Element
8DOF HW011 HW010
Displacements
Strains
Stresses
w
v
N
M
fs
1 1 1
3 3 3
0 0 0
1 1 1
2 1 0
0 0 0
1 1 1
2 1 0
displacement element. In this paper, the e2ects of slip locking in mixed elements will be illustrated by introducing mixed elements based on the 8DOF displacement !eld. Once the displacement !eld is selected, the Limitation Principles [20,21] must be considered for a proper de!nition of the strain and stress approximations. Roughly speaking, the Limitations Principles state that the displacement approach and the three-!eld mixed approach become identical if the strain approximation ”˜ and the stress approximation r˜ are equal or richer than the strain ”˜u and the stress r˜u calculated from the approximation of the displacement !eld. Di2erent mixed elements were tested by the authors using the same shape functions for conjugated strains and stresses (i.e. the matrices E and S are equal). The tests showed that it is necessary to reduce approximation of the slip (and the interface force) in order to reduce locking. As a result, four choices are possible: constant or linear slip (interface force fs ) combined with constant axial strains (axial forces N ) and constant or linear curvature (moment M12 ). However, elements with constant curvature have a too poor behaviour if compared to the element with linear curvature, consequently only the results for elements with linear curvature are discussed in this paper. In Table 5, a comparisons between the shape functions of the 8DOF element and of the three-!eld mixed element with linear curvature and reduced slip (HW010 and HW011) is illustrated. The strain and stress polynomial degree of the displacement-based element are indicated in italics because they are not independent from the displacement shape functions but derived from them (using compatibility and linear constitutive law). It is worth saying that if the displacement–stress two-!eld mixed formulation is adopted and if the same representations is used for the stress both in the three-!eld and two-!eld formulations, then the same results are obtained as a consequence of the Limitation Theorems (two-!eld and three-!eld formulations coincide in the case considered). Similar considerations apply to the displacement– strain two-!eld mixed formulation, previously introduced by the authors [22] before the more general three-!eld formulation. 4.3. Comparisons between mixed and displacement elements This section reports some numerical results in order to compare the behaviour of the approximate solution, furnished by the !nite elements 8DOF, HW010 and HW011, with respect to the variation of the connection sti2ness. As in the previous case, a cantilever with length L=2 is considered. Fig. 10 describes the curvature obtained for four di2erent values of the connection sti2ness parameter ( L = 1, 7, 20, 50). For low sti2ness values ( L = 1 and L = 7), the three elements give
1924
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2
L/2 1.2
1.2 1.0
exact
1.0
HW011 HW010
0.4
χ/χmax
χ/χmax
8DOF
0.6
0.2
8DOF
0.0
0.0 -0.2
1.0
αL= 1
-0.4
0.8
exact
8DOF HW011
χ/χ max
0.6
0.4
8DOF
0.4 0.2
HW011
0.0
0.0
-0.4
αL = 7
1.0
exact
0.6
-0.2
HW010
1.2
0.8
0.2
HW011
0.4 0.2
1.2
χ/χmax
0.6
-0.2 -0.4
exact
0.8
0.8
-0.2
αL = 20
HW010
-0.4 -0.6
HW010 αL = 50
Fig. 10. Comparisons between mixed and displacement elements: curvature in the cantilever.
the same results (both displacement element and mixed elements describe the curvature with a linear function). When sti2ness increases, it is evident that the curvatures computed by the 8DOF and by the HW011 are forced to be constant along the beam axis, while the HW010 linear curvature does not degenerate into a constant. Fig. 11 describes the interface slip obtained for the same four di2erent values of the connection sti2ness parameter ( L = 1, 7, 20, 50). For low sti2ness values, the 8DOF element is locally more accurate than the mixed elements, but when sti2ness increases the HW010 element furnishes lower local errors. Notice that the mixed elements cannot satisfy the boundary conditions. Finally, in Figs. 12 and 13, the same problem is analysed using a two-element discretisation: it can be observed that a more re!ned mesh reduces the local errors but do not eliminate the described problems. In fact, the curvature trends for the 8DOF and HW011 elements are forced to be a constant along each element, while slip oscillations once again arise. The convergence curves for the mixed elements, for the sake of brevity not reported here, are very close to the curves of the 8DOF element.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2 1.4
L/2 1.4
αL = 1
1.2
1.0
HW010
δ/δ max
δ/δmax
0.8
exact
αL = 7
8DOF
1.2
8DOF
1.0
0.6
1925
HW010
0.8
HW011
0.6
exact
HW011 0.4
0.4
0.2
0.2
0.0
0.0
2.5
1.5
HW011
HW011
2.0
HW010
0.5
1.5
δ/δmax
δ/δmax
1.0
exact 0.0
HW010
1.0 0.5 0.0
exact
-0.5
-0.5 -1.0
8DOF αL = 20
8DOF
-1.0 -1.5
αL = 50
Fig. 11. Comparisons between mixed and displacement elements: slip in the cantilever.
5. Conclusions In this paper, the authors discuss problems that may take place in the !nite element analysis of composite beams with deformable shear connection, where slip occurs between beam components while contact is preserved. A coupling between the two axial displacement !elds and the =exural displacement !eld arise in this beam model due to the deformable interface. If the axial and =exural displacement shape functions are not consistent, then locking problems occur for high values of the connection sti2ness (slip locking). Slip locking is comparable to membrane and shear locking that may develop in common beam elements, and it mostly a2ect the description of the curvature and interface slip. This is a problem of particular interest in the non-linear analysis of composite beams, because the non-linear behaviour of the shear connection is characterised by a constitutive law with a sti2ness varying from very a high value when the slip is close to zero, to a very low value when the connection is close to failure. Since the distribution of interface slip varies greatly along the beam, the solution obtained, which depends strongly on the connection sti2ness, is not accurate. A strategy to avoid locking by means of a calibrated choice of the
1926
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
αL/2
L/2
1.2
1.2
exact
1.0
1.0
0.8
exact
0.8
8DOF
HW011 HW010
0.4
χ/χmax
χ/χmax
8DOF 0.6
0.2
HW010
0.0
αL = 1
-0.2
1.2 1.0
HW011
0.4 0.2
0.0 -0.2
0.6
αL = 7
1.2
exact
1.0
0.8
exact
0.8
8DOF
HW011
0.4
HW010
0.2
0.6
HW011 HW010
0.4 0.2
0.0 -0.2
χ/χmax
χ/χmax
8DOF 0.6
0.0
αL = 20
-0.2
αL = 50
Fig. 12. Comparisons between mixed and displacement elements: curvature in the cantilever (two-element mesh).
displacement shape functions is then proposed and a comparison between di2erent displacement-based elements, in terms of local error, global error and convergence rate, is illustrated. The authors also analyse the e2ect of locking on mixed elements. Starting from the displacement !eld of the element with shape functions of the lowest degree required by the problem (8DOF element), three-!eld mixed elements are introduced. Since the 8DOF displacement element is a2ected by locking, the applications intend to clarify if the mixed elements are still a2ected by slip locking. The numerical tests performed show that a mixed element based on a locking-su2ering displacement element can still be a2ected by locking. It is possible to reduce locking by lowering the degree of the shape functions for the slip and interface force, but this leads to poor description of these entities. Locking-free displacement elements, formulated accordingly to the indications of consistence between axial and =exural !elds, and mixed elements based on locking-free displacement shape functions, allow to obtain more accurate and reliable results in both linear and non-linear analyses.
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L/2 1.2
L/2 1.4
αL = 1
8DOF 1.0
exact
δ/δmax
δ/δmax
αL = 7
HW011
1.2
1.0 0.8 0.6
8DOF
0.8
HW010
0.6
exact
HW010
0.4
1927
0.4
HW011
0.2
0.2 0.0
0.0
3.0
2.0
HW011 2.5
1.5
HW010 exact
δ/δmax
δ/δmax
1.0 0.5
1.5
HW010
1.0 0.5
0.0
0.0 -0.5 -1.0
HW011
2.0
8DOF
-0.5
αL = 20
-1.0
exact 8DOF
αL = 50
Fig. 13. Comparisons between mixed and displacement elements: slip in the cantilever (two-element mesh).
Appendix A A.1. The 8DOF element sti>ness matrix
K8DOF = where
KLL
KLR
T KLR
KRR
EA1 L 0 KLL = 0 0
;
(A.1)
0
0
EA2 L
0
0 0
EJ12 L3 EJ12 6 2 L
12
L 0 3 L − 0 3 +k h EJ12 6 2 2 L hL EJ12 4 − L 12
L − 3 L 3 h − 2 hL 12
h 2 h − 2 6h2 5L h2 10
hL − 12 hL 12 ; 2 h 10 2h2 L 15
(A.2)
1928
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
L 0 0 0 3 L EA2 − 0 0 3 L +k h EJ12 EJ12 − 0 12 3 −6 2 2 L L hL EJ12 EJ12 0 −6 2 4 − L L 12 L EA1 − 0 0 0 L 6 L EA − 2 0 0 0 − 6 L +k KLR = h EJ12 EJ12 0 0 −12 6 3 2 2 L L hL EJ12 EJ12 0 0 −6 2 2 L L 12 EA1 L 0 KRR = 0 0
L − 3 L 3 h 2 hL 12
L 6 L 6 h − 2 hL − 12 −
hL − 12 hL 12 ; 2 h − 10 2h2 L 15 h hL − 2 12 h −hL 2 12 : 6h2 h2 − 5L 10 h2 h2 L − − 10 30
h − 2 h 2 6h2 5L h2 − 10
(A.3)
(A.4)
A.2. The HW010 mixed element sti>ness matrix
KV HW010 = where
KLL
KLR
T KLR
KRR
EA1 L 0 KLL = 0 0 EA1 L 0 KRR = 0 0
;
(A.5)
0
0
EA2 L
0
0 0
EJ12 L3 EJ12 6 2 L
12
0
0
EA2 L
0
0 0
EJ12 12 3 L EJ12 −6 2 L
L h L 0 − 4 4 2 L h L 0 − − 0 +k 4 4 2 ; EJ12 2 h h h 6 2 − 0 L 2 2 L EJ12 4 0 0 0 0 L L L h 0 − − 0 4 4 2 L L h 0 − 0 +k 4 4 2 ; EJ12 2 h h h −6 2 0 − L 2 2 L EJ12 4 0 0 0 0 L 0
(A.6)
(A.7)
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
KLR =
−
EA1 L 0
−
0
0
0
EA2 L
0
0
0
0
0
0
EJ12 L3 EJ12 −6 2 L
−12
EJ12 L2 EJ12 2 L 6
L 4 L − +k 4 h 2 0
L 4 L 4 h − 2 0
−
h 2 h 2 h2 − L 0 −
0 0 : 0 0
1929
(A.8)
References [1] D.J. Oehlers, M.A. Bradford, Elementary Behaviour of Composite Steel and Concrete Structural Members, Butterworth-Heinemann, London, 2000. [2] N.M. Newmark, C.P. Siess, I.M. Viest, Tests and analysis of composite beams with incomplete interaction, Proc. Soc. Exp. Stress Anal. 9 (1) (1951) 75–92. [3] A. Dall’Asta, A. Zona, Non-linear analysis of composite beams by a displacement approach, Comput. Struct. 80 (27–30) (2002) 2217–2228. [4] M.R. Salari, E. Spacone, P.B. Shing, D.M. Frangopol, Non-linear analysis of composite beams with deformable shear connectors, J. Struct. Eng. ASCE 124 (10) (1998) 1148–1158. [5] M.R. Salari, E. Spacone, Finite element formulations of one-dimensional elements with bond-slip, Eng. Struct. 23 (7) (2001) 815–826. [6] A. Ayoub, F.C. Filippou, Mixed formulation of nonlinear steel–concrete composite beam element, J. Struct. Eng. ASCE 126 (3) (2000) 371–381. [7] A. Ayoub, A two-!eld mixed variational principle for partially connected composite beams, Finite Elem. Anal. Des. 37 (2001) 929–959. [8] B. Cas, F. Saje, M. Saje, I. Planinc, Nonlinear !nite element analysis of composite planar frames with interlayer slip, in: B.H.V. Topping, Z. Bittnar (Eds.), Proceedings of the Sixth International Conference on Computational Structures Technology, Civil-Comp Press, Stirling, Scotland, 2000, Paper 135. [9] A. Dall’Asta, A. Zona, Three-!eld mixed formulation for the non-linear analysis of composite beams with deformable shear connection, Finite Elem. Anal. Des. 40 (4) (2004) 425– 448. [10] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, New York, 2000. [11] D. Chapelle, Reliability of !nite element methods for thin shells, in: B.H.V. Topping (Ed.), Computational Mechanics for the Twenty-First Century, Saxe-Coburg Publications, Edinburgh, 2000, pp. 99–108. [12] G. Prathap, B.P. Naganarayana, Stress oscillations and spurious load mechanisms in variationally inconsistent assumed strain formulations, Int. J. Numer. Methods Eng. 33 (1992) 2181–2197. [13] M.A. Cris!eld, Non-linear Finite Element Analysis of Solids and Structures, Vol. 1: Essentials, Wiley, New York, 1991. [14] L. Yunhua, Explanation and elimination of shear locking and membrane locking with !eld consistency approach, Comput. Methods Appl. Mech. Eng. 162 (1998) 249–269. [15] J.N. Reddy, On locking-free shear deformable beam !nite elements, Comput. Methods Appl. Mech. Eng. 149 (1997) 113–132. [16] I.J. Blaauwendraad, Realistic analysis of reinforced concrete framed structures, Heron 18 (4) (1972) 1–31. [17] A.K. Gupta, P.S. Ma, Error in eccentric beam formulation, Int. J. Numer. Methods Eng. 11 (1977) 1473–1483. [18] M.A. Cris!eld, The eccentricity issue in the design of plate and shell elements, Commun. Appl. Numer. Methods 7 (1991) 47–56. [19] U.A. Girhammar, D. Pan, Dynamic analysis of composite members with interlayer slip, Int. J. Solids Struct. 30 (6) (1993) 797–823.
1930
A. Dall’Asta, A. Zona / Finite Elements in Analysis and Design 40 (2004) 1907 – 1930
[20] B.M. Fraeijs De Veubeke, Displacement and equilibrium models in the !nite element method, in: O.C. Zienkiewicz, G.S. Hollister (Eds.), Stress Analysis, Wiley & Sons, London, 1965, pp. 145–196 Reprinted in: Int. J. Numer. Methods Eng. 52 (2001) 287–342. [21] H. Stolarski, T. Belytschko, Limitation Principles for mixed !nite elements based on the Hu-Washizu variational formulation, Comput. Methods Appl. Mech. Eng. 60 (1987) 195–216. [22] A. Dall’Asta, A. Zona, Finite elements for the analysis of composite beams with interlayer slip, Proceedings of the XVIII Congresso C.T.A., Vol. 2, Venezia, Italy, 26–28 September 2001, pp. 365–374.