A hybrid Bézier based multi-step method and differential quadrature for 3D transient response of variable stiffness composite plates

Composite Structures 154 (2016) 344–359

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

A hybrid Bézier based multi-step method and differential quadrature for 3D transient response of variable stiffness composite plates Y. Heydarpour, M.M. Aghdam ⇑ Department of Mechanical Engineering, Amirkabir University of Technology, Hafez Ave., Tehran 15875-4413, Iran

a r t i c l e

i n f o

Article history: Received 7 December 2015 Revised 17 June 2016 Accepted 22 July 2016 Available online 25 July 2016 Keywords: Three-dimensional Transient analysis Laminated plates Curvilinear fibers Layerwise-differential quadrature method Multi-step method

a b s t r a c t Three-dimensional (3D) transient analysis of a variable stiffness composite laminated (VSCL) plate with curvilinear fibers is presented. The fiber orientation angle changes linearly with respect to the horizontal coordinate in each layer. Various combinations of arbitrary boundary conditions at edges of the plate are assumed. Layerwise-differential quadrature method (LW-DQM) as an efficient and accurate numerical tool is employed to discretize the obtained governing equations for each layer in the spatial domain. Using the DQM enables one to accurately and efficiently discretize the partial differential equations and also implement the boundary and compatibility conditions in their strong forms. Then, solution of the resulting system of differential equations in the time domain, a novel multi-step method based on the Bézier curves is applied. Much lower computational time of the new multi-step with respect to the Newmark time integration scheme together with numerical stability of the method is presented. Apart from convergence and accuracy of the method, effects of variation in the fiber path, different geometrical parameters and boundary conditions on the transient response of the VSCL plates under dynamic loading are investigated. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Fiber reinforced composite materials have been widely used in different industries including aeronautic and astronautic technology, automobile and other modern technologies, mainly due to the high strength and light weight together with other improved properties. By selecting specific materials for matrix, fibers and also by adjusting the number of plies together with fiber orientation, one may find various opportunities to tailor parts. More recently, in order to improve performance of the composite materials, researchers suggested using variation of the stiffness in the space [1–24]. This category of laminates is identified as Variable Stiffness Composite Laminates (VSCL). The variable stiffness composites may be obtained by varying the fiber volume fraction in the laminate [1], dropping or adding plies to the laminate [4] or using curvilinear fibers [2,3,5–24]. The type of VSCL of interest herein is the one where curvilinear fibers are used. Improving the structural response and achieving additional weight reduction by variation of the fiber orientation [21–24] may consider as advantage of VSCL materials to replace the conventional straightfibers. ⇑ Corresponding author. E-mail address: [email protected] (M.M. Aghdam). http://dx.doi.org/10.1016/j.compstruct.2016.07.060 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

It is obvious that the production of composite laminates with rectilinear fibers is much easier than those with curvilinear fibers. However, previous studies exhibited that the curvilinear fibers as the reinforcements in the matrix may achieve better performance in some VSCL plates [2,3,5–24]. For example, numerical results demonstrated that using VSCL panels may lead to considerable increase in the amplitude of the optimal fundamental frequency [8]. Also, results indicated that using variable stiffness plates in comparison with constant stiffness can meaningfully change natural frequencies and mode shapes of vibration [11]. On the other hand, machines that efficiently arrange fibers with a variable orientation are somehow a recent challenge which involves many unknown issues related to mechanical behavior of this category of variable stiffness materials [14]. Hence, static, buckling, free vibration and transient analysis of structural elements made of VSCL with curvilinear fibers such as beams, plates and shells are essential for engineering design and manufacture. In recent years, the buckling and free vibration analyses of VSCL plates have been investigated by many researchers; see for instance [2,9–17,19,20]. In comparison to the research works in these topics, only few studies can be found on the static analysis of VSCL plates with curvilinear fibers [21–24]. In addition, the transient analysis of constant stiffness composite laminated (CSCL) plates and structures with straight fibers subjected to various

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

345

combinations of loadings and boundary conditions have been studied extensively in previous works [25–36]. Akhavan et al. [21] studied large deflection and stresses of variable stiffness composite laminated plates with a new p-version finite element, using third-order shear deformation theory. They considered several boundary conditions and different loading conditions. Yazdani et al. [22] presented a new p-version finite element, with hierarchic basis functions based on zig-zag layerwise theory to study the non-linear behavior of variable stiffness composite plates. Newton–Raphson method was employed to solve the equations of equilibrium. The plate was assumed to be fully clamped in their edges. In addition, they modeled the problem using the commercial finite element software ABAQUS. Yazdani and Ribeiro [23] developed a p-version finite element model for VSCL plates based on a layerwise plate theory to determine deflection and stress distribution of asymmetric composite laminates. Clamped boundary conditions were considered at all the edges. Also, results were compared with the finite element ABAQUS and ANSYS codes. Akbarzadeh et al. [24] investigated effects of transverse shear deformation and embedded manufacturing defects on the structural response of a variable stiffness plate made by Automated Fiber Placement (AFP). Static bending, buckling, and free vibration were considered. They obtained the governing equations based on the classical and shear deformation theories, which were solved for simply supported plates using hybrid Fourier–Galerkin method. To the authors’ best of knowledge, the three-dimensional transient analysis of VSCL plates subjected to dynamic loading has not been investigated in the open literature. The layerwise theory is a refined theory that can be employed to improve the displacement field and also strains and stresses with minimum computational cost. This theory assumes separate displacement field expansions within each subdivision. Hence, the layerwise theory provides a more accurate representation of the strain field in discrete layers compared to equivalent singlelayer theories. Since the shear strains are discontinuous, this leaves the possibility of the continuous transverse stresses between adjacent layers in the layerwise theory. Besides, layerwise theory provides capability of precise modeling of moderately thick and thick plates. In the present study, a linear fiber variation function, which, in spite of its simplicity, allows one to model plates with a reasonably wide range of properties, is considered. It is assumed that the individual layers are composed of curvilinear fibers surrounded by an isotropic matrix. Due to the intrinsic complexity of the three-dimensional formulation of composite laminated plates with curvilinear fibers, the layerwise-differential quadrature method (LW-DQM) as an efficient and more accurate numerical technique [35,37–44] is utilized to discretize the equations of motion in the spatial domain. This converts the governing partial differential

equations to a system of ordinary differential equations in the time domain. Then, to obtain transient response of the VCSL plate in the time space, the resulting ordinary differential equations (ODE) of motion are solved in the time domain using a newly developed multi-step technique based on the Bézier curves [45]. The present formulation and method of solution are validated and fast rate of convergence plus high accuracy with minimum computational cost is demonstrated. Finally, a detailed parametric study is performed to investigate effects of varying fiber orientation, different geometrical parameters and boundary conditions on the transient response of the VSCL plates which can be used as benchmark in future studies.

Fig. 1. Geometry of the laminated plate.

Fig. 2. Fiber orientation in a variable stiffness lamina.

2. Mathematical modeling Consider a rectangular laminated plate composed of NL perfectly bonded layers with curvilinear fiber which is subjected to dynamic load on the upper surface (Fig. 1). The variations of the field variables across the plate thickness are accurately modeled by dividing the plate into a set of co-axial mathematical layers in the thickness direction. In particular, it is assumed that each physical layer is divided into three mathematical layers. The material properties and thickness of each layer of the plate are assumed to be arbitrary. The length of the plate is denoted as a, width b and total thickness h, which is shown in Fig. 1. Also, the Cartesian coordinate system with coordinate variables (x, y, z) is used to label the material points of the laminated plate in the undeformed reference configuration. The displacement components of an arbitrary material point of the plate are indicated as u, v and w along the x, y and z-directions, respectively. 2.1. Model for VSCL plate with curvilinear fibers Unlike conventional composites with rectilinear fibers, curvilinear fiber paths cannot be described by a single orientation angle. A wide range of categories to describe the fiber path orientation can be considered [5]. It is assumed that, within each VSCL layer, the reference fiber path orientation varies linearly with x from one value T0 at the center to another value T1 at the edges, which is shown in Fig. 2. The orientation of a single curvilinear fiber path is indicated by [hT0|T1i], which represent the positive rotation angle of the principal material axes of the layer from the x-axis in the xy-plane. The orientation of the reference fiber path in each layer of the VSCL plate is given by [11],

