An analytical solution for bending of transversely isotropic thick rectangular plates with variable thickness

Applied Mathematical Modelling 77 (2020) 1582–1602

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

An analytical solution for bending of transversely isotropic thick rectangular plates with variable thickness Faezeh Yekkalam Tash, Bahram Navayi Neya∗ Babol Noshirvani University of Technology, Babol, Mazandaran, Iran

a r t i c l e

i n f o

Article history: Received 4 April 2019 Revised 14 July 2019 Accepted 19 August 2019 Available online 23 August 2019 Keywords: Variable thickness plates Displacement potential functions Bending exact solution Transversely isotropic plates

a b s t r a c t In this study, the bending solution of simply supported transversely isotropic thick rectangular plates with thickness variations is provided using displacement potential functions. To achieve this purpose, governing partial differential equations in terms of displacements are obtained as the quadratic and fourth order. Then, the governing equations are solved using the separation of variables method satisfying exact boundary conditions. The advantage of the purposed method is that there is no limitation on the thickness of the plate or the way the plate thickness is being varied. No simplifying assumption in the analysis process leads to the applicability and reliability of the present method to plates with any arbitrarily chosen thickness. In order to confirm the accuracy of the proposed solution, the obtained results are compared with existing published analytical works for thin variable thickness and thick constant thickness plate. Also, due to the lack of analytical research on thick plates with variable thickness, the obtained results are verified using the finite element method which shows excellent agreement. The results show that the maximum displacement of the plates with variable thickness is moved from the center toward the thinner plate edge. In addition, results exhibit the profound effects of both thickness and aspect ratio on stress distribution along the thickness of the plate. Results also show that varying thickness has not a profound impact on bending and twisting moments in transversely isotropic plates. Five different materials consist of four transversely isotropic and one isotropic, as a special case, are considered in this paper, which it is shown that the material properties have a more considerable impact on higher thickness plate. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Plates are used in a vast majority of engineering applications as structural elements in naval, aerospace, mechanical, civil engineering and now have found a wider range of potential applications. In most engineering structures, such as bridge slabs, pavements, and decks, plates are under transverse loading, which results in bending. Consequently, the analysis carried out on the bends in these plates forms the basis of planner structural designing and considerable emphasis was given to this field. During the past decades, a large list of literature has been devoted to analyzing thin rectangular plates with the uniform thickness, mostly based on classical plate theory, CPT, that are available in [1]. In addition, the plates with varying thickness acquired immense popularity to fulfill the practical requirements. Such plates are used as an optimizing method in view ∗

Corresponding author. E-mail address: [email protected] (B. Navayi Neya).

https://doi.org/10.1016/j.apm.2019.08.017 0307-904X/© 2019 Elsevier Inc. All rights reserved.

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1583

of improving material utilization or lightening the elements, especially in modern structures and aerospace technology. However, with allowance for the thickness variation in one or two (or both) coordinate directions, the mathematical model and solution of the corresponding boundary-value problem become more complex. As a result, while a considerable amount of information is available on analyzing the uniform thickness plates with various boundary conditions and loading, only a limited number of research on thick rectangular plates with variable thickness have been published. The first elastic solution for a rectangular plate of variable thickness was developed by Olsson which is only valid for a linearly varying thickness [2]. In 1963, Pope [3] suggested an approximate analysis of the small-deflection deformation under normal loading of a thin elastic plate with the symmetrical taper in planform with edges either simply supported or clamped. Gangnon et al. [4] proposed a finite strip technique to solve the thickness variation of thick rectangular plate problems, based on Mindlin’s plate theory to take into account shear deformation. Ohga and Shigmatsu [5] applied the combined boundary element with the transfer matrix method to solve the bending of constant and variable thickness rectangular plates. Fertis [6] developed a convenient and general method to solve variable thickness plates with various boundary conditions and loading by using equivalent flat plates. Since the CPT assumptions are used in their method, this proposed method is solely applicable to thin plates. In the late twentieth century, Kashtalyan and Nemish [7] considered a three-dimensional boundary-value bending of a variable thickness orthotropic plate, with nonplanar faces perpendicular to an applied load, which admitted representation by a double trigonometric series. Chaves et al. [8] used the Boundary Element Method, BEM, formulation based on the Kirchhoff’s hypothesis, derived from dealing with the varying plate thicknesses. In 2003, Zenkour [9] presented a solution to the variable thickness plate problem based on the CPT, and its numerical solution was carried out by using the small parameter method and Levy-type approach. His work is not exact due to using CPT and replacing a variable thickness plate by a set of constant thickness plates. More recently, Xu and Zhou [10] presented the three-dimensional elasticity solution for simply supported isotropic rectangular plates with variable thickness. In their solution, the displacements are assumed as trigonometrical functions and the unknown coefficients are approximately determined by using the expansions of double Fourier sinusoidal series to the upper surface and lower surface conditions of the plate [10]. In 2005, Hughes et al. proposed the idea of using Non-Uniform Rational B-Spline basis functions in the finite element analysis procedure and since then, a bunch of papers published based on isogeometric analysis. For instance, Yu et al. [11] used a combination of the isogeometric analysis, the level set and a simple first order shear deformation theory for simulating free vibration and buckling problems of laminated composite plates with cutouts through a new numerical method. Le-Manh et al. conducted nonlinear bending and buckling analysis of variable thickness composite plates in the framework of isogeometric analysis [12]. Lieu et al. also employed the same methodology in order to analyze bending and free vibration of in-plane bi-directional functionally graded plates with variable thickness; in which through some numerical examples the impacts of variable thickness, material property, thickness ratio, boundary condition on bending and free vibration responses are discussed [13]. When selecting a material for a specific application the material properties must satisfy the function and the operating conditions of the component or the structure being designed. Composite materials wide usage in structural applications interested many attentions in the behavior of anisotropic materials. The biggest advantage of modern composite materials is their highest strength-to-weight ratios in structures. One of the most applicable composite materials in modern technology is transversely isotropic materials. As mentioned, several attempts were made to analyze the variable thickness plates, only a few analytical solutions were derived for isotropic rectangular plates with variable thickness and the number of research for anisotropic plates is even less. Accordingly, in the last decades, a large number of researchers have desired to understand the behavior of transversely isotropic mediums. Although, due to the challenge of mathematical dealing with these materials, they have not comprehended much. Moreover, as structural elements, transversely isotropic plates, which have one axis of symmetry and five constant coefficients, are an important class of structures, widely used in the industrial fields. Both high importance and complication have generated considerable interest in applying a variety of methods for analyzing this kind of plates. The potential function method, PFM, is one of the most efficient methods to obtain the exact solution of elasticity problems. Therefore, several different potential techniques have been developed to solve problems within both displacement and stress formulations. Methods related to the displacement formulation include the scalar and vector potentials from the Helmholtz decomposition, Galerkin vector, and Papkovich–Neuber functions [14]. Displacement potential function method, DPF, is a productive method to solve three-dimensional elasticity problems. Particularly, Love [15] presented an applicable DPF for isotropic materials. Lekhnitskii [16] enhanced Love’s solution for transversely isotropic problems with axial-symmetrically. Next, Lekhnitskii’s potential function improved by Hai-Chang [17] and Nowacki [18] for solving three-dimensional elasticity problems with an arbitrary medium. Ultimately, Lekhnitskii–Hu–Nowackii function has been used to analyze various problems, including full-space and half-space. More recently, the senior author of present work, Navayi Neya [19–22] successfully used the DPF method to present an analytical solution to isotropic and transversely isotropic simply supported rectangular plates with any arbitrary but constant thickness. Also, Moslemi et al. [20,21] used this method to solve the elastic buckling of the same plate. Nateghi Babagi et al. [22] modeled exact dynamic responses of a rectangular plate with an arbitrary constant thickness under moving load employing DPF method proposed for dynamic problems. Based on DPF, Eskandari-Ghadi et al. presented solutions to solve different elastostatic [23] and elastodynamic problems [24–30] either in full space or half space of transversely isotropic materials. For instance, they obtained the dynamic response of an elastic transversely isotropic half-space [31].

