- Email: [email protected]
Thin plate element for applied element method
Structures 22 (2019) 1–12
Contents lists available at ScienceDirect
Structures journal homepage: www.elsevier.com/locate/structures
Thin plate element for applied element method D. Lincy Christy , T.M. Madhavan Pillai, Praveen Nagarajan ⁎
T
Department of Civil Engineering, National Institute of Technology Calicut, India
ARTICLE INFO
ABSTRACT
Keywords: Applied element method Thin plate element Stiffness matrix Poisson effect
Applied Element Method (AEM) is an efficient method for numerical analysis of structures. It is capable of analyzing the structure up to collapse. In AEM, the structure is modelled as a group of rigid elements connected by springs. It has the benefit of lesser processing time and memory requirement when compared to Finite Element Method. It is also able to track crack propagation more realistically. To extend the advantages of AEM to the analysis of thin plate structures, thin plate element has to be developed. This is done in this paper. Being a rigid body method, AEM does not consider Poisson's ratio in its basic formulation. As Poisson's ratio plays a major role in the behavior of plate structures, the method to bring in Poisson effect in thin plate element is discussed. The procedure to determine displacements, strains, stresses and stress resultants are explained. The mass matrix for thin plate element is also found. To validate the thin plate element, both static and dynamic analyses are carried out. The analyses gave reliable results which show the efficiency of the thin plate element of AEM.
1. Introduction Applied Element Method (AEM) is an outcome of discrete element method used for structural analysis. Meguro and Tagel-Din introduced this numerical method in 1997 [1]. They advantageously used this method for the fracture analysis of reinforced concrete structures [1,2]. They detailed the procedure to determine crack pattern using AEM [3]. The procedure to consider Poisson effect in two-dimensional (2D) element was also discussed by them [4,5]. With the help of this, they analysed elastic prism subjected to in-plane loading. In 2000, a method to find mass matrix for 2D element in AEM was devised [6]. The geometric nonlinearity was handled by residual force vector generated due to geometric changes instead of geometric stiffness matrix [6,7]. Khalil analysed moment-resisting steel frame using AEM for progressive collapse assessment [8]. Cismasiu et al. studied the failure of reinforced shear wall using AEM in 2017 [9]. In 2018, the study on the behavior of brick masonry wall under in-plane cyclic loading [10] and the convergence study on plain concrete beams [11] were done. As springs are used to model the structure, the material property can be easily varied within the structure in AEM. This is advantageous for the analysis of heterogeneous materials such as functionally graded materials. In 2015, Shakeri and Bargi did static analysis of rectangular column, beam and slab to study the effect of number of springs and elements on the accuracy of results [12]. The Poisson's ratio was taken as 0 for all the cases. The analysis of rectangular slab with transverse loading was
⁎
accomplished using three-dimensional (3D) elements. They plotted the normal stress and shear stress distribution in rectangular beam. Dynamic analysis was also carried out with the help of Newmark-β method. From the literature it is seen that, analysis of thin plate with transverse loading and considering Poisson's ratio has not been done. In this paper, a thin plate element is developed for AEM. The stiffness matrix for the analysis of thin plate is derived from the stiffness matrix of a pair of springs in 3D analysis. The effect of Poisson's ratio on the behavior of plate structures is significant. So, the technique to consider the effect of Poisson's ratio is illustrated using the analysis with thin plate element. The method to calculate displacements, strains, stresses, bending moments, twisting moment and shear force are described. Static and dynamic analyses of plates are performed to verify the thin plate element developed. Plates of different boundary condition and loading are analysed to find the deflection, moments and shear force. The normal stress distribution and shear stress distribution are also plotted. Convergence study with respect to number of elements and number of springs is done. All the results are compared with that of Classical Plate Theory (CPT). The mass matrix of thin plate element is determined to carry out dynamic analysis. Using this, the fundamental natural frequency of plates simply supported on all edges and clamped on all edges are found. The contour plot of the mode shapes are also drawn and verified.
Corresponding author. E-mail addresses: [email protected] (D.L. Christy), [email protected] (T.M.M. Pillai), [email protected] (P. Nagarajan).
https://doi.org/10.1016/j.istruc.2019.07.010 Received 15 May 2019; Received in revised form 22 July 2019; Accepted 23 July 2019 2352-0124/ © 2019 Institution of Structural Engineers. Published by Elsevier Ltd. All rights reserved.
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 1. (a) 3D element. (b) Connection of 3D elements.
2. Applied element method
3. Thin plate element for AEM
AEM is used for the numerical analysis of structures. In this method, the structure is divided into rigid elements connected by pairs of springs. The degrees of freedom (DOFs) are located at the centre of the elements. There are 6 DOFs for a 3D element as shown in Fig. 1(a). They are 3 translations and 3 rotations. In 3D analysis, a pair of springs consists of one spring to carry normal force (A) and two springs to take shear forces (B and C) along the transverse directions (Fig. 1(b)). Every pair of springs has a stiffness matrix corresponding to the DOFs of the elements connected by them. The stiffness matrix of a pair of springs connecting 2D elements is given in [13]. The stiffness matrix for a pair of springs (Ks) located at (x, y, z) from the centroid of the first element shown in Fig. 1(b) is given by Eq. (1).
For the analysis of plates, 2D element is sufficient, if only in-plane loading is applied. Anyway, when transverse loading is applied, thin plate element for transverse loading has to be developed. To arrive at the behavior of thin plate using 3D element, 3 major modifications are to be made. They are discussed in this section.
Ks =
kn 0 0 0 knz -kny -kn 0 ks 0 -ksz 0 ks x 0 0 0 ks ksy -ksx 0 0 0 -ksz ksy ks(y2+ z2) -ksxy -ksxz 0 knz 0 -ksx -ksyx knz2+ksx2 -knyz -knz -kny ksx 0 -kszx -knzy kny2+ksx2 kny
3.1. Modification 1: value of E′ The value of E′ depends upon the constitutive relation of the problem considered. For normal spring, strain along its direction alone is important. Hence, the value of E′ can be evaluated by ignoring other normal strains. The value of E′ for different types of problems is given in
0 0 0 -knz kn y -ks 0 ksz 0 ksx 0 -ks -ksy -ksx 0 ksz -ksy -ks(y2+ z2) -ksxy -ksxz 0 ksx ksyx -knz2+ksx2 knyz -ksx 0 kszx knzy -kny2+ksx2
-kn 0 0 0 -knz kn y kn 0 0 0 knz -kny 0 -ks 0 ksz 0 -ksx 0 ks 0 -ksz 0 -ksx 0 0 -ks -ksy ksx 0 0 0 ks ks y ksx 0 0 ksz -ksy -ks(y2+ z2) ksxy ksxz 0 -ksz ksy ks(y2+ z2) ksxy ksxz -knz 0 -ksx -ksyx -knz2+ksx2 knyz knz 0 ksx ksyx knz2+ksx2 -knyz kny ksx 0 -kszx knzy -kny2+ksx2 -kny -ksx 0 kszx -knzy kny2+ksx2
The stiffnesses of the normal spring (kn) and shear springs (ks) are given by Eq. (2).
kn =
E As GAs ; ks = ;E = l l
(1)
Table 1. Here, E and ν are the modulus of elasticity and Poisson's ratio respectively. From Table 1, it is seen that the stiffness of thin plate is more than that of beam by a factor 1/(1−ν2). The additional rigidity is due to the Poisson's ratio.
(2)
where, G, l and As are the modulus of rigidity, length of the element and area of cross-section of the portion of the element treated as the spring respectively. The normal stress and normal strain along the direction of normal spring are denoted by σ and ɛ respectively. For all the springs, x = l/2. The 3D element is the skeleton of AEM. All other types of elements can be obtained by reducing or modifying the 3D element. Similarly, the stiffness matrix for any analysis can be derived from the stiffness matrix for 3D analysis. In the case of 2D analysis, one of the shear springs is not necessary. Also, the DOFs can be reduced to 3; 2 translations and 1 rotation. The size of stiffness matrix now reduces to 6 × 6 by removing the columns and rows corresponding to the unwanted DOFs. When it comes to one-dimensional (1D) analysis, only normal spring is to be considered and a single translational DOF becomes sufficient. The size of stiffness matrix becomes as small as 2 × 2.
3.2. Modification 2: degrees of freedom required for analysis The basic behavior of thin plates helps to simplify the analysis. According to Kirchhoff's assumptions, the translations along the middle plane is 0. Therefore, if one element is taken along the depth, the first 2 DOFs [u and v in Fig. 1(a)] are not necessary. Element rotation about Zaxis is given by Eq. (3). Table 1 Value of E′ for different problems. 1D/2D E
2
Plane stress/thin plate E (1
2)
Plane strain/3D E (1 ) (1 + )(1 2 )
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 2. (a) Thin plate element (b) connection of thin plate elements.
z
=
1 2
v x
u =0 y
displaced profile of the group of elements when unit deformation is applied to the central element along its DOFs (w, θx and θy) are shown in Fig. 3(c)–(e) respectively. Unit deformation along the depth direction, does not induce any force in the normal springs and hence Poisson effect does not come into picture. However, rotation about X-axis and Y-axis strains the normal springs along Y-axis and X-axis respectively. Due to unit rotation about Y-axis (θy = 1), the line AB on the right face of the central element is displaced to A′B′ as shown in Fig. 3(e). This creates compressive strain and tensile strain above and below the middle plane respectively. The direct forces generated in the direction of X-axis due to these strains are already taken into account in the stiffness matrix. In addition to the direct forces, transverse forces develop along Y-axis due to Poisson effect. These are not considered in the stiffness matrix. The normal stress along Y-axis (σy) is:-
(3)
For thin plates, as u and v are ignored, θz turns out to be 0. Therefore, the 6th DOF also can be neglected. Hence, translation along Z-axis (w), rotation about X-axis (θx) and rotation about Y-axis (θy) are the DOFs necessary for thin plate analysis. The DOFs in thin plate element are shown in Fig. 2(a). All the 3 springs are required to connect the thin plate elements as shown in Fig. 2(b). The stiffness matrix of a pair of springs now reduces to a size of 6 × 6 [Eq. (4)].
Ks =
ks
ks y
ks l/2
ks
ks y
ks l/2
ks y
ks (y 2 + z 2)
ks ly /2
ks y
ks (y 2 + z 2)
ks ly /2
ks l/2 ks
ks ly /2 ks y
ks y
ks (y 2 + z 2)
ks l/2
ks ly /2
kn z 2 + ks l2/4 ks l/2
ks l/2 ks
ks ly /2 ks y
ks ly /2
ks y
ks (y 2 + z 2)
ks ly /2
ks ly /2
kn z 2 + ks l2/4
kn z 2 + ks l2/4 ks l/2
kn z 2 + ks l2/4 ks l/2
(4)
To get the complete stiffness matrix at a face of the plate element, the stiffness matrix of individual pairs of springs should be summed up. If infinite number of springs are adopted in the section, the springs will be of infinitesimal area. Then, the sum of the stiffness matrix of all pairs of springs (K) can be obtained by integrating the stiffness matrix across the section as given by Eq. (5). K=
K=
k1 =
y
=
0 k2 0 0 k2 0
k4 0 k3 k4 0 k5
k1 0 k4 k1 0 k4
0 k2 0 0 k2 0
x
y
Gbh k5 = 2
(
y
+
x)
= z / l and
y
=0
Therefore,
k4 0 k5 k4 0 k3
Gbh Gbh (b2 + h2) bh (Eb2 + 3Gl2) k4 = k2 = k3 = l 12l 12l
2
For the shaded quarter portion in Fig. 3(e), the normal strain along X-axis (ɛx) and normal strain along Y-axis (ɛy) are:-
K s dA k1 0 k4 k1 0 k4
E 1
=
E
z l
2
1
The total compressive force (p) and the total tensile force for a quarter portion [shaded region in Fig. 3(e)] of the element will be equal and is given by Eq. (6). The lever arm (r) between total compressive force and total tensile force is given by Eq. (7). These forces in the quarter part of the plate element creates a moment (m) at the centre of the plate given by Eq. (8).
bh (Eb2 3Gl2) 12l
(5) where, b and h are the width and thickness of the plate element. This stiffness matrix [Eq. (5)] can be directly used at a section instead of summing up the individual stiffness matrices of all pairs of springs.
h /2 l /2
p=
y 0
=
3.3. Modification 3: including the effect of Poisson's ratio
Eh2 2) 16(1
(6)
2 h 3
(7)
m = pr
(8)
r=
Poisson's ratio has a major role in the behavior of plates. Hence, Poisson effect must be invariably considered in the analysis of thin plates. The methodology to bring in Poisson effect in 2D element is discussed by Meguro and Tagel-Din [4,5]. For this, they introduced additional terms into the stiffness matrix. According to the Modification 2, a thin plate element has 3 DOFs which are different from those in 2D element. Hence, a method has to be devised to include Poisson effect in thin plate element. This is discussed in the next section.
dxdz
0
In the case of rectangular plate element, the total compressive force and total tensile force (p) are given by Eq. (9). h /2 b /2
p= 0 0 h /2 l /2
E 2
1 E
z Eh2 dydz = for 2) b 16(1
x
=1
z Eh2 dxdz = for 2) l 16(1
y
=1
4. Consideration of Poisson effect in thin plate element
p=
The procedure to consider Poisson's ratio in square thin plate element is discussed in this section. Additional terms to be introduced into the stiffness matrix to bring in Poisson effect are found by physical approach. For this, consider a plate meshed into square elements. An element and its surrounding elements are shown in Fig. 3(a)–(b). The
From Eqs. (6) and (9), it is seen that the total forces developed (p) are independent of the length of the sides of the plate element. Hence, the additional stiffness coefficients for rectangular plate element will be same as that of square plate element. Hence, further derivations are done for square element.
0
3
0
1
2
(9)
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 3. (a)–(e). (a) Plan of the plate elements along with the springs. (b) Group of elements (springs are not shown for clarity). (c) Unit deformation along w. (d) Unit deformation along θx. (e) Unit deformation along θy.
Additional stiffness (Ka1, Ka2, Ka3 and Ka4) to various elements, corresponding to rotational DOFs are given by Eq. (10).
When unit rotation is applied to element 0 in Fig. 4, forces develop in the surrounding elements and the element itself. There are 9 different cases with respect to the position of element 0 namely interior, left edge, right edge, top edge, top left corner, top right corner, bottom edge, bottom left corner and bottom right corner. The resultant force above and below the middle plane are indicated by solid line and dotted line respectively. The forces developed in the elements for unit rotation about X-axis are shown in Fig. 4(a)–(i) for the 9 cases respectively. The figures show the X-Y plane of the plate. The element numbers are indicated within the elements. If the force is directed towards the centre of the element, it is compressive and otherwise it is tensile. The forces developed in the elements for unit rotation about Y-axis are shown in Fig. 5(a)–(i) for the 9 cases respectively.
K a1 =
(m0 m0 )
K a2 =
( 0m m0 )
K a3 =
(m0
m 0
)
K a4 =
( 0m
m 0
)
(10) These stiffness matrices [Eq. (10)] are to be added to the stiffness matrices of the individual elements. The elements to which these stiffness matrices are to be added, depends upon the location of the displaced element and are given in Table 2. For example, if the element is located in the interior region, the stiffness matrix Ka1 should be added to the stiffness matrix of the elements 5 and 7. The stiffness matrix Ka4 is to be added to the stiffness 4
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 3. (continued)
matrix of elements 6 and 8. Similarly, the additional stiffness due to displacement of the DOFs in all the elements should be determined and added to the global stiffness matrix appropriately.
9th rows of the stiffness matrix in Eq. (1).
5. Displacements, strains, stresses and stress resultants
The normal strain (εx) and shear strains (γxy and γxz) at a point located at (l/2, y, z) with respect to the centre of element 0, which is at the joint of element 0 and element 2 shown in Fig. 7 can be determined using Eqs. (14)–(16).
5.2. Strain
The displacements, strains, stresses and stress resultants such as bending moments, twisting moment and shear forces at various joints are derived in this section. The size of the square element is l. For the derivations, a group of elements shown in Fig. 6 is considered. The element numbers are indicated within the elements.
=
x
=
z(
y2
u x y0 )
(14)
l
5.1. Displacements Displacement along X-direction (u), Y-direction (v) and Z-direction (w) with respect to the centre of element 0 for any point in between section CDD′C′ and GHH′G′ in Fig. 6 are given by Eqs. (11)–(13) respectively.
u=
xz ( l
xz v= ( l
w=
y2
x0
x (w 2 l
=
(11)
x 2)
x = ( l
(12)
xy ( l
x2
x 0)
+x
y0
y2
y0 )
+
=
2z (
xz
y0)
w0 ) +
xy
2z x0
2w
x y x2 )
l u w = + z x
1 (w 2 l
w0 ) +
y ( l
x 0)
x2
+
y0
(15)
Taking average shear strain along X-direction,
(13)
xz
The number in the subscript of the displacements represents the element number. These displacements are derived from the 7th, 8th and
=
1 (w 2 l
w0 ) +
y ( l
x2
x 0)
+
1 ( 2
y0
+
y2 )
(16)
For finding the normal strain along Y-direction (εy) at E, the average of the normal strain at the two nearby joints CE and EG is taken. The 5
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 4. Forces due to unit rotation about X-axis (a) interior (b) left edge (c) right edge (d) top edge (e) top left corner (f) top right corner (g) bottom edge (h) bottom left corner (i) bottom right corner.
normal strains at various joints in the group of elements are shown in Fig. 7. At the joint EF, normal strain in region EK and FK are given by Eq. (17). y
=
=
y1
y3
+
+
y2
2 y4
y1 =
z ( x0
x3 )
l
y2 =
(17) z ( x2
x7 )
l
y3 =
z ( x1
x0 )
l
y4 =
z ( x6
x2 )
l
=
E 1
2
(
x
+
y)
xy
=G
xy
xz
=G
xz
(18)
From Eq. (16), it is seen that the transverse shear strain (γxz) does not vary across the depth. Hence, as per Eq. (18), transverse shear stress (τxz) is also a constant across the depth. So, to obtain parabolic shear stress distribution, complementary shear stress (τzx) should be found. Consider an element shown in Fig. 8. The stresses along X-direction on the faces of the element are indicated in the figure. The top and bottom planes are devoid of stresses. The portion of the element at a distance z1 from the middle plane is cut out and the forces on its faces (F1, F2, F3 and F4) are indicated in Fig. 8. ‘A’ in Fig. 8 represents the area of the shaded section. The complementary shear stress at the shaded section is given by Eq. (19). Using equilibrium of forces,
(at EK)
(at FK)
2
x
where,
5.3. Stress The strains determined in section 5.2 is used to evaluate the normal stress (σx) and shear stresses (τxy, τxz) according to CPT given by Eq. (18).
zx
× A = F1
where, 6
F2 + F3
F4
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 5. Force due to unit rotation about Y-axis (a) interior (b) left edge (c) right edge (d) top edge (e) top left corner (f) top right corner (g) bottom edge (h) bottom left corner (i) bottom right corner. Table 2 Details of elements to which stiffness matrix are to be added. Position of element 0 Interior Left edge Right edge Top edge Top left corner Top right corner Bottom edge Bottom left corner Bottom right corner
Ka1 5, 7 7 5 5 0 5 7 7 0
Ka2 – 3 1 2 2 1 4 3 4
Ka3 – 1 3 4 1 4 2 2 3
F1 = Ka4
F2 =
6, 8 6 8 6 6 0 8 0 8
E 2)
4(1 E
2)
4(1
2(
y2
y0)
+
2(
y0
y4 )
+
2 2
(
x1
x3
+
x6
x7 )
h2 4
z12
(
x1
x3
+
x5
x8 )
h2 4
z12
F3 =
E ( 2(1 + )
y3
y0 )
h2 4
z12
F4 =
E ( 2(1 + )
y0
y1 )
h2 4
z12
zx
=
E 2) l2
2(1 (1
)(
y1
( 2
y2
y0
2 +
y0 y3 )
+
y4 )
h2 4
+
4
(
x5
+
x6
x7
+
x 8)
+
z12 (19)
From Eq. (19), it is seen that the complementary shear stress varies parabolically. 7
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 6. DOFs in a group of elements.
Eh3 ( y2 2) l 12(1 Eh3 = ( y2 2) l 12(1 =
y0 ) y0 )
+
+
2 2
(
(
=
2) l 12(1 Eh3 = ( 2) l 12(1
(
x 3)
x0
x0
x 3)
+
x0 h /2
My = Eh3
x0
dz
+
2 2
(
(
y0
+
y3
y0
+
x 7)
x3
x1
yz
h /2
+
+
x2
x2
+
y8)
y4
y2
y3
x 6)
+
y7 )
(EKK E ) (FKK F )
(ACC A ) (ECC E ) (21)
h /2
Mxy =
h /2
Gh3 = ( 6l
x2
xy z
dz
x 0)
(22)
The shear force in plate at the face 1 (Qx) and face 2 (Qy) are given by Eqs. (23) and (24) respectively. Fig. 7. Strain at element joints.
5.4. Stress resultants
1 = Gh (w 2 l
The bending moment (Mx) on the face 1 in the plate at EKK′E′ and FKK′F′ shown in Fig. 6 are given by Eq. (20). At face 2, the bending moment (My) at ACC′A′ and ECC′E′ are given by Eq. (21). Eq. (22) gives the twisting moment (Mxy) at face 1.
Mx =
h /2 h /2
x z dz
h /2
Qx =
h /2
1 w0 ) + ( 4
1 (w 3 l
w0 )
(20)
Fig. 8. Stresses along X-direction in an element. 8
h /2
1 ( 2
x0
dz x 0)
x2
h /2
Qy = = Gh
xz
yz
+
+
1 ( 2
y0
+
1 ( 4
y3
+
y2 )
(23)
dz x 3)
y0 )
(24)
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 9. (a) Case 1: Simply supported plate with sinusoidal loading. (b) Case 2: Simply supported plate with moments distributed along two opposite edges. (c) Case 3: Plate simply supported on two opposite edges and clamped on other edges subjected to UDL. Table 3 Results of simply supported plate with sinusoidal loading. Center (x = 0.6 m, y = 0.6 m)
W Mx My Mxy Qx Qy
AEM
CPT [14]
10.757 47,770.287 47,770.287 0 1.863 1.863
10.762 47,418.314 47,418.314 0 0 0
Quarter (x = 0.3 m, y = 0.3 m) Percentage difference (%) −0.05 0.74 0.74 0 – –
(a)
AEM
CPT [14]
5.384 24,154.588 24,263.209 12,739.826 94.317 94.316
5.381 23,709.157 23,709.157 12,766.469 95.493 95.493
(b)
Fig. 10. (a) Contour plot of deflection determined using AEM and CPT. (b) Deflected profile of the mid-section.
(a)
(b)
Fig. 11. (a) Normal stress distribution. (b) Complementary shear stress distribution.
9
Percentage difference (%) 0.05 1.88 2.34 −0.21 −1.23 −1.23
Structures 22 (2019) 1–12
D.L. Christy, et al.
Table 4 Results of problems in Case 2. ly/lx
0.5 0.75 1 1.5
w [×M0ls2/D] AEM
CPT [14]
0.0965 0.0623 0.0370 0.0283
0.0964 0.062 0.0368 0.0280
Percentage difference (%)
0.07 0.53 0.64 0.94
Mx [×M0]
Percentage difference (%)
AEM
CPT [14]
0.391 0.431 0.401 0.269
0.387 0.424 0.394 0.264
My [×M0]
1.11 1.77 1.84 1.76
Percentage difference (%)
AEM
CPT [14]
0.782 0.488 0.263 0.047
0.770 0.476 0.256 0.046
1.55 2.57 2.91 3.15
Table 5 Results of problems in Case 3. ly/lx
1 1.2 1.4 2 1/1.2 1/1.4 ½
w [×qls4/D]
Percentage difference (%)
AEM
CPT [14]
0.00190 0.00220 0.00238 0.00257 0.00317 0.00457 0.00826
0.00192 0.00223 0.00240 0.00260 0.00319 0.00460 0.00844
−0.95 −1.30 −0.78 −1.34 −0.67 −0.73 −2.16
Mx [×qls2]
Percentage difference (%)
AEM
CPT [14]
0.0244 0.0214 0.0187 0.0137 0.0378 0.0518 0.0864
0.0244 0.0215 0.0192 0.0142 0.0376 0.0514 0.0869
−0.04 −0.50 −2.44 −3.31 0.62 0.72 −0.61
(a)
My [×qls2]
Percentage difference (%)
AEM
CPT [14]
0.0332 0.0372 0.0395 0.0412 0.039954 0.044379 0.04648
0.0332 0.0375 0.0399 0.0420 0.0400 0.0448 0.0474
−0.07 −0.68 −1.08 −2.02 −0.12 −0.94 −1.94
(b)
(c) Fig. 12. (a) Deflection vs number of elements. (b) Deflection vs number of springs. (c) Percentage difference vs computational time.
element of length l and thickness h is given by Eq. (25). ρ is the mass density of the material.
Table 6 Fundamental natural frequency of plates. ly/lx
Simply supported plate AEM
1 0.8 0.6 0.5
640.314 822.317 1218.517 1620.167
Clamped plate
Analytical [15]
Percentage difference (%)
628.141 804.806 1186.489 1570.353
1.94 2.18 2.70 3.17
AEM
1161.400 1509.249 2330.866 3192.366
Analytical [15]
Percentage difference (%)
1156.683 1501.513 2314.100 3165.108
0.41 0.52 0.72 0.86
pl 2 h 0 0 M= 0 pl4 h/6 0 0 0 pl4 h/6
(25)
7. Verification of numerical model The thin plate element is verified by analyzing plates of different boundary conditions and loading. Both static and dynamic analysis using square elements are performed. All the analyses are done using MATLAB code. The modulus of elasticity and Poisson's ratio of the material used are 2 × 105 N/mm2 and 0.3 respectively. The thickness of the plates is 30 mm. The elements are connected using 100 springs in
6. Mass matrix The mass matrix (M) of an element is determined by lumping the mass of the element at the centre. The mass matrix of a square thin plate 10
Structures 22 (2019) 1–12
D.L. Christy, et al.
(a)
(b)
Fig. 13. Contour plot of first mode of vibration of simply supported plate determined (a) using AEM (b) analytically.
(a)
(b)
Fig. 14. Contour plot of first mode of vibration of clamped plate determined (a) using AEM (b) analytically.
Fig. 15. (a) Second mode of simply supported plate. (b) Third mode of fixed plate.
each face on all the directions.
shows the deflected profile of the mid-section. The distribution of normal stress and complementary shear stress at the section at quarter point of the plate are shown in Fig. 11. In the second case (Fig. 9(b)), a simply supported plate with moments distributed along the edges parallel to X-axis is taken. For the study, a moment of 1 Nm is adopted. The deflection and bending moments at the centre of the plate are determined for various length to width ratios and are presented in Table 4. The smaller dimension is maintained as 1.2 m for all the problems in the second case. The third case (Fig. 9(c)) is a rectangular plate with two opposite edges simply supported and other two edges parallel to the X-axis clamped. The plate is subjected to uniformly distributed load (UDL) of intensity 10 kN/m2. The analysis is done for various length to width ratios keeping the larger dimension as 1.2 m. The deflection and bending moments at the centre of the plate of all the problems are presented in Table 5. For all the cases, the values of deflection, bending moments and
7.1. Static analysis Three cases shown in Fig. 9(a)–(c) are considered. lx and ly are the length of the plate along X-axis and Y-axis respectively. For the cases, the size of the element is approximately 7.45 mm. The first case (Fig. 9(a)) is a square plate of length 1.2 m, simply supported on all edges with sinusoidal loading as per Eq. (26). The intensity of loading at the centre (q0) is 1 N/mm2.
q = q0 sin
y x sin lx ly
(26)
The deflection (w), bending moments, twisting moment and shear forces at centre and quarter point of the plate obtained using AEM are given in Table 3. The contour plot of the deflected shape of the plate determined using AEM and CPT are shown in Fig. 10(a). Fig. 10(b) 11
Structures 22 (2019) 1–12
D.L. Christy, et al.
shear forces according to CPT are also given in Tables 3, 4 and 5. The deflection, moments and shear forces as per CPT in Table 3 are calculated using Eqs. (e), (f) and (g) in section 27 of [14]. ls in Tables 4 and 5 represents the smaller span. The deflection coefficients and bending moment coefficients of CPT given in Tables 4 and 5 are obtained from Tables 28 and 29 of [14]. Figs. 10–11 and Tables 3–5 shows that AEM could predict reliable values of deflection, stresses, bending moments, twisting moment and shear forces.
8. Conclusions The following results are arrived from the study:▪ Thin plate element for AEM is developed and validated. ▪ The method to incorporate Poisson's ratio in thin plate element is illustrated. ▪ The equations for calculating displacements, strains, stresses, bending moments, twisting moment and shear forces are derived. ▪ Thin plate element in AEM predicted the deflection, normal stress distribution, complementary shear stress distribution, bending moments, twisting moment and shear forces with acceptable accuracy. ▪ The deflection converged as the number of elements and the number of springs increased. The percentage error in deflection reduces to the acceptable limit at a reasonable computational time. ▪ The fundamental natural frequency and mode shape of plates are determined with reasonable accuracy, with the help of the thin plate element developed.
7.2. Convergence study To study the convergence of AEM with respect to the number of elements and number of springs, a simply supported plate with sinusoidal loading (Case 1 of Section 7.1) is considered. The variation of deflection due to the increase in number of elements and springs are shown in Fig. 12(a) and (b) respectively. The percentage difference of the central deflection obtained by AEM from that determined analytically using CPT are plotted against computational time in Fig. 12(c). Fig. 12(a) and (b) shows that as the number of elements and springs increases, the deflection converges towards that of CPT. As usual, convergence is mainly affected by the number of elements rather than the number of springs. From Fig. 12(c), it is seen that the percentage error in deflection reduces to the acceptable limit at a reasonable computational time.
The study mainly focussed on the development of thin plate element in AEM. This can be extended to the analysis of plates having different property and geometry. References [1] Meguro K, Tagel-Din H. A new efficient technique for fracture analysis of structures. Bull of Earthq Resistant Struct Res Cent, IIS, Univ of Tokyo 1997 Mar;30:103–16. [2] Meguro K, Tagel-Din H. A new simplified and efficient technique for fracture behavior analysis of concrete structures. Proceedings of the 3rd international conference on fracture mechanics of concrete and concrete structures (FRAMCOS-3). 1998 Oct.. p. 911–20. [3] Tagel-Din H, Meguro K. Nonlinear simulation of RC structures using applied element method. Struct Eng Earthq Eng 2000 Jul 21;17(2):13–24. [4] Tagel-Din H, Meguro K. Consideration of Poisson’s ratio effect in structural analysis using elements with three degrees of freedom. Bull of Earthq Resistant Cent 1998;31:41–50. [5] Meguro K, Tagel-Din H. Applied element method for structural analysis: theory and application for linear materials. Struct. Eng. Earthq. Eng. 2000;17(1):21–35. [6] Tagel-Din H, Meguro K. Applied element method for dynamic large deformation analysis of structures. Struct. Eng. Earthq. Eng. 2000 Oct 21;17(2):215–24. [7] Meguro K, Tagel-Din H. Applied element method used for large displacement structural analysis. J Nat Disaster Sci 2002 Jun 1;24(1):25–34. [8] Khalil AA. Enhanced modeling of steel structures for progressive collapse analysis using the applied element method. J Perform Constr Facil 2012 Dec;26(6):766–79. https://doi.org/10.1061/(ASCE)CF.1943-5509.0000267. [9] Cismasiu C, Ramos AP, Moldovan ID, Ferreira DF, Filho JB. Applied element method simulation of experimental failure modes in RC shear walls. Comput Concr 2017;19(4). [10] Malomo D, Pinho R, Penna A. Using the applied element method for modelling calcium-silicate brick masonry subjected to in-plane cyclic loading. Earthq Eng Struct Dyn 2018 Jun;47(7):1610–30. [11] Christy DL, Pillai TM, Nagarajan P. Analysis of concrete beams using applied element method. IOP conference series: materials science and engineering. 330(1). 2018 Mar. p. 1–7. [12] Shakeri A, Bargi K. Use of applied element method for structural analysis. KSCE J Civ Eng 2015 Jul 1;19(5):1375–84. https://doi.org/10.1007/s12205-015-0625-4. [13] Prashidha Kharel's Blog [Internet]. Formulating the applied element method: linear 2D part(1) March 4 - [cited 2019 April 18]. Available from:. http://prashidha. blogspot.in/2014/03/formulating-applied-element-method.html; 2014. [14] Timoshenko SP, Woinowsky-Krieger S. Theory of plates and shells. McGraw-Hill; 1959. [15] Ventsel E, Krauthammer T. Thin plates and shells: theory: analysis, and applications. CRC press; 2001.
7.3. Dynamic analysis The dynamic analysis of plates simply supported on all edges and clamped on all edges are studied in this section. Plates of different length to width ratios are analysed by keeping the larger dimension as 1.2 m and element size as approximately 11.88 mm. The fundamental natural frequency (ω) obtained analytically [15] are given by Eqs. (27) and (28) for simply supported and clamped plates respectively. The corresponding values determined using AEM are presented in Table 6.
=
2
1 1 + 2 lx 2 ly
= 4 25
D m
1 1 28 + 4 + 2 2 2 2 lx 4 ly 3 7 lx l y
(27)
D h
(28)
Eh3 2) . 12(1
where, D = Here, ρ and m are the mass density and mass per unit area respectively. Table 6 shows that the natural frequency determined using AEM are in line with that derived by analytical means (based on CPT). AEM predicts the fundamental natural frequency with very less error. The contour plot of the first mode of vibration of simply supported plate, determined using AEM and analytically are shown in Fig. 13(a) and (b) respectively. The corresponding plots of clamped plate are shown in Fig. 14(a) and (b). Fig. 15(a) and (b) shows the surface plot of second mode of simply supported plate and third mode of fixed plate respectively. Figs. 13–15 shows that the thin plate element in AEM is capable of predicting the mode shapes of plates.
12
Contents lists available at ScienceDirect
Structures journal homepage: www.elsevier.com/locate/structures
Thin plate element for applied element method D. Lincy Christy , T.M. Madhavan Pillai, Praveen Nagarajan ⁎
T
Department of Civil Engineering, National Institute of Technology Calicut, India
ARTICLE INFO
ABSTRACT
Keywords: Applied element method Thin plate element Stiffness matrix Poisson effect
Applied Element Method (AEM) is an efficient method for numerical analysis of structures. It is capable of analyzing the structure up to collapse. In AEM, the structure is modelled as a group of rigid elements connected by springs. It has the benefit of lesser processing time and memory requirement when compared to Finite Element Method. It is also able to track crack propagation more realistically. To extend the advantages of AEM to the analysis of thin plate structures, thin plate element has to be developed. This is done in this paper. Being a rigid body method, AEM does not consider Poisson's ratio in its basic formulation. As Poisson's ratio plays a major role in the behavior of plate structures, the method to bring in Poisson effect in thin plate element is discussed. The procedure to determine displacements, strains, stresses and stress resultants are explained. The mass matrix for thin plate element is also found. To validate the thin plate element, both static and dynamic analyses are carried out. The analyses gave reliable results which show the efficiency of the thin plate element of AEM.
1. Introduction Applied Element Method (AEM) is an outcome of discrete element method used for structural analysis. Meguro and Tagel-Din introduced this numerical method in 1997 [1]. They advantageously used this method for the fracture analysis of reinforced concrete structures [1,2]. They detailed the procedure to determine crack pattern using AEM [3]. The procedure to consider Poisson effect in two-dimensional (2D) element was also discussed by them [4,5]. With the help of this, they analysed elastic prism subjected to in-plane loading. In 2000, a method to find mass matrix for 2D element in AEM was devised [6]. The geometric nonlinearity was handled by residual force vector generated due to geometric changes instead of geometric stiffness matrix [6,7]. Khalil analysed moment-resisting steel frame using AEM for progressive collapse assessment [8]. Cismasiu et al. studied the failure of reinforced shear wall using AEM in 2017 [9]. In 2018, the study on the behavior of brick masonry wall under in-plane cyclic loading [10] and the convergence study on plain concrete beams [11] were done. As springs are used to model the structure, the material property can be easily varied within the structure in AEM. This is advantageous for the analysis of heterogeneous materials such as functionally graded materials. In 2015, Shakeri and Bargi did static analysis of rectangular column, beam and slab to study the effect of number of springs and elements on the accuracy of results [12]. The Poisson's ratio was taken as 0 for all the cases. The analysis of rectangular slab with transverse loading was
⁎
accomplished using three-dimensional (3D) elements. They plotted the normal stress and shear stress distribution in rectangular beam. Dynamic analysis was also carried out with the help of Newmark-β method. From the literature it is seen that, analysis of thin plate with transverse loading and considering Poisson's ratio has not been done. In this paper, a thin plate element is developed for AEM. The stiffness matrix for the analysis of thin plate is derived from the stiffness matrix of a pair of springs in 3D analysis. The effect of Poisson's ratio on the behavior of plate structures is significant. So, the technique to consider the effect of Poisson's ratio is illustrated using the analysis with thin plate element. The method to calculate displacements, strains, stresses, bending moments, twisting moment and shear force are described. Static and dynamic analyses of plates are performed to verify the thin plate element developed. Plates of different boundary condition and loading are analysed to find the deflection, moments and shear force. The normal stress distribution and shear stress distribution are also plotted. Convergence study with respect to number of elements and number of springs is done. All the results are compared with that of Classical Plate Theory (CPT). The mass matrix of thin plate element is determined to carry out dynamic analysis. Using this, the fundamental natural frequency of plates simply supported on all edges and clamped on all edges are found. The contour plot of the mode shapes are also drawn and verified.
Corresponding author. E-mail addresses: [email protected] (D.L. Christy), [email protected] (T.M.M. Pillai), [email protected] (P. Nagarajan).
https://doi.org/10.1016/j.istruc.2019.07.010 Received 15 May 2019; Received in revised form 22 July 2019; Accepted 23 July 2019 2352-0124/ © 2019 Institution of Structural Engineers. Published by Elsevier Ltd. All rights reserved.
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 1. (a) 3D element. (b) Connection of 3D elements.
2. Applied element method
3. Thin plate element for AEM
AEM is used for the numerical analysis of structures. In this method, the structure is divided into rigid elements connected by pairs of springs. The degrees of freedom (DOFs) are located at the centre of the elements. There are 6 DOFs for a 3D element as shown in Fig. 1(a). They are 3 translations and 3 rotations. In 3D analysis, a pair of springs consists of one spring to carry normal force (A) and two springs to take shear forces (B and C) along the transverse directions (Fig. 1(b)). Every pair of springs has a stiffness matrix corresponding to the DOFs of the elements connected by them. The stiffness matrix of a pair of springs connecting 2D elements is given in [13]. The stiffness matrix for a pair of springs (Ks) located at (x, y, z) from the centroid of the first element shown in Fig. 1(b) is given by Eq. (1).
For the analysis of plates, 2D element is sufficient, if only in-plane loading is applied. Anyway, when transverse loading is applied, thin plate element for transverse loading has to be developed. To arrive at the behavior of thin plate using 3D element, 3 major modifications are to be made. They are discussed in this section.
Ks =
kn 0 0 0 knz -kny -kn 0 ks 0 -ksz 0 ks x 0 0 0 ks ksy -ksx 0 0 0 -ksz ksy ks(y2+ z2) -ksxy -ksxz 0 knz 0 -ksx -ksyx knz2+ksx2 -knyz -knz -kny ksx 0 -kszx -knzy kny2+ksx2 kny
3.1. Modification 1: value of E′ The value of E′ depends upon the constitutive relation of the problem considered. For normal spring, strain along its direction alone is important. Hence, the value of E′ can be evaluated by ignoring other normal strains. The value of E′ for different types of problems is given in
0 0 0 -knz kn y -ks 0 ksz 0 ksx 0 -ks -ksy -ksx 0 ksz -ksy -ks(y2+ z2) -ksxy -ksxz 0 ksx ksyx -knz2+ksx2 knyz -ksx 0 kszx knzy -kny2+ksx2
-kn 0 0 0 -knz kn y kn 0 0 0 knz -kny 0 -ks 0 ksz 0 -ksx 0 ks 0 -ksz 0 -ksx 0 0 -ks -ksy ksx 0 0 0 ks ks y ksx 0 0 ksz -ksy -ks(y2+ z2) ksxy ksxz 0 -ksz ksy ks(y2+ z2) ksxy ksxz -knz 0 -ksx -ksyx -knz2+ksx2 knyz knz 0 ksx ksyx knz2+ksx2 -knyz kny ksx 0 -kszx knzy -kny2+ksx2 -kny -ksx 0 kszx -knzy kny2+ksx2
The stiffnesses of the normal spring (kn) and shear springs (ks) are given by Eq. (2).
kn =
E As GAs ; ks = ;E = l l
(1)
Table 1. Here, E and ν are the modulus of elasticity and Poisson's ratio respectively. From Table 1, it is seen that the stiffness of thin plate is more than that of beam by a factor 1/(1−ν2). The additional rigidity is due to the Poisson's ratio.
(2)
where, G, l and As are the modulus of rigidity, length of the element and area of cross-section of the portion of the element treated as the spring respectively. The normal stress and normal strain along the direction of normal spring are denoted by σ and ɛ respectively. For all the springs, x = l/2. The 3D element is the skeleton of AEM. All other types of elements can be obtained by reducing or modifying the 3D element. Similarly, the stiffness matrix for any analysis can be derived from the stiffness matrix for 3D analysis. In the case of 2D analysis, one of the shear springs is not necessary. Also, the DOFs can be reduced to 3; 2 translations and 1 rotation. The size of stiffness matrix now reduces to 6 × 6 by removing the columns and rows corresponding to the unwanted DOFs. When it comes to one-dimensional (1D) analysis, only normal spring is to be considered and a single translational DOF becomes sufficient. The size of stiffness matrix becomes as small as 2 × 2.
3.2. Modification 2: degrees of freedom required for analysis The basic behavior of thin plates helps to simplify the analysis. According to Kirchhoff's assumptions, the translations along the middle plane is 0. Therefore, if one element is taken along the depth, the first 2 DOFs [u and v in Fig. 1(a)] are not necessary. Element rotation about Zaxis is given by Eq. (3). Table 1 Value of E′ for different problems. 1D/2D E
2
Plane stress/thin plate E (1
2)
Plane strain/3D E (1 ) (1 + )(1 2 )
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 2. (a) Thin plate element (b) connection of thin plate elements.
z
=
1 2
v x
u =0 y
displaced profile of the group of elements when unit deformation is applied to the central element along its DOFs (w, θx and θy) are shown in Fig. 3(c)–(e) respectively. Unit deformation along the depth direction, does not induce any force in the normal springs and hence Poisson effect does not come into picture. However, rotation about X-axis and Y-axis strains the normal springs along Y-axis and X-axis respectively. Due to unit rotation about Y-axis (θy = 1), the line AB on the right face of the central element is displaced to A′B′ as shown in Fig. 3(e). This creates compressive strain and tensile strain above and below the middle plane respectively. The direct forces generated in the direction of X-axis due to these strains are already taken into account in the stiffness matrix. In addition to the direct forces, transverse forces develop along Y-axis due to Poisson effect. These are not considered in the stiffness matrix. The normal stress along Y-axis (σy) is:-
(3)
For thin plates, as u and v are ignored, θz turns out to be 0. Therefore, the 6th DOF also can be neglected. Hence, translation along Z-axis (w), rotation about X-axis (θx) and rotation about Y-axis (θy) are the DOFs necessary for thin plate analysis. The DOFs in thin plate element are shown in Fig. 2(a). All the 3 springs are required to connect the thin plate elements as shown in Fig. 2(b). The stiffness matrix of a pair of springs now reduces to a size of 6 × 6 [Eq. (4)].
Ks =
ks
ks y
ks l/2
ks
ks y
ks l/2
ks y
ks (y 2 + z 2)
ks ly /2
ks y
ks (y 2 + z 2)
ks ly /2
ks l/2 ks
ks ly /2 ks y
ks y
ks (y 2 + z 2)
ks l/2
ks ly /2
kn z 2 + ks l2/4 ks l/2
ks l/2 ks
ks ly /2 ks y
ks ly /2
ks y
ks (y 2 + z 2)
ks ly /2
ks ly /2
kn z 2 + ks l2/4
kn z 2 + ks l2/4 ks l/2
kn z 2 + ks l2/4 ks l/2
(4)
To get the complete stiffness matrix at a face of the plate element, the stiffness matrix of individual pairs of springs should be summed up. If infinite number of springs are adopted in the section, the springs will be of infinitesimal area. Then, the sum of the stiffness matrix of all pairs of springs (K) can be obtained by integrating the stiffness matrix across the section as given by Eq. (5). K=
K=
k1 =
y
=
0 k2 0 0 k2 0
k4 0 k3 k4 0 k5
k1 0 k4 k1 0 k4
0 k2 0 0 k2 0
x
y
Gbh k5 = 2
(
y
+
x)
= z / l and
y
=0
Therefore,
k4 0 k5 k4 0 k3
Gbh Gbh (b2 + h2) bh (Eb2 + 3Gl2) k4 = k2 = k3 = l 12l 12l
2
For the shaded quarter portion in Fig. 3(e), the normal strain along X-axis (ɛx) and normal strain along Y-axis (ɛy) are:-
K s dA k1 0 k4 k1 0 k4
E 1
=
E
z l
2
1
The total compressive force (p) and the total tensile force for a quarter portion [shaded region in Fig. 3(e)] of the element will be equal and is given by Eq. (6). The lever arm (r) between total compressive force and total tensile force is given by Eq. (7). These forces in the quarter part of the plate element creates a moment (m) at the centre of the plate given by Eq. (8).
bh (Eb2 3Gl2) 12l
(5) where, b and h are the width and thickness of the plate element. This stiffness matrix [Eq. (5)] can be directly used at a section instead of summing up the individual stiffness matrices of all pairs of springs.
h /2 l /2
p=
y 0
=
3.3. Modification 3: including the effect of Poisson's ratio
Eh2 2) 16(1
(6)
2 h 3
(7)
m = pr
(8)
r=
Poisson's ratio has a major role in the behavior of plates. Hence, Poisson effect must be invariably considered in the analysis of thin plates. The methodology to bring in Poisson effect in 2D element is discussed by Meguro and Tagel-Din [4,5]. For this, they introduced additional terms into the stiffness matrix. According to the Modification 2, a thin plate element has 3 DOFs which are different from those in 2D element. Hence, a method has to be devised to include Poisson effect in thin plate element. This is discussed in the next section.
dxdz
0
In the case of rectangular plate element, the total compressive force and total tensile force (p) are given by Eq. (9). h /2 b /2
p= 0 0 h /2 l /2
E 2
1 E
z Eh2 dydz = for 2) b 16(1
x
=1
z Eh2 dxdz = for 2) l 16(1
y
=1
4. Consideration of Poisson effect in thin plate element
p=
The procedure to consider Poisson's ratio in square thin plate element is discussed in this section. Additional terms to be introduced into the stiffness matrix to bring in Poisson effect are found by physical approach. For this, consider a plate meshed into square elements. An element and its surrounding elements are shown in Fig. 3(a)–(b). The
From Eqs. (6) and (9), it is seen that the total forces developed (p) are independent of the length of the sides of the plate element. Hence, the additional stiffness coefficients for rectangular plate element will be same as that of square plate element. Hence, further derivations are done for square element.
0
3
0
1
2
(9)
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 3. (a)–(e). (a) Plan of the plate elements along with the springs. (b) Group of elements (springs are not shown for clarity). (c) Unit deformation along w. (d) Unit deformation along θx. (e) Unit deformation along θy.
Additional stiffness (Ka1, Ka2, Ka3 and Ka4) to various elements, corresponding to rotational DOFs are given by Eq. (10).
When unit rotation is applied to element 0 in Fig. 4, forces develop in the surrounding elements and the element itself. There are 9 different cases with respect to the position of element 0 namely interior, left edge, right edge, top edge, top left corner, top right corner, bottom edge, bottom left corner and bottom right corner. The resultant force above and below the middle plane are indicated by solid line and dotted line respectively. The forces developed in the elements for unit rotation about X-axis are shown in Fig. 4(a)–(i) for the 9 cases respectively. The figures show the X-Y plane of the plate. The element numbers are indicated within the elements. If the force is directed towards the centre of the element, it is compressive and otherwise it is tensile. The forces developed in the elements for unit rotation about Y-axis are shown in Fig. 5(a)–(i) for the 9 cases respectively.
K a1 =
(m0 m0 )
K a2 =
( 0m m0 )
K a3 =
(m0
m 0
)
K a4 =
( 0m
m 0
)
(10) These stiffness matrices [Eq. (10)] are to be added to the stiffness matrices of the individual elements. The elements to which these stiffness matrices are to be added, depends upon the location of the displaced element and are given in Table 2. For example, if the element is located in the interior region, the stiffness matrix Ka1 should be added to the stiffness matrix of the elements 5 and 7. The stiffness matrix Ka4 is to be added to the stiffness 4
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 3. (continued)
matrix of elements 6 and 8. Similarly, the additional stiffness due to displacement of the DOFs in all the elements should be determined and added to the global stiffness matrix appropriately.
9th rows of the stiffness matrix in Eq. (1).
5. Displacements, strains, stresses and stress resultants
The normal strain (εx) and shear strains (γxy and γxz) at a point located at (l/2, y, z) with respect to the centre of element 0, which is at the joint of element 0 and element 2 shown in Fig. 7 can be determined using Eqs. (14)–(16).
5.2. Strain
The displacements, strains, stresses and stress resultants such as bending moments, twisting moment and shear forces at various joints are derived in this section. The size of the square element is l. For the derivations, a group of elements shown in Fig. 6 is considered. The element numbers are indicated within the elements.
=
x
=
z(
y2
u x y0 )
(14)
l
5.1. Displacements Displacement along X-direction (u), Y-direction (v) and Z-direction (w) with respect to the centre of element 0 for any point in between section CDD′C′ and GHH′G′ in Fig. 6 are given by Eqs. (11)–(13) respectively.
u=
xz ( l
xz v= ( l
w=
y2
x0
x (w 2 l
=
(11)
x 2)
x = ( l
(12)
xy ( l
x2
x 0)
+x
y0
y2
y0 )
+
=
2z (
xz
y0)
w0 ) +
xy
2z x0
2w
x y x2 )
l u w = + z x
1 (w 2 l
w0 ) +
y ( l
x 0)
x2
+
y0
(15)
Taking average shear strain along X-direction,
(13)
xz
The number in the subscript of the displacements represents the element number. These displacements are derived from the 7th, 8th and
=
1 (w 2 l
w0 ) +
y ( l
x2
x 0)
+
1 ( 2
y0
+
y2 )
(16)
For finding the normal strain along Y-direction (εy) at E, the average of the normal strain at the two nearby joints CE and EG is taken. The 5
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 4. Forces due to unit rotation about X-axis (a) interior (b) left edge (c) right edge (d) top edge (e) top left corner (f) top right corner (g) bottom edge (h) bottom left corner (i) bottom right corner.
normal strains at various joints in the group of elements are shown in Fig. 7. At the joint EF, normal strain in region EK and FK are given by Eq. (17). y
=
=
y1
y3
+
+
y2
2 y4
y1 =
z ( x0
x3 )
l
y2 =
(17) z ( x2
x7 )
l
y3 =
z ( x1
x0 )
l
y4 =
z ( x6
x2 )
l
=
E 1
2
(
x
+
y)
xy
=G
xy
xz
=G
xz
(18)
From Eq. (16), it is seen that the transverse shear strain (γxz) does not vary across the depth. Hence, as per Eq. (18), transverse shear stress (τxz) is also a constant across the depth. So, to obtain parabolic shear stress distribution, complementary shear stress (τzx) should be found. Consider an element shown in Fig. 8. The stresses along X-direction on the faces of the element are indicated in the figure. The top and bottom planes are devoid of stresses. The portion of the element at a distance z1 from the middle plane is cut out and the forces on its faces (F1, F2, F3 and F4) are indicated in Fig. 8. ‘A’ in Fig. 8 represents the area of the shaded section. The complementary shear stress at the shaded section is given by Eq. (19). Using equilibrium of forces,
(at EK)
(at FK)
2
x
where,
5.3. Stress The strains determined in section 5.2 is used to evaluate the normal stress (σx) and shear stresses (τxy, τxz) according to CPT given by Eq. (18).
zx
× A = F1
where, 6
F2 + F3
F4
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 5. Force due to unit rotation about Y-axis (a) interior (b) left edge (c) right edge (d) top edge (e) top left corner (f) top right corner (g) bottom edge (h) bottom left corner (i) bottom right corner. Table 2 Details of elements to which stiffness matrix are to be added. Position of element 0 Interior Left edge Right edge Top edge Top left corner Top right corner Bottom edge Bottom left corner Bottom right corner
Ka1 5, 7 7 5 5 0 5 7 7 0
Ka2 – 3 1 2 2 1 4 3 4
Ka3 – 1 3 4 1 4 2 2 3
F1 = Ka4
F2 =
6, 8 6 8 6 6 0 8 0 8
E 2)
4(1 E
2)
4(1
2(
y2
y0)
+
2(
y0
y4 )
+
2 2
(
x1
x3
+
x6
x7 )
h2 4
z12
(
x1
x3
+
x5
x8 )
h2 4
z12
F3 =
E ( 2(1 + )
y3
y0 )
h2 4
z12
F4 =
E ( 2(1 + )
y0
y1 )
h2 4
z12
zx
=
E 2) l2
2(1 (1
)(
y1
( 2
y2
y0
2 +
y0 y3 )
+
y4 )
h2 4
+
4
(
x5
+
x6
x7
+
x 8)
+
z12 (19)
From Eq. (19), it is seen that the complementary shear stress varies parabolically. 7
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 6. DOFs in a group of elements.
Eh3 ( y2 2) l 12(1 Eh3 = ( y2 2) l 12(1 =
y0 ) y0 )
+
+
2 2
(
(
=
2) l 12(1 Eh3 = ( 2) l 12(1
(
x 3)
x0
x0
x 3)
+
x0 h /2
My = Eh3
x0
dz
+
2 2
(
(
y0
+
y3
y0
+
x 7)
x3
x1
yz
h /2
+
+
x2
x2
+
y8)
y4
y2
y3
x 6)
+
y7 )
(EKK E ) (FKK F )
(ACC A ) (ECC E ) (21)
h /2
Mxy =
h /2
Gh3 = ( 6l
x2
xy z
dz
x 0)
(22)
The shear force in plate at the face 1 (Qx) and face 2 (Qy) are given by Eqs. (23) and (24) respectively. Fig. 7. Strain at element joints.
5.4. Stress resultants
1 = Gh (w 2 l
The bending moment (Mx) on the face 1 in the plate at EKK′E′ and FKK′F′ shown in Fig. 6 are given by Eq. (20). At face 2, the bending moment (My) at ACC′A′ and ECC′E′ are given by Eq. (21). Eq. (22) gives the twisting moment (Mxy) at face 1.
Mx =
h /2 h /2
x z dz
h /2
Qx =
h /2
1 w0 ) + ( 4
1 (w 3 l
w0 )
(20)
Fig. 8. Stresses along X-direction in an element. 8
h /2
1 ( 2
x0
dz x 0)
x2
h /2
Qy = = Gh
xz
yz
+
+
1 ( 2
y0
+
1 ( 4
y3
+
y2 )
(23)
dz x 3)
y0 )
(24)
Structures 22 (2019) 1–12
D.L. Christy, et al.
Fig. 9. (a) Case 1: Simply supported plate with sinusoidal loading. (b) Case 2: Simply supported plate with moments distributed along two opposite edges. (c) Case 3: Plate simply supported on two opposite edges and clamped on other edges subjected to UDL. Table 3 Results of simply supported plate with sinusoidal loading. Center (x = 0.6 m, y = 0.6 m)
W Mx My Mxy Qx Qy
AEM
CPT [14]
10.757 47,770.287 47,770.287 0 1.863 1.863
10.762 47,418.314 47,418.314 0 0 0
Quarter (x = 0.3 m, y = 0.3 m) Percentage difference (%) −0.05 0.74 0.74 0 – –
(a)
AEM
CPT [14]
5.384 24,154.588 24,263.209 12,739.826 94.317 94.316
5.381 23,709.157 23,709.157 12,766.469 95.493 95.493
(b)
Fig. 10. (a) Contour plot of deflection determined using AEM and CPT. (b) Deflected profile of the mid-section.
(a)
(b)
Fig. 11. (a) Normal stress distribution. (b) Complementary shear stress distribution.
9
Percentage difference (%) 0.05 1.88 2.34 −0.21 −1.23 −1.23
Structures 22 (2019) 1–12
D.L. Christy, et al.
Table 4 Results of problems in Case 2. ly/lx
0.5 0.75 1 1.5
w [×M0ls2/D] AEM
CPT [14]
0.0965 0.0623 0.0370 0.0283
0.0964 0.062 0.0368 0.0280
Percentage difference (%)
0.07 0.53 0.64 0.94
Mx [×M0]
Percentage difference (%)
AEM
CPT [14]
0.391 0.431 0.401 0.269
0.387 0.424 0.394 0.264
My [×M0]
1.11 1.77 1.84 1.76
Percentage difference (%)
AEM
CPT [14]
0.782 0.488 0.263 0.047
0.770 0.476 0.256 0.046
1.55 2.57 2.91 3.15
Table 5 Results of problems in Case 3. ly/lx
1 1.2 1.4 2 1/1.2 1/1.4 ½
w [×qls4/D]
Percentage difference (%)
AEM
CPT [14]
0.00190 0.00220 0.00238 0.00257 0.00317 0.00457 0.00826
0.00192 0.00223 0.00240 0.00260 0.00319 0.00460 0.00844
−0.95 −1.30 −0.78 −1.34 −0.67 −0.73 −2.16
Mx [×qls2]
Percentage difference (%)
AEM
CPT [14]
0.0244 0.0214 0.0187 0.0137 0.0378 0.0518 0.0864
0.0244 0.0215 0.0192 0.0142 0.0376 0.0514 0.0869
−0.04 −0.50 −2.44 −3.31 0.62 0.72 −0.61
(a)
My [×qls2]
Percentage difference (%)
AEM
CPT [14]
0.0332 0.0372 0.0395 0.0412 0.039954 0.044379 0.04648
0.0332 0.0375 0.0399 0.0420 0.0400 0.0448 0.0474
−0.07 −0.68 −1.08 −2.02 −0.12 −0.94 −1.94
(b)
(c) Fig. 12. (a) Deflection vs number of elements. (b) Deflection vs number of springs. (c) Percentage difference vs computational time.
element of length l and thickness h is given by Eq. (25). ρ is the mass density of the material.
Table 6 Fundamental natural frequency of plates. ly/lx
Simply supported plate AEM
1 0.8 0.6 0.5
640.314 822.317 1218.517 1620.167
Clamped plate
Analytical [15]
Percentage difference (%)
628.141 804.806 1186.489 1570.353
1.94 2.18 2.70 3.17
AEM
1161.400 1509.249 2330.866 3192.366
Analytical [15]
Percentage difference (%)
1156.683 1501.513 2314.100 3165.108
0.41 0.52 0.72 0.86
pl 2 h 0 0 M= 0 pl4 h/6 0 0 0 pl4 h/6
(25)
7. Verification of numerical model The thin plate element is verified by analyzing plates of different boundary conditions and loading. Both static and dynamic analysis using square elements are performed. All the analyses are done using MATLAB code. The modulus of elasticity and Poisson's ratio of the material used are 2 × 105 N/mm2 and 0.3 respectively. The thickness of the plates is 30 mm. The elements are connected using 100 springs in
6. Mass matrix The mass matrix (M) of an element is determined by lumping the mass of the element at the centre. The mass matrix of a square thin plate 10
Structures 22 (2019) 1–12
D.L. Christy, et al.
(a)
(b)
Fig. 13. Contour plot of first mode of vibration of simply supported plate determined (a) using AEM (b) analytically.
(a)
(b)
Fig. 14. Contour plot of first mode of vibration of clamped plate determined (a) using AEM (b) analytically.
Fig. 15. (a) Second mode of simply supported plate. (b) Third mode of fixed plate.
each face on all the directions.
shows the deflected profile of the mid-section. The distribution of normal stress and complementary shear stress at the section at quarter point of the plate are shown in Fig. 11. In the second case (Fig. 9(b)), a simply supported plate with moments distributed along the edges parallel to X-axis is taken. For the study, a moment of 1 Nm is adopted. The deflection and bending moments at the centre of the plate are determined for various length to width ratios and are presented in Table 4. The smaller dimension is maintained as 1.2 m for all the problems in the second case. The third case (Fig. 9(c)) is a rectangular plate with two opposite edges simply supported and other two edges parallel to the X-axis clamped. The plate is subjected to uniformly distributed load (UDL) of intensity 10 kN/m2. The analysis is done for various length to width ratios keeping the larger dimension as 1.2 m. The deflection and bending moments at the centre of the plate of all the problems are presented in Table 5. For all the cases, the values of deflection, bending moments and
7.1. Static analysis Three cases shown in Fig. 9(a)–(c) are considered. lx and ly are the length of the plate along X-axis and Y-axis respectively. For the cases, the size of the element is approximately 7.45 mm. The first case (Fig. 9(a)) is a square plate of length 1.2 m, simply supported on all edges with sinusoidal loading as per Eq. (26). The intensity of loading at the centre (q0) is 1 N/mm2.
q = q0 sin
y x sin lx ly
(26)
The deflection (w), bending moments, twisting moment and shear forces at centre and quarter point of the plate obtained using AEM are given in Table 3. The contour plot of the deflected shape of the plate determined using AEM and CPT are shown in Fig. 10(a). Fig. 10(b) 11
Structures 22 (2019) 1–12
D.L. Christy, et al.
shear forces according to CPT are also given in Tables 3, 4 and 5. The deflection, moments and shear forces as per CPT in Table 3 are calculated using Eqs. (e), (f) and (g) in section 27 of [14]. ls in Tables 4 and 5 represents the smaller span. The deflection coefficients and bending moment coefficients of CPT given in Tables 4 and 5 are obtained from Tables 28 and 29 of [14]. Figs. 10–11 and Tables 3–5 shows that AEM could predict reliable values of deflection, stresses, bending moments, twisting moment and shear forces.
8. Conclusions The following results are arrived from the study:▪ Thin plate element for AEM is developed and validated. ▪ The method to incorporate Poisson's ratio in thin plate element is illustrated. ▪ The equations for calculating displacements, strains, stresses, bending moments, twisting moment and shear forces are derived. ▪ Thin plate element in AEM predicted the deflection, normal stress distribution, complementary shear stress distribution, bending moments, twisting moment and shear forces with acceptable accuracy. ▪ The deflection converged as the number of elements and the number of springs increased. The percentage error in deflection reduces to the acceptable limit at a reasonable computational time. ▪ The fundamental natural frequency and mode shape of plates are determined with reasonable accuracy, with the help of the thin plate element developed.
7.2. Convergence study To study the convergence of AEM with respect to the number of elements and number of springs, a simply supported plate with sinusoidal loading (Case 1 of Section 7.1) is considered. The variation of deflection due to the increase in number of elements and springs are shown in Fig. 12(a) and (b) respectively. The percentage difference of the central deflection obtained by AEM from that determined analytically using CPT are plotted against computational time in Fig. 12(c). Fig. 12(a) and (b) shows that as the number of elements and springs increases, the deflection converges towards that of CPT. As usual, convergence is mainly affected by the number of elements rather than the number of springs. From Fig. 12(c), it is seen that the percentage error in deflection reduces to the acceptable limit at a reasonable computational time.
The study mainly focussed on the development of thin plate element in AEM. This can be extended to the analysis of plates having different property and geometry. References [1] Meguro K, Tagel-Din H. A new efficient technique for fracture analysis of structures. Bull of Earthq Resistant Struct Res Cent, IIS, Univ of Tokyo 1997 Mar;30:103–16. [2] Meguro K, Tagel-Din H. A new simplified and efficient technique for fracture behavior analysis of concrete structures. Proceedings of the 3rd international conference on fracture mechanics of concrete and concrete structures (FRAMCOS-3). 1998 Oct.. p. 911–20. [3] Tagel-Din H, Meguro K. Nonlinear simulation of RC structures using applied element method. Struct Eng Earthq Eng 2000 Jul 21;17(2):13–24. [4] Tagel-Din H, Meguro K. Consideration of Poisson’s ratio effect in structural analysis using elements with three degrees of freedom. Bull of Earthq Resistant Cent 1998;31:41–50. [5] Meguro K, Tagel-Din H. Applied element method for structural analysis: theory and application for linear materials. Struct. Eng. Earthq. Eng. 2000;17(1):21–35. [6] Tagel-Din H, Meguro K. Applied element method for dynamic large deformation analysis of structures. Struct. Eng. Earthq. Eng. 2000 Oct 21;17(2):215–24. [7] Meguro K, Tagel-Din H. Applied element method used for large displacement structural analysis. J Nat Disaster Sci 2002 Jun 1;24(1):25–34. [8] Khalil AA. Enhanced modeling of steel structures for progressive collapse analysis using the applied element method. J Perform Constr Facil 2012 Dec;26(6):766–79. https://doi.org/10.1061/(ASCE)CF.1943-5509.0000267. [9] Cismasiu C, Ramos AP, Moldovan ID, Ferreira DF, Filho JB. Applied element method simulation of experimental failure modes in RC shear walls. Comput Concr 2017;19(4). [10] Malomo D, Pinho R, Penna A. Using the applied element method for modelling calcium-silicate brick masonry subjected to in-plane cyclic loading. Earthq Eng Struct Dyn 2018 Jun;47(7):1610–30. [11] Christy DL, Pillai TM, Nagarajan P. Analysis of concrete beams using applied element method. IOP conference series: materials science and engineering. 330(1). 2018 Mar. p. 1–7. [12] Shakeri A, Bargi K. Use of applied element method for structural analysis. KSCE J Civ Eng 2015 Jul 1;19(5):1375–84. https://doi.org/10.1007/s12205-015-0625-4. [13] Prashidha Kharel's Blog [Internet]. Formulating the applied element method: linear 2D part(1) March 4 - [cited 2019 April 18]. Available from:. http://prashidha. blogspot.in/2014/03/formulating-applied-element-method.html; 2014. [14] Timoshenko SP, Woinowsky-Krieger S. Theory of plates and shells. McGraw-Hill; 1959. [15] Ventsel E, Krauthammer T. Thin plates and shells: theory: analysis, and applications. CRC press; 2001.
7.3. Dynamic analysis The dynamic analysis of plates simply supported on all edges and clamped on all edges are studied in this section. Plates of different length to width ratios are analysed by keeping the larger dimension as 1.2 m and element size as approximately 11.88 mm. The fundamental natural frequency (ω) obtained analytically [15] are given by Eqs. (27) and (28) for simply supported and clamped plates respectively. The corresponding values determined using AEM are presented in Table 6.
=
2
1 1 + 2 lx 2 ly
= 4 25
D m
1 1 28 + 4 + 2 2 2 2 lx 4 ly 3 7 lx l y
(27)
D h
(28)
Eh3 2) . 12(1
where, D = Here, ρ and m are the mass density and mass per unit area respectively. Table 6 shows that the natural frequency determined using AEM are in line with that derived by analytical means (based on CPT). AEM predicts the fundamental natural frequency with very less error. The contour plot of the first mode of vibration of simply supported plate, determined using AEM and analytically are shown in Fig. 13(a) and (b) respectively. The corresponding plots of clamped plate are shown in Fig. 14(a) and (b). Fig. 15(a) and (b) shows the surface plot of second mode of simply supported plate and third mode of fixed plate respectively. Figs. 13–15 shows that the thin plate element in AEM is capable of predicting the mode shapes of plates.
12