hðxÞ ¼

2 ðT 1  T 0 Þ jxj þ T0 a

ð1Þ

There are two approaches to create the other fiber paths; i.e. parallel and shifted method [11]. The parallel method produces curvilinear fiber paths so that each path is defined as a set of points

346

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

Table 1 Comparison between the new multi-step formula and the Newmark’s scheme for  square VSCL [h30°, 60°i] plate ðt ¼ 0:01 s; a=h ¼ 5; Nm = 1, N x ¼ N y ¼ N z ¼ 7 . Boundary conditions

CCCC

SCSC

Nt

New multi-step formula

Newmark’s method

W

CPU time (s)

W

CPU time (s)

100 200 500 1000 1500 2000 2500

0.300 0.500 0.636 0.682 0.692 0.692 0.691

0.002155 0.004514 0.012062 0.031097 0.036592 0.061463 0.089006

0.333 0.515 0.644 0.686 0.696 0.694 0.692

21.9466 43.88788 109.8422 239.8415 329.8407 450.8151 680.2514

100 200 500 700 1000 1500 2000

0.852 0.988 1.119 1.151 1.150 1.149 1.149

0.00199 0.004359 0.011575 0.015553 0.02263 0.031692 0.058634

0.810 0.970 1.105 1.140 1.146 1.148 1.148

23.834767 44.375243 117.767857 156.065791 223.932962 317.563542 437.097345

limited and should not be larger than a deterministic value [11]. The curvature of the fiber orientation more than a specified value makes significant gaps and overlaps in VSCL plates [21]. 2.2. Effective material properties and the 3D constitutive relations Based on the 3D linear theory of elasticity, the constitutive relations, at an arbitrary material point of the kth physical layer can be summarized as [34],

2

8 ðkÞ 9 > r > > > > xx > > < rðkÞ > =

3 ðkÞ 8 ðkÞ 9 Q 16 > > e > > > 7> xx ðkÞ ðkÞ ðkÞ > < ðkÞ > = Q 22 Q 23 Q 26 7 7 eyy ; 7 ðkÞ ðkÞ ðkÞ 7> ðkÞ > Q 23 Q 33 Q 36 5> e zz > > > > ðkÞ > : ; ðkÞ ðkÞ ðkÞ 2exy Q 26 Q 36 Q 66 #( ðkÞ ) ðkÞ Q 45 2eyz

ðkÞ

ðkÞ

Q 11

Q 12

6 ðkÞ 6Q 6 12 ¼ 6 ðkÞ ðkÞ > > 6Q > > r zz > > > > 4 13 : ðkÞ ; ðkÞ rxy Q 16 ( ðkÞ ) " ðkÞ Q 44 ryz ¼ ðkÞ ðkÞ Q 45 Q 55 rxzðkÞ yy

ðkÞ

Q 13

ð2a; bÞ

2exz

ðkÞ

where rij and eij (i, j = x, y, z) are the stress and strain components, ðkÞ

ðkÞ

ðkÞ

respectively; and Q ij ði; j ¼ 1; 2; . . . ; 6Þ are the transformed compolying a constant distance from the reference curve, which changes the tow-orientation in different fibers and leads to more difficult computations. However, the shifting of the reference fiber path along the y axis to generate the remaining paths is more straight-forward than the use of parallel paths, i.e. the twoorientation is equal in all fibers of the lamina, leading to a relatively lower computations efforts. At each point along the fiber reference path, the orientation of the fibers should be chosen from an acceptable range, imposing a manufacturing constraint. Therefore, the curvature of fibers is

2.00

N t=50 Nt=200

nents of the layer stiffness [34]. The strain–displacement relations, based on the linear threedimensional theory of elasticity, at an arbitrary material point of the kth physical layer are: ðkÞ ðkÞ ðkÞ eðkÞ eðkÞ eðkÞ xx ¼ @u =@x; yy ¼ @ v =@y; zz ¼ @w =@z;

  1 @ v ðkÞ @wðkÞ þ ; 2 @z @y  ðkÞ  1 @u @ v ðkÞ ¼ þ 2 @y @x

eðkÞ yz ¼ eðkÞ xy

2.05

N t=400 N t=500

1.32

0.58

0.59

Nt=50 Nt=200

  1 @uðkÞ @wðkÞ þ ; 2 @z @x ð3a-fÞ

N t=400 N t=500

W

W

1.29

exzðkÞ ¼

-0.13

-0.14

-0.84

-0.87

-1.55

-1.60

0

2.15

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o Nt=50 N t=200

0.01

0

2.17

N t=400 N t=500

1.41

0.65

0.65

t (s) (b) T1 = 45o Nt=50 N t=200

0.01

N t=400 N t=500

W

W

1.40

0.002 0.004 0.006 0.008

-0.10

-0.11

-0.85

-0.87

-1.60

-1.63 0

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0.01

0

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

Fig. 3. Convergence of the transient response from Bezier method for different time steps for clamped square VSCL plates (a/h = 5, Nm = 3, T0 = 45°, Nx = Nz = Nz = 7).

347

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

2.05

Nz=3 N z=5

2.05

Nz=7 Nz=9

1.30

0.59

0.55

Nz=7 Nz=9

W

W

1.32

Nz=3 Nz=5

-0.14

-0.20

-0.87

-0.95

-1.70

-1.60 0

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o

2.20

Nz=3 Nz=5

0

0.01

0.002 0.004 0.006 0.008

t (s) (b) T1 = 45o

2.25

Nz=7 Nz=9

1.42

0.60

0.59

Nz=7 Nz=9

W

W

1.40

Nz=3 Nz=5

0.01

-0.20

-0.24

-1.00

-1.07

-1.90

-1.80 0

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0

0.01

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

Fig. 4. Convergence of the transient response against the number of DQ grid points for clamped square VSCL plates (a/h = 5, Nm = 3, T0 = 45°, Nx = Ny = Nz, Nt = 400).

2.3. Governing equations and solution procedure In order to accurately model the variation of the field variables across the plate thickness, the plate is divided into N m ðP N L Þ mathematical layers in the thickness direction. The governing differential equations of the eth mathematical layer together with the related external boundary conditions and also the compatibility conditions at the interface of layers (e) and (e + 1) are presented. The boundary conditions at the top and bottom surfaces of the VSCL plates are as, ðeÞ At z ¼ h=2 : rxz ¼ 0; ðeÞ ¼ 0; At z ¼ h=2 : rxz

ðeÞ ryz ¼ 0; rðeÞ zz ¼ 0 for e ¼ 1

ðeÞ ryz ¼ 0; rðeÞ zz ¼ qðtÞ for e ¼ N m

rðeÞ yy ðx; y; z; tÞ ¼ 0;

v ðeÞ ðx; y; z; tÞ ¼ 0 or rxyðeÞ ðx; y; z; tÞ ¼ 0; ðeÞ ryz ðx; y; z; tÞ ¼ 0

ð7a-fÞ

ðeÞ

ðeÞ

ðeÞ

2 ðeÞ 2 ðeÞ dQ 66 @uðeÞ @ 2 uðeÞ dQ 16 @uðeÞ ðeÞ @ u ðeÞ @ u þ Q þ Q þ þ 26 45 @x2 dx @x @y2 dx @y @z2 ðeÞ ðeÞ 2 ðeÞ 2 ðeÞ 2 ðeÞ ðeÞ ðeÞ dQ dQ @ v @ v @ v @ v ðeÞ ðeÞ ðeÞ @ v 66 26 þ Q 66 þ Q þ Q þ þ 22 44 @x2 dx @x @y2 dx @y @z2 ðeÞ 2 ðeÞ  @ 2 uðeÞ dQ 36 @wðeÞ  ðeÞ @ v ðeÞ ðeÞ þ þ Q 12 þ Q 66 þ 2Q 26 dx @z @x@y @x@y   @ 2 wðeÞ   @ 2 wðeÞ @ 2 v ðeÞ ðeÞ ðeÞ ðeÞ ðeÞ ¼0 þ Q 36 þ Q 45 þ Q 23 þ Q 44  qðeÞ @x@z @y@z @t2 ð10Þ

ðeÞ

ð6a-fÞ