1584

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

Fig. 1. Geometry of the considered variable thickness plate.

The purpose of this paper is developing a general analytical method to analyze the responses of simply supported thick transversely isotropic rectangular plates with variable thickness. To achieve this goal, by applying the Lekhnitskii– Hu–Nowacki displacement potential functions, the governing equations are obtained in terms of two partial differential equations. The governing partial differential equations are solved by using separation of variables method satisfying the exact boundary conditions. The main novelty of the presented method is using no simplifying assumption in regard to displacement distribution in the plate and strain or stress distribution in the plate thickness, which leads to the applicability of the results for thin, moderately thick and thick plates. Furthermore, the DPF method, used for variable thickness plates for the first time in this paper, is extremely useful in analyzing different kinds of composite materials. Accuracy and comprehension are the merits of this paper, taking into account the transverse shear, hence the present work results can be used as a benchmark for related studies. 2. Theoretical formulation 2.1. Governing equations Consider a transversely isotropic rectangular plate of length a, width b and variable thickness subjected to an arbitrary load on its surfaces as shown in Fig. 1. The origin of a Cartesian coordinate system is placed on the corner of mid-plane of the plate and the z-axis is along thickness and coincides with the symmetry axis of the material. In addition, the plate is simply supported on four edges and contains linear elastic materials. Material reactions under stresses can be described by a set of constitutive equations which provides a link between stresses and strains in the mediums. The constitutive equations in linear elasticity are represented by the generalized Hooke’s law and can be written as form [14]:

σi j = C εi j i, j = x, y, z

(1)

where σ ij is a second order tensor known as stress tensor and its elements are the stress components. ε ij is another second order tensor known as strain tensor and its elements are the strain components. C is a fourth order tensor known as stiffness tensor. For a transversely isotropic material in rectangular coordinates the elasticity stiffness matrix can be written as [14]:

⎡

A11 ⎢A12 ⎢A13 C=⎢ ⎢0 ⎣ 0 0

A12 A11 A13 0 0 0

A13 A13 A33 0 0 0

0 0 0 A44 0 0

0 0 0 0 A55 0

⎤

0 0 ⎥ 0 ⎥ ⎥ 0 ⎥ ⎦ 0 A66

(2)

in which A11 – A66 represent the elasticity constants of a material. The number of independent elasticity constants in a transversely material element is five as A55 and A66 are dependent to other elastic constants; A55 = A44 and 2A66 = A11 − A12 . The engineering constants are defined as modulus of elasticity, E, shear modulus, G, and Poisson’s ratio,ν , in the isotropy plane, and transverse modulus of elasticity, shear modulus and Poisson’s ratio, as E , G , and ν  , respectively. Elasticity constants are related to engineering constants as follows:



A11 =

E 1−



E E

ν2

(1 + ν ) 1 − ν −



2 EE

ν

2

 , A13 = A12 =

Eν 1 − ν − 2 EE ν  2



F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

A33 =



E  (1 − ν ) 1 − ν − 2 EE ν  2

 , A44 = G , A66 =

E =G 2 (1 + ν )

1585

(3)

By inverting the matrix C, and using Eq. (3), the elasticity constants of the materials can be achieved.

⎡

C −1

1/E ⎢ −ν /E ⎢−ν  /E =⎢ ⎢ 0 ⎣ 0 0

−ν /E 1/E −ν  /E 0 0 0

−ν  /E  −ν  /E  1/E  0 0 0

0 0 0 1/G 0 0

0 0 0 0 1/G 0

⎤

0 0 ⎥ ⎥ 0 ⎥ ⎥ 0 ⎦ 0 2(1 + ν )/E

(4)

Also, The elastic strain-displacement relationship for the material in Cartesian coordinates (x, y, z) can be written in the following form [14]:

εi j =

1 2



∂ ui ∂ u j + ∂j ∂i

(5)

The stress components introduced previously must satisfy the following differential equations of equilibrium in terms of stresses for a body in the absence of body forces as:

⎧ ∂ σx ∂ τxy ∂ τxz ⎪ + + =0 ⎪ ⎪ ∂x ∂y ∂z ⎪ ⎪ ⎨ ∂ τyx ∂ σy ∂ τyz + + =0 ∂y ∂z ⎪ ∂x ⎪ ⎪ ⎪ ⎪ ∂ τzx + ∂ τzy + ∂ σz = 0 ⎩ ∂x ∂y ∂z

(6)

By substituting Eq. (1) into relation (6), replacing strains by displacement components as Eq. (5), and simplifying the outcome, the equation of equilibrium in terms of displacements is represented as follows:



 

  2

 ∂ 2u ∂ 2v ∂ 2w A11 − A12 ∂ 2 u ∂ 2v ∂w ∂ u + A + A + + + A + =0 12 13 44 ∂ x∂ y ∂ x∂ z 2 ∂ x2 ∂ y 2 ∂ x∂ y ∂ z 2 ∂ x∂ z 

    2

 A11 − A12 ∂ 2 u ∂ 2v ∂ 2u ∂ 2v ∂ 2w ∂ v ∂ 2w + + A12 + A11 2 + A13 + A44 + =0 2 ∂ x∂ y ∂ x 2 ∂ x∂ y ∂ y∂ z ∂y ∂ z 2 ∂ y∂ z  2

 

   ∂ 2w ∂ 2u ∂ 2v ∂ 2w ∂ u ∂ 2w ∂ 2v A44 + + A + + A + A + A =0 44 13 13 33 ∂ x∂ z ∂ x 2 ∂ y∂ z ∂ y 2 ∂ x∂ z ∂ y∂ z ∂ z2 A11

(7)

where u, v, and w, are displacement components in the x, y, and z coordinate directions, respectively. Displacement components of any point of a transversely isotropic rectangular plate can be presented in terms of the Lekhnitskii-Hu-Nowackii DPF, F(x, y, z) and χ (x, y, z), as Eq. (8) [23]:

⎧ A13 + A44 ∂ 2 F ∂χ ⎪ ⎪ u=− − ⎪ ⎪ A44 ∂ x∂ z ∂ y ⎪ ⎪ ⎨ A + A44 ∂ 2 F ∂χ v = − 13 + A44 ∂ y∂ z ∂x ⎪ ⎪ 2

⎪ 2 ⎪ A ∂ ∂ 2F ∂ ⎪ 11 ⎪ + F + ⎩w = A44 ∂ x2 ∂ y2 ∂ z2

(8)

By substituting Eq. (8) into Eq. (7), two independent partial differential equations in terms of the Lekhnitskii-Hu-Nowackii DPFs are achieved as follows:

∂2 ∂ 2χ ∂2 + χ + A44 2 = 0 2 2 2 ∂x ∂y ∂x 4

4



2  ∂4 ∂4 ∂4 ∂ ∂ ∂4 A33 A44 F − A13 + 2A13 A44 − A11 A33 + F + A11 A44 + +2 2 2 F =0 ∂ z4 ∂ x2 ∂ z 2 ∂ y2 ∂ z 2 ∂ x4 ∂ y4 ∂x ∂y (A11 − A12 )



(9) (10)

Since the coefficients in the partial differential Eq. (10) are constants, its solution is in the form of exponential function, esz , and thus the characteristic equation is determined as [19]:





A33 A44 s4 − A213 + 2A13 A44 − A11 A33 s2 + A11 A44 = 0

(11)

It is assumed that s21 and s22 are the roots of the Eigen-value Eq. (11) so the two relations (12) and (13) can be written



s21 + s22 = −

A213 + 2A13 A44 − A11 A33 A33 A44



(12)

1586

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

s21 s22 =

A11 A33

(13)

Substituting the Eqs. (12) and (13) into Eq. (10) gives:



∂2 ∂2 ∂2 + + 2 2 2 2 ∂x ∂y s1 ∂ z



∂2 ∂2 ∂2 + + F =0 ∂ x2 ∂ y2 s22 ∂ z2

(14)

Also, the first governing differential equation could be written as:



∂2 ∂2 ∂2 + + χ =0 ∂ x2 ∂ y2 s20 ∂ z2

(15)

where s20 = A66 /A44 . To solve governing Eqs. (14) and (15), let use the separation of variables method, so each of the functions F and χ can be presented as the multiple of three independent functions in terms of x, y, and z as relations (16) and (17).

F (x, y, z ) = f (x )g(y )h(z )

(16)

χ (x, y, z ) = f1 (x )g1 (y )h1 (z )

(17)

2.2. Boundary conditions The differential equations which must be solved, require appropriate boundary conditions. The common types of boundary conditions for elasticity problems, normally formulated by specifying either the displacements or tractions at boundaries [14]. Applying static and kinematic boundary conditions on the edges of a plate leads to find the related unknown coefficient and remaining unknown coefficients can be achieved by applying traction boundary conditions on upper and lower surfaces of the plate. 2.2.1. Edges boundary conditions Static and kinematic boundary conditions of a simply supported rectangular plate are according to relation (18).



x = 0, a y = 0, b

w = 0 at

(18a)

Mx = 0 at x = 0, a

(18b)

My = 0 at y = 0, b

(18c)

where Mx and Mx represent bending moments in the x and y directions, respectively. Substituting Eqs. (16) and (17) in Eqs. (14) and (15), and satisfying boundary conditions of a simply supported plate on all its four edges, the two potential functions F(x, y, z) and χ (x, y, z) will be achieved as Eqs. (19) and (20): ∞  ∞ 

F (x, y, z ) =

sin

 mπ

 nπ

a

b

m=1 n=1

x sin

y [C1 sinh (s1 αmn z ) + C2 cos h(s1 αmn z ) + C3 sinh (s2 αmn z )

+ C4 s2 cosh (s2 αmn z )]

χ (x, y, z ) =

∞  ∞  m=1 n=1

cos

(19)

 mπ

 nπ

a

b

x cos

y [C5 sinh (s0 αmn z ) + C6 cos h(s0 αmn z )]

(20)

2 = (mπ /a )2 + (nπ /b)2 and in which m, and n are the number of cycles of the nth harmonic in the analysis interval. αmn C1 –C6 are coefficients that can be determined by applying the upper and lower surface boundary conditions.

2.2.2. Traction boundary conditions As surfaces of the plate are not perpendicular to Cartesian coordinates system, the stress and strain components must follow a new coordinate system, the transformed one. Let the l1 , m1 , n1 , l2 , m2 , and n2 be the direction-cosines of the main and transformed coordinates, which are the cosines of the angles made by the plate upper and lower faces and the three coordinate axes. According the transformed system, the status of the upper surface of the plate can be written as [14]:



l1 (x, y )σx + m1 (x, y )τyx + n1 (x, y )τzx = q1 (x, y ) m1 (x, y )σy + n1 (x, y )τzy + l1 (x, y )τxy = q2 (x, y ) n1 (x, y )σz + l1 (x, y )τxz + m1 (x, y )τyz = q3 (x, y )

(21)

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1587

The direction cosines in upper surface satisfy

l12 (x, y ) + m21 (x, y ) + n21 (x, y ) = 1

(22)

And for lower surface of the plate as



l2 (x, y )σx + m2 (x, y )τyx + n2 (x, y )τzx = q4 (x, y ) m2 (x, y )σy + n2 (x, y )τzy + l2 (x, y )τxy = q5 (x, y ) n2 (x, y )σz + l2 (x, y )τxz + m2 (x, y )τyz = q6 (x, y )

(23)

where

l22 (x, y ) + m22 (x, y ) + n22 (x, y ) = 1

(24)

As mentioned earlier, the plate is subjected to an arbitrary load on its surfaces. The expression of the applied load, q(x, y), has to be sought in the form of an infinite Fourier series as follows [1]:

q(x, y ) =

∞  ∞ 

qmn sin

m=1 n=1

mπ x nπ y sin , (m, n = 1, 3, 5, . . .) a b

(25)

where qmn represents coefficient to be determined with regard to the way the load is applied; in this study, qmn = 16/mnπ 2 . In Eqs. (21) and (23), qi (x,y) is the component of external loading, so there is not any limitation in exposed load, and it could be perpendicular or not, so the presented method is applicable for a large group of plate problems. By replacing displacement components in relation (5) and putting the obtained elastic strains in the relation (1) the stress relations regarding the stress-displacement and displacement potential functions could be written as:

     3

 A + A44 ∂ 3 F A13 + A44 ∂ 3 F ∂ 2χ ∂ 2χ A11 ∂ 3F ∂ 3F ∂ F σx = A11 − 13 − + A − + + A + + 12 13 A44 A44 A44 ∂ x2 ∂ z ∂ x 2 ∂ z ∂ x∂ y ∂ y 2 ∂ z ∂ x∂ y ∂ y2 ∂ z ∂ z3     

 A + A44 ∂ 3 F A13 + A44 ∂ 3 F ∂ 2χ ∂ 2χ A11 ∂ 3F ∂ 3F ∂ 3F − + A − + + A + + σy = A12 − 13 11 13 A44 A44 A44 ∂ x2 ∂ z ∂ x 2 ∂ z ∂ x∂ y ∂ y 2 ∂ z ∂ x∂ y ∂ y2 ∂ z ∂ z3     

 A13 + A44 ∂ 3 F A13 + A44 ∂ 3 F ∂ 2χ ∂ 2χ A11 ∂ 3F ∂ 3F ∂ 3F − + A13 − + + A33 + + σz = A13 − A44 A44 A44 ∂ x2 ∂ z ∂ x 2 ∂ z ∂ x∂ y ∂ y 2 ∂ z ∂ x∂ y ∂ y2 ∂ z ∂ z3 

 A + A44 ∂ 3 F ∂ 2χ ∂ 2χ + − τxy = A66 −2 13 A44 ∂ x∂ y∂ z ∂ x2 ∂ y2    3

 3 2 A13 + A44 ∂ F ∂ χ A11 ∂ 3F ∂ 3F ∂ F + + + + τyz = A44 − A44 A44 ∂ x2 ∂ y ∂ y ∂ z 2 ∂ x∂ z ∂ y3 ∂ y∂ z 2   

 A13 + A44 ∂ 3 F ∂ 2χ A11 ∂ 3 F ∂ 3F ∂ 3F − + + + (26) τxz = A44 − A44 A44 ∂ x3 ∂ x∂ z 2 ∂ y ∂ z ∂ x∂ y 2 ∂ x∂ z 2