ðeÞ

2 ðeÞ 2 ðeÞ @ 2 uðeÞ dQ 11 @uðeÞ dQ 16 @uðeÞ ðeÞ @ u ðeÞ @ u þ Q 66 þ Q 55 þ þ 2 2 @x dx @x @y dx @y @z2 ðeÞ ðeÞ 2 ðeÞ 2 ðeÞ 2 ðeÞ ðeÞ ðeÞ dQ 16 @ v dQ 12 @ v ðeÞ @ v ðeÞ @ v ðeÞ @ v þ Q 16 þ þ þ Q 26 þ Q 45 2 2 @x dx @x @y dx @y @z2 ðeÞ 2 ðeÞ 2 ðeÞ   ðeÞ dQ 13 @w ðeÞ @ u ðeÞ ðeÞ @ v þ þ 2Q 16 þ Q 12 þ Q 66 dx @z @x@y @x@y 2 ðeÞ   2 ðeÞ @ 2 uðeÞ ðeÞ ðeÞ @ w ðeÞ ðeÞ @ w ¼ 0 ð9Þ þ Q 36 þ Q 45  qðeÞ þ ðQ 13 þ Q 55 Þ @x@z @y@z @t2

ðeÞ

Q 16

At y ¼ b=2; b=2 :

either wðeÞ ðx; y; z; tÞ ¼ 0 or

ð8a-fÞ

t¼0

Q 11

rðeÞ xx ðx; y; z; tÞ ¼ 0; ðeÞ ðeÞ ðx; y; z; tÞ ¼ 0; either v ðx; y; z; tÞ ¼ 0 or rxy

either

t¼0

ð5a-cÞ

Either uðeÞ ðx; y; z; tÞ ¼ 0 or

Either uðeÞ ðx; y; z; tÞ ¼ 0 or

t¼0

ð4a-cÞ

At x ¼ a=2; a=2 :

rxzðeÞ ðx; y; z; tÞ ¼ 0

uðeÞ ðx; y; z; 0Þ ¼ v ðeÞ ðx; y; z; 0Þ ¼ wðeÞ ðx; y; z; 0Þ ¼ 0;    @uðeÞ  @ v ðeÞ  @wðeÞ  ¼ ¼ ¼ 0: @t  @t  @t 

where again e = 1, 2, . . ., Nm. Using Eqs. (2) and (3), the three-dimensional equations of motion of the VSCL plate can be derived in terms of the displacement components as,

Hereafter, the superscript ‘e’ is used to indicate the material properties and field variables of the eth mathematical layer. Also, the external boundary conditions at the edges of the VSCL plates areas,

either wðeÞ ðx; y; z; tÞ ¼ 0 or

where e = 1, 2, . . ., Nm. In this study, it is also assumed that before any loading, the VSCL plate is at rest in the undeformed configuration. Hence,

348

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

2.25

Nz=3

2.30

Nz=7 Nz=9

N z=5

1.48

0.63

0.66

Nz=7 Nz=9

W

W

1.44

Nz=3 Nz=5

-0.18

-0.16

-0.99

-0.98

-1.80

-1.80 0

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o

2.35

Nz=3

0.01

0

2.35

Nz=7 Nz=9

Nz=5

1.48

0.65

0.61

t (s) (b) T1 = 45o Nz=3 Nz=5

0.01

Nz=7 Nz=9

W

W

1.50

0.002 0.004 0.006 0.008

-0.20

-0.26

-1.05

-1.13

-2.00

-1.90 0

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0

0.01

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

Fig. 5. Convergence of the transient response against the number of DQ grid points for SCSC square VSCL plates (a/h = 5, Nm = 3, T0 = 45°, Nx = Ny = Nz, Nt = 400).

Table 2 Convergence of the non-dimensional center deflections against the number of mathematical layer for square VSCL [h30°, 60°i] plate ðt ¼ 0:01 s; a=h ¼ 5Þ. Boundary conditions

Nm

3

5

7

9

CCCC

1 3 5 7

0.237 0.688 0.691 0.692

0.593 0.691 0.692 0.692

0.692 0.692 0.692 –

0.692 0.692 – –

1 3 5 7

0.675 1.148 1.148 1.150

1.083 1.149 1.150 1.150

1.150 1.150 1.150 –

1.150 1.150 – –

SCSC

Nx = Ny = Nz

  ðeÞ  ðeÞ  2 ðeÞ dQ 55 @uðeÞ @wðeÞ dQ 45 @ v ðeÞ @wðeÞ ðeÞ @ w þ þ þ þ Q 55 dx @z @x dx @z @y @x2 2 ðeÞ 2 ðeÞ 2 ðeÞ   ðeÞ @ w ðeÞ @ w ðeÞ ðeÞ @ u þ Q 44 þ Q 33 þ Q 13 þ Q 55 @y2 @z2 @x@z   2 ðeÞ   2 ðeÞ ðeÞ ðeÞ @ u ðeÞ ðeÞ @ v þ Q 36 þ Q 45 þ Q 36 þ Q 45 @y@z @x@z 2 ðeÞ   2 ðeÞ @ 2 wðeÞ ðeÞ ðeÞ @ v ðeÞ @ w þ Q 23 þ Q 44 ¼0 þ 2Q 45  qðeÞ @y@z @x@y @t 2

ð11Þ

where Eqs. (9)–(11) are the equations of motion along the x, y and zdirections, respectively. The geometrical and natural compatibility conditions at the interface of two adjacent layers of the VSCL plate are:

Table 3 Convergence and accuracy of maximum non-dimensional deflections of simply-supported laminated composite plates under uniformly distributed loading (a/b = 1). a/h

Nx = Nz

Analytical [34]

5

7

9

11

13

Orthotropic plate

10 20 100

0.9681 0.7331 0.6568

0.9548 0.7299 0.6566

0.9548 0.7299 0.6566

0.9548 0.7299 0.6566

0.9548 0.7299 0.6566

0.9519 0.7262 0.6528

(0/90)

10 20 100

1.9754 1.7935 1.7354

1.9608 1.7735 1.7167

1.9535 1.7696 1.7102

1.9498 1.7625 1.7035

1.9498 1.7625 1.7035

1.9468 1.7582 1.6980

(0/90/0)

10 20 100

1.1640 0.8013 0.6767

1.0784 0.7813 0.6745

1.0265 0.7604 0.6728

1.0258 0.7596 0.6723

1.0258 0.7596 0.6723

1.0219 0.7572 0.6697

349

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359 Table 4 Accuracy of the results for (0/90)s square laminated composite plate under a uniform load (a/h = 50). Boundary Conditions P ^ xx ð0; 0; h=2Þ

Levy’s method

SemiAnalytical [46]

SSSS

0.8228

0.8228

0.8228

CSCS

0.3166

0.3166

0.3165

SSSS

0.3560

0.3560

0.3560

CSCS

0.0457

0.0457

0.0456

^ zz ða=4; a=4; 0Þ R

SSSS

0.5004

0.5005

0.5002

CSCS

0.5005

0.5005

0.5005

^ xy ða=4; a=4; h=4Þ R

SSSS

0.0079

0.0079

0.0078

CSCS

0.0021

0.0021

0.0021

P ^ yy ð0; 0; h=4Þ

^ xz ða=4; a=4; 0Þ R ^ yz ða=4; a=4; 0Þ R

SSSS

0.2347

0.2347

0.2346

CSCS

0.2968

0.2968

0.2966

SSSS

0.0300

0.0299

0.0303

CSCS

0.0350

0.0352

0.0345

    ðeÞ ðeþ1Þ uðeÞ z2 ¼ uðeþ1Þ z1 ;     ðeÞ ðeþ1Þ ; wðeÞ z2 ¼ wðeþ1Þ z1

  ðeÞ   ðeÞ ðeÞ ðeÞ @ v ðeÞ  ðeÞ @u ðeÞ @ v ðeÞ @w ðeÞ @u Q 13 þ Q 23 þ Q 33 þ Q 36 þ  ðeÞ @x @y @z @y @x z¼z2  ðeþ1Þ ðeþ1Þ ðeþ1Þ ðeþ1Þ @u ðeþ1Þ @ v ðeþ1Þ @w ¼ Q 13 þ Q 23 þ Q 33 @x @y @z  ðeþ1Þ   ðeþ1Þ  @v ðeþ1Þ @u  ð12a-fÞ þ þQ 36  ðeþ1Þ : @y @x z¼z 1