By replacing the relations (26) into Eqs. (21) and (23), and solving the equations by using MATLAB software, the remaining coefficients could be obtained. 3. Numerical calculation and discussion In this study, the variation of plate thickness is arbitrary which can be occurred in x or y direction, even both. Thus z1 (x,y) and z2 (x,y) are considered as the functions describing the variation of upper and lower surfaces of the plate, respectively.



 y

z1 (x, y ) = −h0 1 + α1 2



b

 y

z2 (x, y ) = h0 1 + α2 2

b

i 

−1

k 

−1

 x

1 + β1 2

a

 x

1 + β2 2

a

j

−1

l 

−1

(27)

which i, j, k, and l are the degree of non-uniformity in x and y directions, α 1 , α 2 , β 1 , and β 2 are the slope of the thickness variation of the plate, and h0 is the half of the plate central thickness. It can be derived from relation (27) that the plate can have any shape, symmetric or antisymmetric. The study results are provided for five different materials, including steel as isotropic, magnesium and zinc, as metallic substances, graphite-epoxy and E-glass-epoxy, as composite substances, having the status of transversely isotropic materials. The elastic constants of these materials are given in Table 1. Also, the elasticity constants of materials are achieved using relation (4) as Table 2. For brevity of notation, non-dimensional parameters η as thickness ratio (h/a), δ as aspect ratio (b/a), ξ as relative position in the y direction (y/b), and γ as relative position in the x direction (x/a), are defined. In addition, the proposed thickness ratio for thin, moderately thick and thick plates are considered as 0.01, 0.20 and 0.50, respectively. In this context, the deflections, moments, stresses, and forces are non-dimensionalized using the relations presented in

1588

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

Fig. 2. Geometry of the rectangular plate with linear thickness variation in y direction.

Fig. 3. Comparison of the PW with the FEM results for a thick square E-glass epoxy plate with linearly varying thickness in the y direction.

Table 1 Material properties (GPa) [12]. Material

A11

A12

A33

A44

A55

Magnesium Zinc Graphite-epoxy E Glass-epoxy Steel

59.700 160.900 8.280 14.930 269.230

61.700 61.000 86.800 47.270 269.230

21.700 50.070 0.285 5.244 115.380

26.200 33.500 2.767 6.567 115.380

16.400 38.300 4.147 4.745 76.920

Table 2 Elasticity constants of the materials. Material

Ex (GPa)

Ez (GPa)

ν xy

ν xz

Gxy (GPa)

Magnesium Zinc Graphite-epoxy E Glass-epoxy Steel

45.4540 11.9330 73.5290 11.8480 200.0000

50.7350 35.2110 82.9560 44.6420 200.0000

0.3568 0.2577 0.3338 0.4170 0.3000

0.2526 0.0632 0.0260 0.2455 0.3000

14.3260 38.2990 4.1470 4.7460 76.9230

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1589

Fig. 4. Non-dimensional central mid-plane deflection (w∗ ) of a thin square plate versus relative position (ξ ).

Fig. 5. Non-dimensional central mid-plane deflection (w∗ ) of a moderately thick square plate versus relative position (ξ ).

dimensionless forms as (28):

100wD0 q0 a4 10 Mx Mx∗ = q0 a2 10My My∗ = q0 a2 10h2 σx σx∗ = q0 10 h 2 σy σy∗ = q0 Q y Qy∗ = q0 w∗ =

10Mxy q0 a2 10Mxy = q0 a2

∗ Mxy = ∗ Mxy

τxy∗ = τxz∗ = τyz∗ = ∗ Qxy

τxy q0

τxz q0

τyz

q0 Qxy = q0

where D0 = Eh3 /12(1 − ν 2 ) and h is the plate central thickness.

(28)

1590

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

Fig. 6. Non-dimensional central mid-plane deflection (w∗ ) of a thick square plate versus relative position (ξ ).

Fig. 7. Non-dimensional central mid-plane deflection (w∗ ) of a square plate versus thickness ratio (η). Table 3 Comparison of non-dimensional results of isotropic square plates with constant thickness under uniform distributed load.

Present work

Nematzadeh et al. [19]

η

w∗

Mx∗

My∗

∗ Mxy

Qx∗

Qy∗

0.01 0.10 0.20 0.01 0.10 0.20

0.4064 0.4248 0.4803 0.4064 0.4249 0.4804

0.4790 0.4805 0.4851 0.4790 0.4805 0.4851

0.4790 0.4805 0.4851 0.4790 0.4805 0.4851

0.3244 0.3244 0.3244 0.3244 0.3244 0.3244

0.3325 0.3325 0.3325 0.3325 0.3325 0.3325

0.3325 0.3325 0.3325 0.3325 0.3325 0.3325

3.1. Verifications The plate studied in this paper has major features, including arbitrary variable thickness and transversely isotropic materials. The present work results, for a special case, when the thickness of the plate is constant but arbitrary, are compared with results presented by Nematzadeh et al. [19]. Accordingly, the non-dimensional displacements and internal moments of a simply supported square plate with different thickness ratios have been tabulated in Table 3. As could be seen from Table 3, the results have a great agreement with the results published by Nematzadeh et al. On the other hand, to verify

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1591

Fig. 8. Non-dimensional central bending moments (My∗ ) of a square plate versus thickness ratio (η).

∗ Fig. 9. Non-dimensional corner twisting moments (Mxy ) of a square plate versus thickness ratio (η) for γ = ξ = 1.0.

the employed method for the variable thickness plates (Fig. 2), the obtained results in this paper and result presented by Zenkour [9] tabulated in Table 4. For the sake of comparison, all contents presented for a linearly varying lower surface in y direction by taper ratios as α 1 = 0.0 and α 2 = 0.2, where the plate substance is isotropic and the engineering constants are chosen as E = 2.1 × 1011 N/m2 and ν = 0.3. Regarding Table 4, the present method is applicable for isotropic plates as an especial case of the transversely isotropic materials. Although, the minor distinction between the present work and Zenkour results is due to simplifications used by Zenkour. Moreover, due to the lack of analytical results for thick varying thickness plate, based on the best knowledge of authors, a finite element method, FEM, using ABAQUS software has been carried out to validate the present work, PW. Evidenced by a square plate (a = b = 1m) with symmetrical linear variation in the y direction. The plate is made of E-glass epoxy substance and subjected to a uniform load on its upper surface. The deflection distributions through the width direction are plotted and presented in Fig. 3. It follows that there is a very good agreement between analytical solution (PW) and the FEM, however some difference between exact and numerical works is inevitable. 3.2. Results and discussion This section attempts to demonstrate how various geometrical parameters such as thickness ratio, aspect ratio, thickness taper ratio, and material properties can have an impact on the behavior of transversely isotropic rectangular plates with variable thicknesses. To achieve this, a linearly variable thickness plate (see Fig. 2) with the constant length a = 1.0m and in three thickness ratios namely thin, moderately thick and thick plates is considered. The plate is under a uniform distributed

1592

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

Fig. 10. Variation of central normal stress (σy∗ ) for different thickness ratio in a square E-glass epoxy plate and γ = ξ = 0.5.

∗ Fig. 11. Distribution of non-dimensional shear stress (τxy ) along the thickness of a square E-glass epoxy plate.

Table 4 Comparison of non-dimensional results of isotropic square plates with varying thickness under uniform distributed load.

δ

ξ

Zenkour [9] w

1.0

1.5

2.0

0.25 0.50 0.75 0.25 0.50 0.75 0.25 0.50 0.75

∗

0.3391 0.4100 0.2625 0.6904 0.7795 0.4914 0.9759 1.0229 0.6453

Present work Mx∗

My∗

w∗

Mx∗

My∗

0.4186 0.4772 0.2727 0.7363 0.8115 0.4666 0.9745 1.0194 0.5908

0.5160 0.4628 0.2130 0.6033 0.4774 0.2317 0.6154 0.4427 0.2405

0.3458 0.4100 0.2566 0.6731 0.7775 0.4993 0.9154 1.0186 0.6789

0.4175 0.4772 0.2701 0.7359 0.8253 0.4670 0.9585 1.0214 0.5910

0.5146 0.4635 0.2100 0.6143 0.4567 0.2317 0.6234 0.4342 0.2315

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1593

Fig. 12. Variation of non-dimensional lateral displacement (w∗ ) for different aspect ratio in a thick E-glass epoxy plate on the line γ = 0.5.

Fig. 13. Variation of non-dimensional bending moment (My∗ ) for different aspect ratio in a thick E-glass epoxy plate on the line γ = 0.5.

load on its upper surface. As a result of this study, the efficiency and effectiveness of the current method in solving the bending problems of a thick plate that vary in the thickness are shown. 3.2.1. The effect of the plate thickness In order to investigate the effect of the thickness ratio on the static responses of a square plate, Figs. 4–6 represent the normalized deflection distribution along the y direction. Regarding Fig. 4, there are minor distinctions in non-dimensional deflections of the thin plates with different substances, these differences become more noticeable as the plate thickness ratio increases, so these differences in the range of the thick plates are the highest. Figs. 5 and 6, exhibit the rises in the distinction between non-dimensional deflections as the plate thickness rising. The magnitude of E-glass epoxy deflections is more compared to other studied materials. Fig. 7 shows the non-dimensional central deflections at the mid-plane of a variable thickness plate. As seen, increasing the thickness ratio raises the value of deflections in all substances, apart from zinc. This deflection rise is more prominent in the range of moderately thick and thick plates in comparison with thin plates. As mentioned, this trend is different about zinc, where the deflection increases by increasing the thickness ratio in the range of 0 < η  0.4, and it rapidly reaches a peak at η = 0.4. Then, in the thickness ratio η > 0.4, the curve turns downward and declines slightly. Figs. 8 and 9 show the distribution of the bending and twisting moments in different materials by changing the plate thickness, respectively. As it is evident from Fig. 8, the bending moments of all materials, except zinc, have significantly

1594

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

∗ Fig. 14. Variation of non-dimensional twisting moment (Mxy ) for different aspect ratio in a thick E-glass epoxy plate on the line γ = 1.0.

Fig. 15. Distribution of non-dimensional normal stress (σy∗ ) for different aspect ratio in a thick E-glass epoxy plate at the section γ = 0.5, z = z1 .

low variations. Two composite materials, graphite epoxy, and E-glass epoxy, have the lowest variations along the changing thickness ratio. The variation of the central bending moment of graphite epoxy remained 0.492 approximately. In the case of E-glass epoxy, this amount differed from 0.5215 to 0.5521, in range of 0 < η  1. For two metallic substances, steel and magnesium, there is an upward trend by increasing the thickness ratio. The only exception is zinc, as mentioned, whose central bending moment varied from 0.3475 to 0.8327 by increasing the thickness ratio, which is about 130% increase. Fig. 9 presents the non-dimensional corner twisting moments, demonstrating that the trends are analogous to the central bending moment, described in advance. Moreover, by increasing the thickness ratio, the corner twisting moments in zinc, magnesium, and steel fall and then vanish and reverse signs at higher thickness ratios; thus Kirchhoff’s supplementary forces at the corners decrease and their directions are then reversed. This trend in zinc, as the only material giving complex characteristic roots, is much more plain than the other materials. Based on Fig. 9, the twisting moment decreases as the thickness ratio grows in the range of 0  η  1, but in different changing rate. Also, in all examined materials except composite ones, there is a point reaching zero and gain reverse sign. By way of illustration, Table 5 presents the taper rates of a steel plate response (TR) by increasing the thickness. The data in Table 5 indicates that by increasing the thickness ratio from 0.01 to 0.1 which is 10 times more, the deflection and bending moment of the plate increase and the twisting moment decreases, slightly. On the other hand, increasing the thickness ratio from 0.2 to 0.5 (150% increase) leads to over a 72% rise in the deflection of the plate which is quite significant. Also, the bending and twisting moment change is not noticeable.

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1595

Fig. 16. Distribution of non-dimensional normal stress (σy∗ ) for different aspect ratio in a thick E-glass epoxy plate at the section γ = 0.5, z = z2 .