ðkÞ z1

ðeÞ

ðkÞ z2

where and (k = e, e + 1) denote the thickness coordinate of the lower and upper surfaces of the kth mathematical layer, respectively. Different types of the classical boundary conditions at the edges of the plate can be assumed which can be achieved by combining the conditions stated in Eqs. (6) and (7). For instance, for the VSCL plate with clamped and simply supported at the edges x = a/2 and a/2, one has Clamped edge: u(e)(x, z, t) = 0, v(e)(x, z, t) = 0, w(e)(x, z, t) = 0.

    ðeÞ ðeþ1Þ z2 ¼ v ðeþ1Þ z1 ;

 ðeÞ   @w @wðeÞ  ðeÞ @u ðeÞ @ v Q 55 þ þ Q 45 þ  ðeÞ @z @x @z @y z¼z2   ðeþ1Þ     ðeþ1Þ ðeþ1Þ @w @wðeþ1Þ  ðeþ1Þ @u ðeþ1Þ @ v ¼ Q 55 þ þ Q 45 þ  ðeþ1Þ ; @z @x @z @y z¼z



1

ðeÞ ðeÞ @uðeÞ ðeÞ @ v ðeÞ @w þ Q 12 þ Q 13 @x @y @z  ðeÞ  @ v ðeÞ ðeÞ @u ¼ 0; þ Q 16 þ @y @x ðeÞ

Simply supported edge :Q 11

v ðeÞ ðx; z; tÞ ¼ 0;

1.2

wðeÞ ðx; z; tÞ ¼ 0:

1.0

ð13a; bÞ

Present Analytical [29]

0.8

w (in)

0.6

0.4

-0.1

0.2

-0.2 Present Analytical [29]

-0.5 0.000

0.002

0.004

0.006

t (s) (a) step loading

0.008

-0.6 0.000

0.002

0.004

0.006

t (s) (b) triangular loading

1.0

0.7

0.008

Present Analytical [29]

0.6

0.5

w (in)

ðeÞ

w (in)



v ðeÞ

1

w (in)



Present

  ðeÞ   ðeÞ   @wðeÞ @wðeÞ  ðeÞ @u ðeÞ @ v Q 45 þ þ Q 44 þ  ðeÞ @z @x @z @y z¼z2   ðeþ1Þ     ðeþ1Þ ðeþ1Þ @w @wðeþ1Þ  ðeþ1Þ @u ðeþ1Þ @ v þ Q 44 ¼ Q 45 þ þ  ðeþ1Þ ; @z @x @z @y z¼z

0.3

0.2

0.1

-0.2 Present

-0.2 0.000

Analytical [29]

0.002

0.004

0.006

t (s) (c) sine loading

0.008

-0.6 0.000

0.002

0.004

0.006

t (s) (d) explosive blast loading

0.008

Fig. 6. Comparison of the center deflection for the simply supported [0/90/0] laminated plates subjected to various loadings.

350

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

37.0

28.8

25.8

17.6

14.6

Σ xx

Σ xx

40.0

6.4

3.4

-4.8

-7.8 Present

Present

Analytical [29]

-16.0 0.000

0.002

0.004

t (s)

Analytical [29]

0.006

0.008

-19.0 0.000

(a) step loading

0.004

t (s)

0.006

0.008

(b) triangular loading 36.0

16.5

24.8

11.2

13.6

Σ xx

Σ xx

21.8

5.9

2.4

0.6

-8.8

Present Analytical [29]

-4.7 0.000

0.002

0.002

0.004

t (s)

0.006

0.008

-20.0 0.000

(c) sine loading

Present Analytical [29]

0.002

0.004

t (s)

0.006

0.008

(d) explosive blast loading

P Fig. 7. Comparison of the non-dimensional in-plane normal stress ð xx ¼ rxx =q0 Þ at center of the top surface of the plate for the simply supported [0/90/0] laminated plates subjected to various loadings.

Table 5 Different loading function Q(t) [29]. Loading Q(t)

Step

1 t 6 t1 0 t > t1

Triangular

1  t=t1 t 6 t 1 0 t > t1

It is evident that presenting exact solution for the resulted system of differential equations with variable coefficients is very difficult if it is not impossible to obtain such a solution. Therefore, an approximate method should be used to obtain numerical results. On the other hand, the differential quadrature method (DQM) as a simple, efficient and accurate numerical method has been confirmed for solving different structural problems in recent years; see for example [35,37–44]. Thus, in the present study, a layerwise-differential quadrature method (LW-DQM) is employed to discretize the governing differential equations of motion and the related boundary and compatibility conditions of the mathematical layers in the spatial domain. Here, three mathematical layers are assumed in each physical layer. Based on the differential quadrature method, each mathematical layer is discretized into a set of Nx, Ny and Nz grid points along the x, y and z directions, respectively. Using the DQ rules, the strong form of the equations of motion within the solution domain together with boundary conditions along boundary grid points are discretized. Furthermore, the compatibility conditions of the layers at the interface of two adjacent layers are properly discretized [40].

Sine

sinðpt=t 1 Þ 0

Explosive blast t 6 t1 t > t1

ec1 t

This procedure transforms the governing partial differential equations and related boundary and compatibility conditions into a system of second-order ordinary differential equations (ODE) in temporal domain, which can be demonstrated in the matrix form as:

½M

( 2 ) d D dt

2

þ ½KfDg ¼ fFðtÞg

ð14Þ

where [M] and [K] are the mass and stiffness matrices while [F] is the load vector. Also, {D} refers to displacement components vector of grid points. Based on the element arrangement of this vector, the elements of the coefficient matrices [M], [K] and the load vector [F] are determined from the discretized form of the governing differential equations and the related boundary and compatibility conditions. Various numerical methods are presented for the solution of the system of ordinary differential equations (ODEs) including single step approaches such as Euler and modified Euler methods, Taylor series and Runge–Kutta method together with multi-step methods such as Adams-Moulton and Adams-Bashforth.

351

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

2.18

o

T1=60 o T1=90

o

T1=45

1.41

2.20

o

T1=30

o

o

T1=30

T1=60 o T1=90

o

T1=45

1.44

0.64

W

W

0.68

-0.13

-0.08

-0.90

-0.84

-1.67

-1.60

0

0.002 0.004 0.006 0.008

t (s)

0.01

0

o (a) T0 = 30

2.05

o

T1=45

1.32

2.00

0.01

o

o

T1=30

T1=60 o T1=90

o

t (s)

o (b) T0 = 45

o

T1=30

0.002 0.004 0.006 0.008

T1=60 o T1=90

o

T1=45

1.30

0.59

W

W

0.60

-0.14

-0.10

-0.87

-0.80

-1.60

-1.50

0

0.002 0.004 0.006 0.008

t (s)

0.01

(c) T0 = 60

o

0

0.002 0.004 0.006 0.008

t (s)

0.01

(d) T0 = 90

o

Fig. 8. Influence of varying the fiber orientation on the transient response of clamped square VSCL plates (a/h = 5, Nm = 3, Nx = Ny = Nz = 7, Nt = 400).

In this study, a novel multi-step method based on the Bézier curves is employed to solve the system of differential Eq. (14) subjected to related initial conditions, as a system of initial value problems, which is adopted from the work of Aghdam et al. [45]. They used the general approach of the multi-step techniques together with various orders of Bézier curves to develop a new multi-step method. They proved the consistency, stability, and convergence of the new technique by using the various stability and error analysis theorems of multi-step methods. They also provided various initial value physical problems in which significant improvement, mainly in terms of computational cost, in comparison with the well-known Adams–Moulton method is demonstrated. In this study, the idea of using Bézier curves is extend to a rather complex system of ODE’s in a structural problem and the efficiency of the method is examined through comparisons with the well-known Newmark algorithm. In order to apply this Bézier based multi-step method to solve the Eq. (14), the system of second-order ordinary differential equation (ODE) should be transformed to a system of first order system of ODEs by definition new variables {g1} = {D} and fg2 g ¼ dtd fDg as:

(

g

g

d f 1g ¼ f 2g dt ½M dtd f 2 g þ ½Kf 1 g

g

g

¼ fFðtÞg

ð15a; bÞ

In the presented work, a four-step formula [45] is employed to solve the Eq. (15), see Appendix A. Therefore, the information of the first four points is required to start solution procedure. Since the first point is the known initial condition, the different step formulas can be to determine the next three steps, which are

presented in the Appendix A. Finally, by solving the system of ordinary differential equations (15), one may obtain the displacement components at any time step for all spatial grid points. 3. Numerical results As the first step, the formulation and solution method is evaluated by investigating the fast rate of convergence for transient analysis of the VSCL plates. In addition, the accuracy of the approach is demonstrated by performing comparison studies with other available data in the limit cases. Then, parametric studies for transient response of the VSCL plates with variable fibers angle for a single ð½hT 0 jT 1 iÞ and symmetric three ð½hT 0 jT 1 i; hT 0 þ 90 jT 1 þ 90 i; hT 0 jT 1 i) layers are carried out. Various dynamic loads can be considered at the upper surface of the plate. Without loss of generality, it is assumed that the VSCL plate is subjected to a conventional blast, i.e. combination of triangular and blast loads, while the load is considered to be uniformly distributed on the plate, which can be expressed as [32],

  t a1 t=t2 e qðtÞ ¼ q0 1  t2

ð16Þ

where t 2 ¼ 0:004 s and a1 ¼ 1:98 [32]. It should be noted that in this study, damping of the system is ignored in the governing equations and therefore, the energy is not dissipated within the cycles. However, decaying behavior in the response, after an initial excitation, is just an illusion due to the transient force. Also, except in the comparison studies, the non-dimensional parameters are defined as,

352

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

2.70

o

T1=60 o T1=90

o

T1=45

1.74

2.40

o

T1=30

o

o

T1=30

T1=60 o T1=90

o

T1=45

1.52

0.78

W

W

0.64

-0.18

-0.24

-1.14

-1.12

-2.10

-2.00

0

0.002 0.004 0.006 0.008

t (s)

0.01

0

(a) T0 = 30

o

2.35

t (s)

0.01

(b) T0 = 45

o

o

2.20

o

T1=30

T1=45

o

o

T1=30

T1=60 o T1=90

o

1.50

0.002 0.004 0.006 0.008

T1=60 o T1=90

o

T1=45

1.40

0.60

W

W

0.65

-0.20

-0.20

-1.05

-1.00

-1.80

-1.90 0

0.002 0.004 0.006 0.008

t (s)

0.01

(c) T0 = 60

0

0.002 0.004 0.006 0.008

t (s)

0.01

(d) T0 = 90

o

o

Fig. 9. Influence of varying the fiber orientation on the transient response of SCSC square VSCL plates (a/h = 5, Nm = 3, Nx = Ny = Nz = 7, Nt = 400).

3

W ¼ 100 X ii

X yz

E2 h

4

wð0; 0; 0Þ;

n¼

2z þ h ; 2h

q b 0  h h ¼ rii 0; 0; ði ¼ x; y; zÞ 2 q0 b   b h ¼ ryz 0;  ; 0 2 q0 b

ð17a-dÞ

In this research work, the individual layers are taken to be of equal thickness. Also, in the examples, material properties for the transient analysis of VSCL plate are chosen from the work of Reddy [34], which are,

E2 ¼ 1 GPa; E1 ¼ 25E2 ;

G12 ¼ G13 ¼ 0:5E2 ;

G23 ¼ 0:2E2 ;

q ¼ 1389:297 kg=m ; m12 ¼ m13 ¼ m23 ¼ 0:25 3

For the sake of brevity, a counter-clock wise notation starting from the edge x = a/2 (see Fig. 2) is used to identify the boundary conditions of the plate. Therefore, the symbolism ‘SCFC’, for instance, refers to a plate with simply supported edge at x = a/2, clamped at y = b/2, free at x = a/2 and clamped at y = b/2. As the first example, a comparison between the rate of convergence of the presented formula based on the Bézier curves and the Newmark algorithm in conjunction with Galerkin scheme for a single layer variable stiffness composite plate with two different boundary conditions is tabulated in Table 1. It is obvious that the CPU time requirement of the presented multi-step technique is considerably lower than the Newmark method by a factor of four digits. For instance, as an example, the computational times for analyses of clamped VSCL plate with Nt = 100, which Nt is the

number of time steps, are 0.002155 and 21.9466s for the Bézier based multi-step and Newmark methods, respectively. Apart from computational cost of the presented method which is significantly lower than the Newmark algorithm, the method shows high numerical stability. Furthermore, from accuracy point of view, presented results reveal that the error of the presented method is much less than that of the Newmark algorithm. In Fig. 3, the convergence of the non-dimensional central deflection parameter of the VSCL plate for various fiber orientations and different time steps are shown. In this Figure, very fast rate of convergence and numerical stability of the presented method based on the Bézier curves are quite evident. The convergence of the present approach for transient analysis of a VSCL plate with two different boundary conditions against the number of DQ grid points along the x, y and z-direction are presented in Figs. 4 and 5. The influence of different fiber orientation on the convergence behavior of the approach is also investigated in these Figures. Also, convergence of the method against the number of mathematical layers is exhibited in Table 2. In this Table, transient response of a single layer variable stiffness composite plate with two different boundary conditions is presented. It is noticeable that in all cases, the non-dimensional central deflection of the variable stiffness composite plate with curvilinear fibers converges rapidly without any numerical instability by decreasing the number of time steps and increasing the number of DQ grid points or otherwise the number of mathematical layers. In order to show the accuracy of the formulation and solution method, the maximum non-dimensional deflections for the static analysis of simply-supported laminated composite square plates under uniformly distributed loading, are compared with the

353

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

3.49

o

T1=30

o

T1=45

2.29

o

T1=60 o T1=90

o

0.35

o

T1=30

T1=60 o T1=90

o

T1=45

0.21

Σ xx

Σ xx

1.09

0.07

-0.11

-0.07

-1.31

-0.21

-2.51

-0.35 0

0.004 0.008 0.012 0.016

t (s)

0.02

0

0.004 0.008 0.012 0.016

t (s)

o (a) T0 = 30

1.50

o

T1=30

o

T1=45

0.96

0.02

o (b) T0 = 90 o

3.95

o

T1=60 o T1=90

o

T1=45

2.56

0.42

o

T1=30

T1=60 o T1=90

Σ yy

Σ yy

1.17

-0.12

-0.22

-0.66

-1.61

-1.20

-3.00 0

0.004 0.008 0.012 0.016

t (s)

0.02

0

0.004 0.008 0.012 0.016

o (c) T0 = 30

t (s)

0.02

o (d) T0 = 90

Fig. 10. Influence of varying the fiber orientation on the non-dimensional normal stress of SCSC square VSCL plates (a/h = 5, Nm = 3, N x ¼ N y ¼ 9; N z ¼ 7; Nt = 400).

analytical solution of Reddy [34] in Table 3. He formulated the problem based on the first-order shear deformation theory (FSDT) of plates. The results are obtained for three different values of the length-to-thickness ratio. In all cases, fast rate of convergence of the method for determination of the non-dimensional deflections and excellent agreement with analytical solution [34] is obvious. Table 4 compares the results of the present method with the semi-analytical solutions given by Naserian Nik and Tahani [46] for the [0°/90°]s square laminated composite plate subjected to a uniform load on its top surface. Also, they used Levy solution to obtain the stress components. In this example, material properties and non-dimensional stress components are [46],

cross-ply laminates using a higher-order plate theory (HSDT). The geometrical and material properties of the plates are [29],

E2 ¼ 12 GPa; E1 ¼ 25E2 ;

The value of the transverse load is taken to be q0 ¼ 68:9476 MPa. Also, it is assumed that the load varies with time according to one of the expressions given in Table 5, where t1 ¼ 0:006 s and for explosive blast loading c1 ¼ 330 s1 [29]. One can see that the central deflection and non-dimensional in-plane P normal stress ð xx Þ at center of the top surface of the plate are in very good agreement with the exact solutions [29] for all values of the time interval. In these examples, efficiency of the presented hybrid solution method consist of DQ and Bézier based multi-step formulation in dynamic analysis of laminated plates subjected to various transient loads including step, triangular, sine, and explosive blast loads is confirmed. After validation of the presented hybrid solution method, some parametric studies are carried out to exhibit the influences of varying fiber orientation, thickness, aspect ratio, different boundary conditions on the transient behavior of the composite laminated plates with curvilinear fibers subjected to dynamic load.