∗ ) for different aspect ratio in a thick E-glass epoxy plate at the section γ = 0.0, z = z1 . Fig. 17. Distribution of non-dimensional shearing stress (τxy

Table 5 Taper rate (TR) of the non-dimensional central results for steel by changing the plate thickness.

η

TR

w∗

0.01 0.1 0.2 0.5

900%

0. 0. 0. 0.

150%

4405 4589 5143 8852

TR

My∗

TR

∗ Mxy

4.17%

0.04784 0.04798 0.04839 0.05142

0.29%

0.03235 0.03215 0.03154 0.02760

72.11%

6.26%

TR −0.62% −12.49%

Distribution of the non-dimensional normal stress (σy∗ ) through the thickness of a thick square E-glass epoxy plate is presented in Fig. 10. The graph is provided along the plate centerline, where γ = 0.5 and ξ = 0.5. As indicated in Fig. 10, the distribution of normal stress remains linear, as long as the thickness of the plate is in the range of thin and moderately thick plates (η ≤ 0.2). By increasing the thickness of the plate and entering the range of the thick plates, this non-linearity becomes more and more. Fig. 11 illustrates how the thickness ratio has an impact on the distribution of shear stress along the z direction when the thickness varies from 0.01 to 1.0. Regarding this figure, an increase in plate thickness leads to a sharp decrease in the maximum value of the non-dimensional shear stress. Also, the maximum shear stress in thin plates, as in the case of a rectangular cross-section of a beam, occurs at z = 0.0, but by increasing the thickness of the plate, it moves from

1596

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

Fig. 18. Variation of non-dimensional transverse shear force (Qy∗ ) for different aspect ratio in a thick square E-glass epoxy plate on the line γ = 0.0.

Fig. 19. Non-dimensional mid-plane deflection (w∗ ) for varying the thickness taper ratio in a thick square E-glass epoxy plate.

the mid-plane to upper surface of the plate, which is loaded. It should be noted that the total area under each curve is invariant. As the transverse shear force of the plate is the integral (the total area below the graph) of shear stress, it is obvious why the transverse shear force is independent of the plate thickness [19]. 3.2.2. The effect of the aspect ratio of a plate The aspect ratio can be another factor that impacts on the responses of a plate structure. Hence, the effect of a thick plate aspect ratio with unchanged varying thickness as the rate of α 1 = α 2 = 0.2 is studied in this section. The plate aspect ratios (δ ) have been chosen 0.4, 1.0, and 2.5. In Figs. 12–14, the deflection, bending moment, and twisting moment of a thick plate are plotted against the relative position in the y direction (ξ ), for E-glass epoxy as an example. According to Fig. 12, by increasing the aspect ratio of the plate, the non-dimensional displacement decreases. It conforms to the fact that at lower aspect ratios, the end effects become dominant as restrict the deflection of the plate; here the displacements are zero at both ends and the maximum values of displacement (w∗ ) are increasing as the aspect ratio increases. From the distribution of ∗ ) versus to ξ (Figs. 13 and 14), can be concluded that non-dimensional bending and twisting moment components (My∗ , Mxy at higher aspect ratios, the graphs become irregular and sinusoidal. Furthermore, by comparing the graphs with δ = 1.0, and 2.5, it is obvious that these graphs are closer compared to the graphs with δ = 0.4, and 1.0. Overall, concerning Figs. 12–14, the plate aspect ratio profoundly influences the plate behavior and a large aspect ratio makes the plate analysis inefficient. The effect of the aspect ratio on the distribution of the normal stress (σy∗ ) along the y direction of a plate is the focus of Figs. 15 and 16. These figures are plotted on two different surfaces of the plate (z = z1 , z = z2 ) for various aspect ratios.

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1597

Fig. 20. Variation of non-dimensional mid-plane deflection (w∗ )for varying the thickness (α 1 = α 2 ) in two section of a thick square plate.

Fig. 21. Variation of non-dimensional bending moment (My∗ ) by varying the thickness taper ratio in a thick square E-glass epoxy plate.

These graphs show that the aspect ratio has a remarkable effect on the distribution of the normal stress (σy∗ ) in both upper and lower surfaces of the plate, and the stresses increase by increasing the aspect ratio, so the normal stresses and aspect ratio have a positive correlation. Also, by increasing the aspect ratio on both the upper surface and lower surface of the plate, the maximum normal stress shifts to the thinner edge of the plate. Fig. 17 shows the non-linear shearing stress distribution of a thick plate on the corner of the top surface of a plate along the y direction (x = 0.0). Based on this figure, the shear stress increases with the increasing aspect ratio (δ ), which reveals the significant effect of aspect ratio. Regarding Fig. 18 the transverse shear force (Qy∗ ) increases by increasing values of aspect ratio from 0.4 to 1.0, whereas it decreases with the increasing value of the aspect ratio from 1.0 to 2.5, also this graph fluctuates slightly.

3.2.3. The effect of variation in the plate thickness The non-dimensional deflection versus y direction of a square plate containing a transversely isotropic material (E-glass epoxy) has been shown in Fig. 19. According to Eq. (25), the slope of thickness variation rise leads to the fall in the plate thickness in the first half and increase in the latter half; the middle thickness remains constant. Fig. 19 shows when the thickness is constant (α 1 = α 2 = 0.0), the maximum deflection occurs at the plate centerline (ξ = 0.5) and by increasing the thickness variation rate, the maximum deflection moves from the center to the plate’s thinner edge.

1598

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

∗ Fig. 22. Variation of non-dimensional twisting moment (Mxy ) by varying the thickness taper ratio in a thick square E-glass epoxy plate.

Fig. 23. Variation of non-dimensional normal stress (σy∗ ) by different taper ratio in thickness for a thick square E-glass epoxy plate on the line γ = 0.5, z = z1 .

As seen earlier on Figs. 4–6, the higher the thickness ratio is, the more effective material properties are. Thus in order to study the transversely isotropic plates, only the figures related to the thick plates are presented. In addition, Fig. 20 shows the effect of the slope of thickness variation on the non-dimensional deflection in two sections of the thick plate, where ξ = 0.25 and ξ = 0.75. It can be observed that the influence of changing the slope of the plate thickness variation on the deflection in section ξ = 0.25 is more dominant in comparison with section ξ = 0.75. Evidently, the effect of the material properties is insignificant and the trend is the same in all the examined materials here. Figs. 21 and 22 represent the variation of the bending and twisting moments versus the y direction for a thick square plate with the different slopes of thickness variation, respectively. From Fig. 21, it can be concluded that the bending moment of a transversely isotropic thick plate is insensitive to the variation of thickness. In addition, Fig. 22 shows that the magnitude of the twisting moment in ξ = 0.0 increases and in ξ = 1.0 decreases by increasing the slope of thickness variation. Table 6 presents the position of maximum deflection of a thick square plate in terms of different variation slope and different materials. According to the contents of this table, when the thickness of the plate is constant, the maximum deflection occurs precisely in the middle of the plate. However, the maximum displacement of a plate with variable thickness moves from the center toward the thinner plate edge. Fig. 23 represents the effect of the slope of thickness variation on the distribution of normal stress (σy∗ ) along the y direction. This figure is provided for a transversely isotropic material (E-glass epoxy), on the upper surface of a thick square

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1599

∗ Fig. 24. Variation of non-dimensional shear stress (τxy ) by different taper ratio in thickness for a thick square E-glass epoxy plate on the line γ = 0.0, z = z1 .

∗ ) by different taper ratio in thickness for a thick square E-glass epoxy plate on the line γ = 0.0, Fig. 25. Variation of non-dimensional shear stress (τxy z = z2 .

Table 6 The relative position of maximum deflection of a thick plate in terms of different variation ratios and different materials in the y direction.

α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5

Steel

Magnesium

Zinc

Graphic epoxy

E-glass epoxy

0.5000 0.4461 0.3921 0.3376 0.2823 0.2258

0.5000 0.4297 0.3963 0.3435 0.2892 0.2328

0.5000 0.4297 0.3647 0.3067 0.2539 0.2037

0.5000 0.4594 0.3959 0.3423 0.2870 0.2297

0.5000 0.4614 0.4011 0.3499 0.2964 0.2398

plate. It can be concluded that in a uniform thickness plate, the variation of σy∗ is symmetric with respect to the middle of the plate. By increasing the slope of variation, the magnitude increase of normal stress (σy∗ ) in the first half is more remarkable than the magnitude decrease of normal stress (σy∗ ) of the second half of the plate. An analysis of the results indicates that the deflection of thin, moderately thick and thick plates has less sensitivity to the amount of m, and n devoted in the analysis process, compared to stress or strain distributions. This sensitivity in stress

1600

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

Table 7 Maximum deflection w∗max of a thick plate in terms of different variation ratios and different materials.

Steel

Magnesium

Zinc

Graphic Epoxy

E-Glass epoxy

α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5

β 1 = β 2 = 0.0

β 1 = β 2 = 0.1

β 1 = β 2 = 0.2

β 1 = β 2 = 0.3

β 1 = β 2 = 0.4

β 1 = β 2 = 0.5

0.6988 0.7033 0.7175 0.7426 0.7816 0.8402 0.7876 0.7921 0.8063 0.8315 0.8706 0.9291 0.4680 0.4741 0.4927 0.5250 0.5736 0.6437 0.7630 0.7671 0.7800 0.8031 0.8392 0.8939 0.8947 0.8990 0.9125 0.9365 0.9738 1.0298

0.7033 0.7173 0.7416 0.7785 0.8323 0.9108 0.7921 0.8061 0.8305 0.8676 0.9215 0.9999 0.4741 0.4924 0.5234 0.5690 0.6329 0.7229 0.7672 0.7799 0.8022 0.8364 0.8866 0.9607 0.8990 0.9123 0.9356 0.9711 1.0228 1.0982

0.7175 0.7416 0.7775 0.8286 0.9009 1.0058 0.8063 0.8305 0.8666 0.9179 0.9904 1.0950 0.4927 0.5234 0.5676 0.6280 0.7104 0.8262 0.7800 0.8022 0.8355 0.8833 0.9516 1.0515 0.9125 0.9356 0.9702 1.0195 1.0894 1.1906

0.7426 0.7785 0.8286 0.8979 0.9950 1.1361 0.8315 0.8676 0.9179 0.9874 1.0846 1.2253 0.5250 0.5690 0.6280 0.7066 0.8133 0.9649 0.8031 0.8364 0.8833 0.9487 1.0414 1.1772 0.9365 0.9711 1.0195 1.0866 1.1808 1.3178

0.7816 0.8322 0.9009 0.9950 1.1272 1.3213 0.8706 0.9214 0.9904 1.0846 1.2166 1.4102 0.5736 0.6329 0.7104 0.8133 0.9546 1.1589 0.8392 0.8866 0.9516 1.0414 1.1687 1.3573 0.9738 1.0227 1.0894 1.1808 1.3096 1.4989

0.8402 0.9108 1.0058 1.1361 1.3213 1.5981 0.9291 0.9999 1.0950 1.2253 1.4102 1.6863 0.6437 0.7229 0.8262 0.9649 1.1589 1.4454 0.8939 0.9607 1.0515 1.1772 1.3573 1.6285 1.0298 1.0982 1.1906 1.3178 1.4989 1.7704

Table 8 Maximum relative deflection w∗max /w∗c of a thick plate in terms of different variation ratios and different materials.

Steel

Magnesium

Zinc

Graphic epoxy

‘E-Glass Epoxy

α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5 α 1 = α 2 = 0.0 α 1 = α 2 = 0.1 α 1 = α 2 = 0.2 α 1 = α 2 = 0.3 α 1 = α 2 = 0.4 α 1 = α 2 = 0.5

β 1 = β 2 = 0.0

β 1 = β 2 = 0.1

β 1 = β 2 = 0.2

β 1 = β 2 = 0.3

β 1 = β 2 = 0.4

β 1 = β 2 = 0.5

1.0000 1.0065 1.0267 1.0627 1.1185 1.2024 1.0000 1.0058 1.0238 1.0558 1.1055 1.1798 1.0000 1.0130 1.0529 1.1219 1.2257 1.3755 1.0000 1.0054 1.0223 1.0525 1.0999 1.1716 1.0000 1.0049 1.0199 1.0467 1.0885 1.1511

1.0065 1.0264 1.0612 1.1140 1.1910 1.3034 1.0058 1.0236 1.0546 1.1016 1.1700 1.2696 1.0130 1.0521 1.1185 1.2159 1.3524 1.5446 1.0054 1.0221 1.0513 1.0962 1.1620 1.2592 1.0049 1.0197 1.0457 1.0854 1.1432 1.2275

1.0267 1.0612 1.1126 1.1858 1.2893 1.4393 1.0238 1.0546 1.1004 1.1656 1.2576 1.3904 1.0529 1.1185 1.2128 1.3419 1.5181 1.7653 1.0223 1.0513 1.0950 1.1577 1.2472 1.3781 1.0199 1.0457 1.0844 1.1395 1.2176 1.3308

1.0627 1.1140 1.1858 1.2849 1.4239 1.6259 1.0558 1.1016 1.1656 1.2538 1.3772 1.5559 1.1219 1.2159 1.3419 1.5098 1.7379 2.0617 1.0525 1.0962 1.1577 1.2434 1.3648 1.5429 1.0467 1.0854 1.1395 1.2146 1.3199 1.4729

1.1185 1.1910 1.2893 1.4239 1.6130 1.8909 1.1055 1.1700 1.2576 1.3772 1.5449 1.7907 1.2257 1.3524 1.5181 1.7379 2.0398 2.4762 1.0998 1.1620 1.2472 1.3648 1.5317 1.7789 1.0885 1.1432 1.2176 1.3199 1.4637 1.6754

1.2024 1.3034 1.4393 1.6259 1.8909 2.2870 1.1798 1.2696 1.3904 1.5559 1.7907 2.1412 1.3755 1.5446 1.7653 2.0617 2.4762 3.0885 1.1716 1.2592 1.3781 1.5429 1.7789 2.1344 1.1511 1.2275 1.3308 1.4729 1.6754 1.9789

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

1601

distribution increases by increasing the plate thickness; in the thick plates is the most. Evidently, the oscillatory behavior in a thick plate, Figs. 15–17, and 23, is due to the numbers allotted to m and n, which is not an infinite amount. ∗ ) along the width of a square thick E-glass Moreover, Figs. 24 and 25 provide information on shear stress distribution (τxy epoxy plate and on the upper and lower surfaces, respectively. According to these Figs. an increase in the taper ratio leads to a decrease in the thickness in the first half of the plate, so the maximum value of shear stress increases noticeably. However, in the second half, the maximum shear stress declines slightly. Tables 7 and 8 present the results of the current method for a square plate with linearly varying thickness subjected to a uniform distributed load for five different materials considered herein. Table 7 demonstrates how the variation of the thickness impacts on the deflection of the plate. Based on this Table, increasing the slope of thickness variation in any direction leads to an increase in maximum deflection of the plate. The maximum deflection, regardless of material type, when the thickness varies in both x, and y directions occurs in the highest taper ratio. In order to detect the trend of the deflection changing, Table 8 has been provided. Deflection of a constant thickness plate termed as w∗c , and all the table contents represent as w∗max /w∗c proportions. According to Table 8, the intensification trend towards the increasing the slope of thickness variation in zinc is remarkable compared to other substances. 4. Conclusion In this paper, the bending behavior of transversely isotropic rectangular plates with varying thickness have been studied using the displacement potential functions. Governing differential equations solved by applying the separation of variables method and satisfying exact plate boundary conditions. Due to no simplifying assumptions, the present work results are accurate and applicable to any plate with arbitrary thickness ratio. The validity of the proposed method is proven by comparing the obtained results with the available published papers and with the numerical method which an excellent agreement was achieved. The most significant feature of the present analysis is that the solving procedure is independent of the thickness of the plate. The main conclusion to be drawn here is that by increasing the variation rate of thickness, the maximum displacement position of the plate shifts toward the thinner edge. Moreover, the results show that the distribution of the twisting moments by changing the thickness ratio in different transversely isotropic materials is unequal. Two composite materials, graphite epoxy and E-glass epoxy, have the lowest variations along the changing thickness ratio. However, the corner twisting moments in zinc, magnesium, and steel fall and then vanish and reverse signs at higher thickness ratios; hence Kirchhoff’s supplementary forces at the corners decline and then their directions are reversed. This trend in zinc, as the only material giving complex characteristic roots, is much more clear than the other studied materials. References [1] R. Szilard, Theories and Applications of Plate Analysis: Classical, Numerical and Engineering Methods, John Wiley & Sons, Inc, 2004. [2] R.G. Olsson, Biegung der rechteckplatte bei linear varanderlicher biegung steifigkeit, Arch. Appl. Mech. 5 (1934) 363–373. [3] G.G. Pope, The Bending Under Normal Loading of Plates Tapered in planform, Aeronautical Research Council Report and Memorandam, Minstry of Aviation, London, 1963, p. 3325. [4] P. Gagnon, C. Gosselin, L. Cloutier, A finite strip method for the analysis of variable thickness rectangular thick plates, Comput. Struct. 63 (1997) 349–362. [5] M. Ohga, T. Shigematsu, Bending analysis of plates with variable thickness by boundary element-transfer matrix method, Comput. Struct. 28 (1988) 635–640. [6] D.G. Fertis, Nonlinear Mechanics, CRC Press, 1998. [7] M.Y. Kashtalyan, Y.N. Nemish, Three-dimensional bending stress-strain state of rectangular variable-thickness orthotropic plates, Int. Appl. Mech. 30 (1994) 952–961. [8] E.W.V. Chaves, G.R. Fernandes, W.S. Venturini, Plate bending boundary element formulation considering variable thickness, Eng. Anal. Bound. Elem. 23 (1999) 405–418. [9] A.M. Zenkour, An exact solution for the bending of thin rectangular plates with uniform, linear, and quadratic thickness variations, Int. J. Mech. Sci. 45 (2003) 295–315. [10] Y.P. Xu, D. Zhou, Three-dimensional elasticity solution of transversely isotropic rectangular plates with variable thickness, Iran. J. Sci. Technol. 34 (2010) 353–369. [11] T. Yu, S. Yin, T. Quoc-Bui, S. Xia, S. Tanaka, S. Hirose, NURBS-based isogeometric analysis of buckling and free vibration problems for laminated composites plates with complicated cutouts using a new simple FSDT theory and level set method, Thin Wall. Struct. 101 (2016) 141–156. [12] T. Le-Manh, Q. Huynh-Van, T.D. Phan, H.D. Phand, H. Nguyen-Xuan, Isogeometric nonlinear bending and buckling analysis of variable-thickness composite plate structures, Compos. Struct. 159 (2017) 818–826. [13] Q.X. Lieu, S. Lee, J. Kang, J. Lee, Bending and free vibration analyses of in-plane bi-directional functionally graded plates with variable thickness using isogeometric analysis, Compos. Struct. 192 (2018) 434–451. [14] M.H. Sadd, Elasticity: Theory, Applications, and Numerics, Academic Press, 2009. [15] A.E.H. Love, A Treatise on the Mathematical Theory of Elasticity, Cambridge University Press, 1892. [16] S.G. Lekhnitskii, Theory of Elasticity of an Anisotropic Body, Holden-Day, 1963. [17] H. Hai-Chang, On the three-dimensional problems of the theory of elasticity of a transversely isotropic body, Acta Phys. Sin. 9 (1953) 130–148. [18] W. Nowacki, The stress function in three-dimensional problems concerning an elastic body characterized by transverse isotropy, Appl. Mech. 2 (1954) 21–25. [19] M. Nematzadeh, M. Eskandari-Ghadi, B. Navayi-Neya, An analytical solution for transversely isotropic simply supported thick rectangular plates using displacement potential functions, J. Strain Anal. Eng. 46 (2011) 121–142. [20] A. Moslemi, B. Navayi-Neya, J. Vaseghi-Amiri, 3-D elasticity buckling solution for simply supported thick rectangular plates using displacement potential functions, Appl. Math. Model. 40 (2016) 5717–5730. [21] A. Moslemi, B. Navayi-Neya, J. Vaseghi-Amiri, Benchmark solution for buckling of thick rectangular transversely isotropic plates under biaxial load, Int. J. Mech. Sci. 131-132 (2017) 356–367. [22] P. Nategh-Babagi, B. Navayi-Neya, M. Dehestani, Three dimensional solution of thick rectangular simply supported plates under a moving load, Meccanica 52 (2017) 3675–3692.

1602

F. Yekkalam Tash and B. Navayi Neya / Applied Mathematical Modelling 77 (2020) 1582–1602

[23] M. Eskandari-Ghadi, A complete solution of the wave equations for transversely isotropic media, J. Elast. 81 (2005) 1–19. [24] M. Eskandari-Ghadi, A. Ardeshir-Behrestaghi, Forced vertical vibration of rigid circular disc buried in an arbitrary depth of a transversely isotropic half space, Soil Dyn. Earthq. Eng. 30 (2010) 547–560. [25] M. Eskandari-Ghadi, A. Ardeshir-Behrestaghi, B. Navayi-Neya, Mathematical analysis for an axissymmetric disc-shaped crack in transversely isotropic half-space, Int. J. Mech. Sci. 68 (2013) 171–179. [26] M. Eskandari-Ghadi, M. Fallahi, A. Ardeshir-Behrestaghi, Forced vertical vibration of rigid circular disc on a transversely isotropic half-space, J. Eng. Mech. 136 (2009) 913–922. [27] M. Eskandari-Ghadi, A. Mirzapour, A. Ardeshir-Behrestaghi, Rocking vibration of a rigid circular disc in a transversely isotropic full-space, Int. J. Numer. Anal. Met. 35 (2011). [28] M. Eskandari-Ghadi, R.Y.S. Pak, A. Ardeshir-Behrestaghi, Elastostatic Green’s functions for an arbitrary internal load in a transversely isotropic bi-material full-space, Int. J. Eng. Sci. 47 (2009) 631–641. [29] M. Eskandari-Ghadi, S. Sattar, Axisymmetric transient waves in transversely isotropic half-space, Soil Dyn. Earthq. Eng. 29 (2009) 347–355. [30] M. Eskandari-Ghadi, S. Sture, R.Y.S. Pak, A. Ardeshir-Behrestaghi, A tri-material elastodynamic solution for a transversely isotropic full-space, Int. J. Sol. Struct. 46 (2009) 1121–1133. [31] M. Eskandari-Ghadi, R.Y.S. Pak, A. Ardeshir-Behrestaghi, Transversely isotropic elastodynamic solution of a finite layer on an infinite subgrade under surface loads, Soil Dyn. Earthq. Eng. 28 (2008) 986–1003.