m12 ¼ 0:25 



G12 ¼ G13 ¼ 0:5E2 ;

^ xx ; R ^ yy ; R ^ xy ¼ ðrxx ; ryy ; rxy Þ R





^ xz ; R ^ yz ¼ ðrxz ; ryz Þ R

h ; q0 b

h

G23 ¼ 0:2E2

2

q0 b

2

^ zz ¼ R

;

rzz q0

ð18a-fÞ

Again, very good agreement can be seen between the results of the present method and those of the Levy method and semianalytical solution [46] for two different boundary conditions. Applicability and efficiency of the presented approach for dynamic analysis of laminated composite plates is verified by the transient solutions of a symmetric [0°/90°/0°] laminated composite plate with simply supported edges for various loadings in Figs. 6 and 7, which are compared with the work of Khdeir and Reddy [29]. They obtained exact solutions for transient responses of

a ¼ b ¼ 5 h; h ¼ 0:1524 m; E2 ¼ 6:895 GPa; E1 ¼ 172:369 GPa; G12 ¼ G13 ¼ 3:448 GPa; G23 ¼ 1:379 GPa;

q ¼ 1603:03 kg=m3 ; m12 ¼ m13 ¼ m23 ¼ 0:25 The plate is subjected to a distributed sinusoidal load in spatial domain as,

qðx; y; tÞ ¼ q0 sin

p 2

þ

p x a

sin

p 2

þ

p y b

Q ðtÞ

ð19Þ

354

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

0.50

o

o

T1=30

T1=60 o T1=90

o

T1=45

0.34

o

o

T1=30

T1=60 o T1=90

o

T1=45

0.29

0.16

Σ yz

Σ yz

0.18

0.42

0.02

0.03

-0.14

-0.10

-0.30

-0.23 0

0.004 0.008 0.012 0.016

t (s)

0.02

0

0.004 0.008 0.012 0.016

t (s)

(a) T0 = 30 0.35

o

T1=30

o

T1=45

0.24

0.02

(b) T0 = 45

o

o

o

T1=60 o T1=90

o

0.25

T1=30

o

T1=45

o

T1=60 o T1=90

0.17

Σ yz

Σ yz

0.13

0.09

0.02

0.01

-0.09

-0.07

-0.20

-0.15 0

0.004 0.008 0.012 0.016

t (s)

0.02

0

0.004 0.008 0.012 0.016

t (s)

o (c) T0 = 60

o (d) T0 = 90

Fig. 11. Influence of varying the fiber orientation on the non-dimensional transverse shear stress ( Nt = 400).

The effects of the fiber orientation angles, as an important parameter, on the transient response of the three layered variable stiffness composite plates are studied. Results for different fiber orientation angles of fully clamped and SCSC boundary conditions are respectively, demonstrated in Figs. 8 and 9. It can be seen that the fiber orientation angles has not significant effect on the nondimensional central deflection of the fully clamped VSCL plate, except for T0 = 30°. Unlike clamped VSCL plate, it is obvious from Fig. 9 that this parameter significantly changes the transient response of VSCL plates with SCSC boundary conditions. The SCSC square VSCL plates with T0 = 90° can experience less deflection than CSCL in the range of 12–25%. In addition, it is observable from this figure that increasing the fiber angle at plate ends (T1), the non-dimensional central deflection increases for all values of the time level, which implies reduction of total plate stiffness. It is interesting to note that the effects of the fiber orientation angles for both boundary conditions at T0 = 30° are considerable. The Fig. 9 show that for the different values of T0, the differences of maximum deflection are varies between 7.8% and 18%. Also, for different values of the fiber angle at plate ends (T1), increasing T0 leads increase of the overall stiffness of the plate. For instance, at T1 = 60°, the maximum deflection for VSCL plate with T0 = 90° is reduced about 13.5% in comparison with T0 = 30°. In both Figures, no significant qualitative difference is observed in the central deflection of the plate by increasing the value of T0, while in quantitative terms the differences are evident. In another study, transient non-dimensional normal and transverse shear stress components with various fiber orientation angles are plotted in Figs. 10 and 11, respectively. It is noticeable

0.02

P

yz)

of SCSC square VSCL plates (a/h = 5, Nm = 3, N x ¼ N y ¼ 9; N z ¼ 7;

that with changes in the curvilinear fiber paths, the stress components vary substantially. Also, it is obvious that the effects of the fiber orientation angles at T1 = 90° are more noticeable than other values of the fiber angle at plate ends. Moreover, in all cases, maximum normal and transverse shear stress components increase by increasing the fiber angle at plate ends (T1). As mentioned earlier, the system does not show any dissipation of energy as any type of damping is ignored in the governing equations. In this section, the response to impact loads of a VSCL plate is analyzed in more detail. Fig. 2 shows that; As shown in Figs. 8–11, varying T0 and T1 have significant influence on the values of the deflection and stress components. The reason is that composite plates with rectilinear fibers in each layer have uniform stiffness and strength whereas VSCL laminates with curvilinear fibers have non-uniform stiffness and strength [18]. The variable stiffness and strength affect the stress re-distribution and cause damage initiation. Therefore using VSCL instead of CSCL plates may decrease deflection and reduce possibility of damage onset and change in the position of damage initiation. Consequently, improved performance may achieve in some cases for VSCL plate. The influence of thickness of the VSCL plate on the time history of non-dimensional central deflection of the VSCL plate with SCSC boundary condition is shown in Fig. 12. It is obvious that by increasing this geometrical parameter, the non-dimensional central deflection reduces for all values of the time level. In addition, it is evident that with increasing the values of the fiber angle at plate ends (T1), the non-dimensional displacement increases for all values of the thickness of the VSCL plate.

355

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

7.15

7.00 Plate thickness 0.5h h 2h

4.40

Plate thickness 0.5h h 2h

4.52

1.89

W

W

1.80

-0.80

-0.74

-3.40

-3.37

-6.00

-6.00 0

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o

7.15

0

0.01

7.25

Plate thickness 0.5h h 2h

4.52

t (s) (b) T1 = 45o

0.01

Plate thickness 0.5h h 2h

4.60

1.95

W

W

1.89

0.002 0.004 0.006 0.008

-0.74

-0.70

-3.37

-3.35

-6.00

-6.00 0

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0.01

0

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

i 3 4 Fig. 12. Influence of the thickness on the transient response of SCSC square VSCL plates ½a=h ¼ 5; T0 = 45°, Nm = 3, Nx = Ny = Nz = 7, Nt = 400, W ¼ 100wðE2 h =q0 b Þ .

2.90

a/b=0.5 a/b=1

2.90

a/b=1.5

1.82

0.78

0.74

a/b=1.5

W

W

1.84

a/b=0.5 a/b=1

-0.28

-0.34

-1.34

-1.42

-2.50

-2.40 0

2.90

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o a/b=0.5 a/b=1

0

0.01

2.90

a/b=1.5

1.80

0.74

0.70

t (s) (b) T1 = 45o a/b=0.5 a/b=1

0.01

a/b=1.5

W

W

1.82

0.002 0.004 0.006 0.008

-0.34

-0.40

-1.42

-1.50

-2.50

-2.60 0

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0.01

0

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

Fig. 13. Influence of the aspect ratio on the transient response of SCSC square VSCL plates (a/h = 5, T0 = 45°, Nm = 3, Nx = Ny = Nz = 7, Nt = 400).

356

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

1.00

0.85 0.62

0.24

0.15

W

W

0.50

-0.14

-0.20

CCCC

CCCC

-0.55 -0.90 0

SCCC SSCC

-0.52

SSSC SSSS

-0.90

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o

SCCC SSCC SSSC SSSS

0

0.01

1.15

0.69

0.68

0.28

0.21

t (s) (b) T1 = 45o

0.01

W

W

1.10

0.002 0.004 0.006 0.008

-0.26

-0.13

CCCC

CCCC SCCC SSCC

-0.54

SSSC SSSS

-0.95 0

SCCC SSCC

-0.73

SSSC SSSS

-1.20

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0

0.01

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

Fig. 14. Influence of different boundary conditions on the transient response of square VSCL plates (a/h = 10, T0 = 45°, Nm = 3, Nx = Ny = Nz = 7, Nt = 400).

2.05

2.15

One layer Three layers

One layer Three layers

1.38

0.55

0.61

W

W

1.30

-0.20

-0.16

-0.95

-0.93

-1.70

-1.70 0

0.002 0.004 0.006 0.008

t (s) (a) T1 = 30o

0.01

0

2.30

0.002 0.004 0.006 0.008

t (s) (b) T1 = 45o

0.01

2.35 One layer

One layer

Three layers

Three layers

0.70

0.73

W

1.54

W

1.50

-0.10

-0.08

-0.90

-0.89

-1.70

-1.70 0

0.002 0.004 0.006 0.008

t (s) (c) T1 = 60 o

0.01

0

0.002 0.004 0.006 0.008

t (s) (d) T1 = 90o

0.01

Fig. 15. Comparison of the transient response for clamped square VSCL plates with single and three layers (a/h = 5, T0 = 45°, Nx = Ny = Nz = 7, Nt = 400).

357

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

2.60

2.60

1.75

One layer

One layer

Three layers

Three layers

1.76

0.90

Σ yy

Σ yy

0.92

0.05

0.08

-0.80

-0.76

-1.65

-1.60 0

0.004 0.008 0.012 0.016

0.02

t (s)

(a) T1 = 30 0.33

0

0.004 0.008 0.012 0.016

(b) T1 = 90 0.38

One layer

0.02

t (s)

o

o

One layer

Three layers

Three layers

0.25

0.09

0.12

Σ yz

Σ yz

0.21

-0.03

-0.01

-0.15

-0.14

-0.27

-0.27 0

0.004 0.008 0.012 0.016

0.02

t (s)

(c) T1 = 30

0

0.004 0.008 0.012 0.016

0.02

t (s)

(d) T1 = 90

o

o

Fig. 16. Comparison of the non-dimensional stress components of SCSC square VSCL plates with single and three layers (a/h = 5, T0 = 45°, N x ¼ N y ¼ 9; N z ¼ 7; Nt = 400).

o

T1=60 o T1=90

o

T1=45

o

0.011

o

T1=30

0.026

o

T1=30

T1=60 o T1=90

o

T1=45

0.002

Σ xz

Σ yz

0.007

-0.007

-0.012

-0.016

-0.031

-0.025

-0.050 0

0.2

0.4

0.6

ξ

0.8

0

1

0.2

0.4

ξ

0.6

0.8

1

0.00

Σ zz

-0.05

-0.10

-0.15 o

o

T1=30

T1=60 o T1=90

o

T1=45

-0.20 0

0.2

0.4

ξ

0.6

0.8

1

Fig. 17. Influence of varying fiber orientation on the non-dimensional stress of SCSC square VSCL plates [a/h = 5, t ¼ 0:01 ðsÞ; T0 = 45°, Nm = 3, Nx = Ny = Nz = 9, Nt = 400].

358

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359

The effect of the aspect ratio as an important geometrical parameter on the time history of the results for the VSCL plate with SCSC boundary condition is given in Fig. 13. It can be concluded that for the VSCL plates with a given thickness, non-dimensional central deflection increases by increasing aspect ratio of VSCL plates for all values of T1. In addition, as with previous results, by increasing values of the fiber angle at plate ends, overall stiffness of the VSCL plates decreases. The influence of the different boundary conditions on the time history of the transient response for the VSCL plates is performed in Fig. 14. For the VSCL plates subjected to dynamic load, different boundary conditions have significant effects on the transient response. Furthermore, it can be seen that by increasing the boundary constraints, the non-dimensional central deflection reduces for all values of the time levels. Also, in all types of boundary conditions, increasing the fiber angle at the plate ends decreases the stiffness of the VSCL plates. Comparisons between the transient non-dimensional central deflection, normal and transverse shear stress components of variable stiffness composite plates with single and three layers are depicted in Figs. 15 and 16. From Fig. 15, it can be seen that with the same total thickness, the VSCL plates composed of one layer have lower non-dimensional central deflection than the symmetric VSCL plates with three layers at each time interval ‘‘t”. For instance, the maximum difference between deflection magnitudes of one layer and three layers VSCL plates at T1 = 30° is approximately 10.5%. Also, it is evident from Fig. 16 that behavior of transient normal and transverse shear stress components for single and three layers are different. It is interesting to note that by increasing the number of layers of VSCL plate from one to three, with the same total thickness, the normal stress decrease and transverse shear stress increase. In order to take more advantage of the presented threedimensional analysis, the last example is dedicated study through thickness variation of various stress components at the point (x = 0, y = 0) with various fiber orientation angles for the SCSC boundary condition. Predictions of the normal and transverse shear stress components of the VSCL plate are exhibited in Fig. 17. Results reveal that the changes in the fiber path lead to considerable effect on the transverse shear stress components of the VSCL plate. However, the non-dimensional transverse normal stress (Rzz) for different fiber orientations remains almost the same. This change in transverse shear stresses magnitudes may be exploited to improve damage resistance in particular applications.

Then, parametric studies are carried out to investigate transient response of VCSL plates subjected to dynamic load with detailed discussions on the results. Results clearly demonstrate that performance of the VSCL is not always better than that of CSCL. However, it is shown that for specified boundary conditions, for example SCSC plate, the changes in the fiber curvature have significant effects on the transverse displacement and stress components. In all cases by increasing the value of T0, lower displacements are obtained, especially when the fiber angle at plate ends (T1) is small. Other interesting achievement of this study is that the maximum normal and transverse shear stress components increase by increasing the value of T1. This is due to the fact that conventional composite with straight fibers (CSCL) have uniform stiffness and strength whereas VSCL laminates have non-uniform stiffness and strength. Different transient responses of the VSCL plate are obtained by varying geometrical parameters; i.e. thickness and aspect ratio, boundary conditions and fiber orientations. The maximum decrease of non-dimensional deflection magnitudes is 13.5% for the VSCL plate at the same values of T1 = 60° with T0 = 90° in comparison with T0 = 30°. The maximum difference of dimensionless deflection for the VSCL plates with SCSC boundary condition and T0 = 90° can reach to 25% lower than conventional composite plates. In addition, one may conclude from the overall results that using curvilinear fibers in composite materials which provides variable stiffness, it is possible to reach a better structural response in comparison with conventional composites, however, this should be integrated and compromised with technological and economic features of the VSCL materials. Finally, owing to the practical significance of the 3D transient analysis of VSCL plates with curvilinear fibers and also lack of information in the open literature, the presented results can be used as benchmark in the future experimental and theoretical studies. Appendix A Different degrees of Bézier curves lead to different step new multi-step formulas [45].

h Two-step formula : ynþ1 ¼ yn þ ð3f n  f n1 Þ 2 Three-step formula : ynþ1 ¼ yn þ Four-step formula : ynþ1 ¼ yn þ

4. Conclusions Based on the three-dimensional theory of elasticity, transient analysis of VSCL plates subjected to dynamic load using a hybrid differential quadrature and Bézier based multi-step method is investigated. Differential quadrature, as a powerful and stable numerical technique, together with layer-wise theory is employed to accurately model the 3D variations of the displacement components in spatial domain. Consequently, the governing partial differential equations of the problem are transform to a system of ordinary differential equations in temporal domain. Employing a new multi-step formula based on the Bézier curves provides accurate and stable solution to the resulting system of ordinary differential equations. Results revealed that predictions of the presented multi-step formulation are in excellent agreement with those of the Newmark scheme. It is also shown that time consuming of the presented method is significantly lower than Newmark method. Moreover, fast rate of convergence and high accuracy of the method are demonstrated by comparing the results with existing analytical solutions for simple problems in the open literature.

h ð9f n  8f n1 þ f n2 Þ 12

h ð175f n  81f n1 þ 15f n2  f n3 Þ 108

General m þ 1-step formula : ynþ1 " Z mþ1 # m X m m mi i ð1  uÞ u du f nmþi : ¼ yn þ h m i 1 i¼0

ðA:1-4Þ

References [1] Leissa AW, Martin AF. Vibration and buckling of rectangular composite plates with variable fiber spacing. Compos Struct 1990;14:339–57. [2] Hyer M, Lee H. The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes. Compos Struct 1991;18:239–361. [3] Hyer MW, Charette RF. Use of curvilinear fiber format in composite structure design. AIAA J 1991;29:1011–5. [4] Curry JM, Johnson ER, Starnes JH. Effect of dropped plies on the strength of graphite–epoxy laminates. AIAA J 1992;30:449–56. [5] Gürdal Z, Olmedo RA. In-plane response of laminates with spatially varying fiber orientations: variable stiffness panel concept. AIAA J 1993;31:751–8. [6] Parthasarathy VN, Kodiyalam S, Davis JE. Optimization of tow fiber paths for composite design. In: Proceedings of the 36th structural dynamics and materials conference, New Orleans; 1995, p. 1031–41.

Y. Heydarpour, M.M. Aghdam / Composite Structures 154 (2016) 344–359 [7] Tatting BF, Gürdal Z. Design and manufacture of elastically tailored tow placed plates. NASA/CR-2002-211919. [8] Abdallah MM, Setoodeh SH, Gurdal Z. Design of variable stiffness composite panels for maximum fundamental frequency using lamination parameters. Compos Struct 2007;81:283–91. [9] Honda S, Oonishi Y, Narita Y, Sasaki K. Vibration analysis of composite rectangular plates reinforced along curved lines. J Syst Des Dyn 2008;2:76–86. [10] Honda S, Narita Y, Sasaki K. Maximizing the fundamental frequency of laminated composite plates with optimally shaped curvilinear fibers. J Syst Des Dyn 2009;3:867–76. [11] Akhavan H, Ribeiro P. Natural modes of vibration of variable stiffness composite laminates with curvilinear fibers. Compos Struct 2011;93:3040–7. [12] Wu Z, Weaver PM, Raju G, Kim BC. Buckling analysis and optimization of variable angle tow composite plates. Thin-Walled Struct 2012;60:163–72. [13] Raju G, Wu Z, Kim BC, Weaver PM. Prebuckling and buckling analysis of variable angle tow plates with general boundary conditions. Compos Struct 2012;94:2961–70. [14] Ribeiro P, Akhavan H. Nonlinear vibrations of variable stiffness composite laminated plates. Compos Struct 2012;94:2424–32. [15] Honda S, Narita Y. Natural frequencies and vibration modes of laminated composite plates reinforced with arbitrary curvilinear fiber shape paths. J Sound Vib 2012;331:180–91. [16] Ribeiro P. Nonlinear free periodic vibrations of variable stiffness composite laminated plates. Nonlinear Dyn 2012;70:1535–48. [17] Houmat A. Nonlinear free vibration of laminated composite rectangular plates with curvilinear fibers. Compos Struct 2013;106:211–24. [18] Akhavan H, Ribeiro P, de Moura MFSF. Damage onset on tow-placed variable stiffness composite laminates. Compos Struct 2014;113:419–28. [19] Houmat A. Nonlinear free vibration analysis of variable stiffness symmetric skew laminates. Eur J Mech-A/Solids 2015;50:70–5. [20] Yazdani S, Ribeiro P. A layerwise p-version finite element formulation for free vibration analysis of thick composite laminates with curvilinear fibres. Compos Struct 2015;120:531–42. [21] Akhavan H, Ribeiro P, de Moura MFSF. Large deflection and stresses in variable stiffness composite laminates with curvilinear fibres. Int J Mech Sci 2013;73:14–26. [22] Yazdani S, Ribeiro P, Rodrigues JD. A p-version layerwise model for large deflection of composite plates with curvilinear fibres. Compos Struct 2014;108:181–90. [23] Yazdani S, Ribeiro P. Geometrically non-linear static analysis of unsymmetric composite plates with curvilinear fibres: p-version layerwise approach. Compos Struct 2014;118:74–85. [24] Akbarzadeh AH, Arian Nik M, Pasini D. The role of shear deformation in laminated plates with curvilinear fiber paths and embedded defects. Compos Struct 2014;118:217–27. [25] Lee YC, Reismann H. Dynamics of rectangular plates. Int J Eng Sci 1969;7:93–113. [26] Reddy JN. Dynamic (transient) analysis of layered anisotropic compositematerial plates. Int J Numer Meth Eng 1983;19:237–55. [27] Bhimaraddi A. Static and transient response of rectangular plates. Thin Walled Struct 1987;5:125–43. [28] Mallikarjuna M, Kant T. Dynamics of laminated composite plates with a higher order theory and finite element discretization. J Sound Vib 1988;126:463–75.

359

[29] Khdeir AA, Reddy JN. Exact solutions for the transient response of symmetric cross-ply laminates using a higher-order plate theory. Compos Sci Technol 1989;34:205–24. [30] Khdeir AA. Transient response of refined cross-ply laminated plates for various boundary conditions. J Acoust Soc Am 1995;97:1664–9. [31] Ray MC, Bhattacharya R, Samanta B. Exact solution for dynamic analysis of composite plates with distributed piezoelectric layers. Compos Struct 1998;66:737–43. [32] Wang YY, Lam KY, Liu GR. A strip element method for the transient analysis of symmetric laminated plates. Int J Solids Struct 2001;38:241–59. [33] Shen HS, Zheng JJ, Huang XL. Dynamic response of shear deformable laminated plates under thermomechanical loading and resting on elastic foundations. Compos Struct 2003;60:57–66. [34] Reddy JN. Mechanics of laminated composite plates and shells theory and analysis. London: CRC Press; 2004. [35] Maleki S, Tahani M, Andakhshideh A. Transient response of laminated plates with arbitrary laminations and boundary conditions under general dynamic loadings. Arch Appl Mech 2012;82:615–30. [36] Marjanovic´ M, Vuksanovic´ D, Meschke G. Geometrically nonlinear transient analysis of delaminated composite and sandwich plates using a layerwise displacement model with contact conditions. Compos Struct 2015;122:67–81. [37] Aghdam MM, Farahani MR, Dashty M, Rezaeiniya SM. Application of generalized differential quadrature method to the bending of thick laminated plates with various boundary conditions. Appl Mech Mater 2006;5–6:407–14. [38] Setoodeh AR, Tahani M, Selahi E. Hybrid layerwise-differential quadrature transient dynamic analysis of functionally graded axisymmetric cylindrical shells subjected to dynamic pressure. Compos Struct 2011;93:2663–70. [39] Tornabene F, Liverani A, Caligiana G. Laminated composite rectangular and annular plates: a GDQ solution for static analysis with a posteriori shear and normal stress recovery. Compos B Eng 2012;43:1847–72. [40] Malekzadeh P, Heydarpour Y, Golbahar Haghighi MR, Vaghefi M. Transient response of rotating laminated functionally graded cylindrical shells in thermal environment. Int J Press Vessels Pip 2012;98:43–56. [41] Heydarpour Y, Aghdam MM, Malekzadeh P. Free vibration analysis of rotating functionally graded carbon nanotube-reinforced composite truncated conical shells. Compos Struct 2014;117:187–200. [42] Heydarpour Y, Malekzadeh P, Aghdam MM. Free vibration of functionally graded truncated conical shells under internal pressure. Meccanica 2014;49:267–82. [43] Malekzadeh P, Heydarpour Y. Mixed Navier-layerwise differential quadrature three-dimensional static and free vibration analysis of functionally graded carbon nanotube reinforced composite laminated plates. Meccanica 2015;50:143–67. [44] Heydarpour Y, Aghdam MM. Transient analysis of rotating functionally graded truncated conical shells based on the Lord-Shulman model. Thin-Walled Struct 2016;104:168–84. [45] Aghdam MM, Fallah A, Haghi P. Nonlinear initial value ordinary differential equations. Nonlinear Approach Eng Appl 2015:117–36. [46] Naserian Nik AM, Tahani M. Analytical solutions for bending analysis of rectangular laminated plates with arbitrary lamination and boundary conditions. J Mech Sci Technol 2009;23:2253–67.