Dynamic analysis of functionally graded porous structures through finite element analysis

Engineering Structures 165 (2018) 287–301

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Dynamic analysis of functionally graded porous structures through finite element analysis

T

⁎

Di Wua,b, Airong Liua, , Youqin Huanga, Yonghui Huanga, Yonglin Pia,b, Wei Gaob a

Guangzhou University-Tamkang University Joint Research Centre for Engineering Structure Disaster Prevention and Control, Guangzhou University, Guangzhou 510006, China b Centre for Infrastructure Engineering and Safety, School of Civil and Environmental Engineering, The University of New South Wales, Sydney, NSW 2052, Australia

A R T I C LE I N FO

A B S T R A C T

Keywords: Functionally graded porous structures Euler-Bernoulli beam Timoshenko beam Dynamic analysis Finite element analysis

A finite element method (FEM) analysis framework is introduced for the free and forced vibration analyses of functionally graded porous (FGP) beam type structures. Within the proposed computational scheme, both EulerBernoulli and Timoshenko beam theories have been adopted such that the explicit stiffness and mass matrices for 2-D FGP beam element through both beam theories are explicitly expressed. Both Young’s modulus and material density of the FGP beam element are simultaneously considered as grading through the thickness of the beam. The material constitutive law of a FGP beam is governed by the typical open-cell metal foam. Furthermore, the damping effects of the FGP structures can be also incorporated within the proposed FEM analysis framework through the Rayleigh damping model. Consequently, the proposed approach establishes a more unified analysis framework which can investigate simple FGP beams as well as complex FGP structural systems involving mixture of different materials. In order to demonstrate the applicability, accuracy, as well as the efficiency of the proposed computational scheme, both FGP beams and frame structures with multiple porosities have been rigorously explored.

1. Introduction From the continuous development of modern society, there has always been a high demand for new building materials for various innovative engineering structures. As a general trend, high cost-effective materials are more desirable than ever before, so safer, cheaper, and more importantly, greener engineering projects can be anticipated [1–3]. Among numerous ingeniously invented engineering materials, the functionally graded porous material (FGPM) has become one of the research foci across many engineering disciplines [4–6]. As one of the advanced new generation composite materials, the FGPM has inherited the advantages from both traditional functionally graded materials and metallic porous materials by having attractive characteristics such as superb stiffness-to-weight ratio, energy dissipation, mechanical damping, as well as designable vibrational frequency [7]. Examples of the implementation of FGPMs in engineering applications are including sandwich panels with steel foam cores, balcony platforms, parking floor slabs, and race car crash absorbers, etc. [7]. Consequently, the FGPM has offered an opportunity for next-generation structures (i.e., at both macro-and micro-scales) possessing robust level of safety but with

⁎

lighter weight and less carbon footage. Since the emerging of such engineered composite material, the research work on FGPM has attracted numerous global research attentions. In 1964, Biot [8] proposed a poroelasticity theory for buckling analysis of fluid-saturated porous slab through the framework of thermodynamics of irreversible processes. Magnucki and Stasiewicz [9] investigated the elastic buckling problem for porous beams subjected to static compressive forces through an energy approach. MagnuckaBlandzi [10] studied the dynamic buckling problem for a metal foam circular plate through an analytical approach by assuming the middle plane is being symmetric. Jasion et al. [11], assessed the static stability of sandwich plates with metal foam cores. Both global and local buckling issues of the sandwich plates were investigated through the analytical and numerical approaches, and the results on the critical buckling loads were also compared with experiments for plates with various boundary conditions. Moreover, the buckling analysis of cylindrical shells made of porous materials subjected to the combination of axial loads and external pressure was investigated by Belica and Magnucki [12]. The geometric nonlinearity was also incorporated and the static critical buckling load was determined through an analytical analysis framework. Jabbari et al. [13] proposed a buckling analysis for

Corresponding author. E-mail address: [email protected] (A. Liu).

https://doi.org/10.1016/j.engstruct.2018.03.023 Received 30 August 2017; Received in revised form 25 February 2018; Accepted 11 March 2018 0141-0296/ © 2018 Elsevier Ltd. All rights reserved.

Engineering Structures 165 (2018) 287–301

D. Wu et al.

porous circular plate with piezoelectric layers. Subsequently, a new analytical approach was proposed in [14] for the static buckling analysis of porous soft ferromagnetic FG circular plates under the influence of magnetic fields. In addition to the buckling analysis of plates made of various porous materials, Chen et al. [15–17] investigated the problem of static bending, elastic buckling, free vibration, forced vibration, as well as nonlinear free vibration for FGP beams with the consideration of shear effects through an analytical analysis framework. Furthermore, Kitipornchai et al. [18] proposed an analytical approach for free vibration and static buckling analyses for FGP beams with nano-graphene platelet reinforcements. Soon after that, Chen et al. [19] extended their analytical approach to investigate the problems of nonlinear vibration and postbuckling for graphene reinforced FGP beams. Even though there are some studies have been reported on various structural analyses for FGP beams, it is noticed that all aforementioned works regarding the dynamic analysis (free and forced vibration) of FGP beams have been proposed through the analytical analysis framework. Moreover, another common feature that can be realized from the previous studies is that only FGP beams with various boundary conditions subjected to different loading regimes have been investigated. However, for real-life engineering applications, complex structural systems (e.g., frames and trusses) are often being implemented to serve various engineering purposes. Therefore, it is essential to develop other types of analysis framework which should have the same level of accuracy as the analytical approach does, but with more extensive applicability. In this study, a finite element method (FEM) analysis framework is proposed for free and forced vibration analyses of FGP beam type structures. Within the proposed computational scheme, both EulerBernoulli and Timoshenko beam theories have been adopted such that the explicit stiffness and mass matrices for 2-D FGP beam element through both beam theories are explicitly expressed. In addition to the capability on performing valid dynamic analysis for FGP beams, the proposed approach is also applicable to FGP structural systems (e.g., frames and trusses) with multiple types of functionally graded materials. Also, the damping effects of the FGP structural systems have been incorporated within the proposed computational scheme through the implementation of the Rayleigh damping model. Consequently, a more unified analysis framework is proposed in this study which rigorously maintains the accuracy level of analytical approaches, but at same time possesses the freedom of solving more complex structural systems involving diversity of engineered composite materials. In order to achieve a more effective presentation, this paper is organized as follows. Section 2 introduces the proposed FEM for static analysis of FGP frame structures. Subsequently, the dynamic analysis of FGP frame structures through FEM is presented in Section 3. In order to demonstrate the applicability, accuracy and efficiency of the proposed computational scheme, three distinctive numerical examples have been thoroughly investigated in Section 4. Finally, the conclusion of this investigation is drawn in Section 5.

( ) ⎤⎦ + E cos ( ) = E ⎡⎣1−e cos ( ) ⎤⎦ ( ) ⎤⎦ + G cos ( ) = G ⎡⎣1−e cos ( ) ⎤⎦ ( ) ⎤⎦ + ρ cos ( ) = ρ ⎡⎣1−ρ cos ( ) ⎤⎦

⎧ E i (z ) = E1i ⎡1−cos ⎣ ⎪ ⎪ i i SMCR G (z ) = G1 ⎡1−cos ⎨ ⎣ ⎪ i i ⎪ ρ (z ) = ρ1 ⎡1−cos ⎣ ⎩

πz hi

i 0

πz hi

i 1

i 0

πz hi

πz hi

i 0

πz hi

i 1

i 0

πz hi

πz hi

πz hi

i 0

i 1

i m

πz hi

(1)

and the MMCR model of the ith FGP beam element can be expressed as [15–17]:

(

)

(

)

(

)

⎧ E i (z ) = E1i ⎡1−cos πzi + π ⎤ + E0i cos πzi + π 4 4 ⎦ 2h 2h ⎣ ⎪ ⎪ πz π = E1i ⎡1−e0i cos i + 4 ⎤ ⎪ 2h ⎣ ⎦ ⎪ ⎪ Gi (z ) = G i ⎡1−cos πz + π ⎤ + G i cos πz + π 1 0 ⎪ 4 ⎦ 4 2hi 2hi ⎣ MMCR ⎨ πz π = G1i ⎡1−e0i cos i + 4 ⎤ ⎪ 2h ⎣ ⎦ ⎪ πz i ⎪ ρi (z ) = ρ ⎡1−cos i + π ⎤ + ρ i cos πzi + π 1 0 4 ⎦ 4 2h 2h ⎣ ⎪ ⎪ πz π i i = ρ1 ⎡1−emcos i + 4 ⎤ ⎪ 2h ⎣ ⎦ ⎩

(

(

)

)

(

(

)

(

E1i

)

(

)

)

(2)

E0i

where and are the maximum and minimum Young’s moduli of the ith FGP beam element, respectively; G1i = E1i/(2 + 2νi ) and G0i = E0i /(2 + 2νi ) are the maximum and minimum shear moduli of the ith FGP beam element, respectively; νi denotes the Poisson’s ratio of the ith FGP beam element; Without loss of generality, the Poison’s ratio (i.e., νi ) of the ith FGP beam element is assumed to be a constant [15–17]. ρ1i and ρ0i are the maximum and minimum densities of the ith FGP beam element, respectively; e0i = 1−E0i / E1i and emi = 1−ρ0i / ρ1i are the porosity coefficients of the relative Young’s modulus and density of the ith FGP beam element which are also governed by the relationship of an open-cell foam as [20,21]: 2

ρi E0i = ⎜⎛ 0i ⎟⎞ ⇒ emi = 1− 1−e0i i E1 ⎝ ρ1 ⎠ hi

hi

(3)

and − 2 ⩽ z ⩽ 2 , hi denotes the thickness of the ith FGP beam element. In order to achieve a more effective illustration of the effects of different porosity distributions, the variation of the Young’s modulus of a generic FGP beam with a rectangular cross-section under two different material constitutive relationships are depicted in Fig. 1. Within the SMCR model, the material properties of the FGP beam can have a symmetric profile about the middle plane of the thickness direction. The maximum material properties are being located at the two extreme layers (e.g., top and bottom layers) of the beam and the minimum material properties located at the neutral axis of the FGP beam. On the other hand, the MMCR model is governing the material properties in a monotonic increasing fashion from the bottom layer to the top layer. That is, within the MMCR model, the minimum material properties are located at the bottom layer and the maximum properties are graded at the top layer of the beam. Consequently, such capability

2. Static analysis of FGP frame structures through FEM 2.1. Material models for FGPM Within the scope of this study, it is presumed that the material properties (i.e., Young’s modulus, shear modulus, and material density) of a generic FGP beam element with rectangular cross section are continuously varying along the beam thickness direction. In order to achieve a more generalized analysis framework, two distinctive constitutive laws of the FGPM are considered which are including the symmetric and monotonic material constitutive relationships (i.e., SMCR and MMCR). Specifically, the SMCR model of the ith FGP beam element can be expressed as [15–17]:

Fig. 1. The Young’s modulus of FGP beam with (a) SMCR and (b) MMCR.

288

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Fig. 2. Geometric layout of a typical FGP beam element.

Eq. (6), the stiffness matrix of the ith FGP beam element (i.e., KiT ) can be explicitly formulated as:

on manipulating the material properties along the thickness direction of the beam has certainly provides additional design freedom for engineers to satisfy real-life situationally dependent engineering design criteria.

KiT = Kiε + K iγ

(7)

where 2.2. Static analysis of FGP frame structures based on the Timoshenko beam model

i K εi,13 K εi,14 0 0 ⎡ K ε,11 ⎢ i K ε,22 K εi,23 0 K εi,25 ⎢ ⎢ K εi,33 K εi,34 K εi,35 Kiε = ⎢ ⎢ K εi,44 0 ⎢ Symmetric K εi,55 ⎢ ⎢ ⎣

The general structural layout of a generic FGP beam element i is illustrated in Fig. 2. Within the linear analysis framework of FEM, the strain field of a FGP beam i based on the Timoshenko beam theory can be formulated as: ∂u0i

i ⎧ ⎪ εxx = ⎨ i ⎪ γxz = ⎩

∂x

−z

∂w0i ∂x

∂ϕxi ∂x

−ϕxi

0 0 ⎡0 K γi ,22 K γi ,23 ⎢ ⎢ K γi ,33 ⎢ K iγ = ⎢ ⎢ ⎢ Symmetric ⎢ ⎣

(4)

i where and γxz denote the normal and shear strains of the ith FGP beam element; u0i and w0i denote the displacements of a generic point on i εxx

the neutral axis of the ith FGP beam element in x- and z-directions, respectively; ϕxi denotes the bending rotation of the cross-section at any point on the neutral axis. Without loss of generality, combinations of Lagrange and Hermite interpolation polynomials have been implemented for the displacement field [22–24]. From the Hooke’s law, the stress field of the ith FGP beam element can be formulated as: i i i ⎧ σxx = E (z ) εxx i i ⎨ τxz = Ksi Gi (z ) γxz ⎩ i σxx

The components of

=

1 2

(5)

∫0 ∫A σxxi εxxi dAidx + 12 ∫0 ∫A τxzi γxzi dAidx Li

i

∫0 ∫

+

Ksi 2

Li

i

2 i 2 i ⎛ ∂u0 −z ∂ φx ⎞ dAidx i ⎟ i E (z ) ⎜ A ∂x ∂x 2

1 2

⎝

∫0 ∫A

i

⎠

1 2

+

Li

∫0 ∫A

i

Ksi 2

{

∂u0i

∂2φxi

∂x

∂x 2

}

i

{

∂w0i ∂x

(10)

∂w Gi (z ) − Gi (z ) ⎤ ⎧ ∂x0 ⎫ i dA dx φxi ⎡ ⎢− Gi (z ) Gi (z ) ⎥ ⎨ i ⎬ ⎣ ⎦ ⎩ φx ⎭

(11)

Consequently, the potential energy of the ith FGP beam based on the Euler-Bernoulli beam theory can be determined as:

⎧ ∂u0 ⎫ i i ⎡ E (z ) − zE (z ) ⎤ ⎪ ∂x ⎪ dAidx ⎢− zE i (z ) z 2E i (z ) ⎥ ⎨ ∂2φxi ⎬ ⎣ ⎦⎪ 2 ⎪ ⎩ ∂x ⎭

}

are presented in Appendix A.

i i σxx = E i (z ) εxx

Ui =

i

Li

∫0 ∫A

(9)

and the normal stress can be formulated as:

∂w i Gi (z ) ⎛⎜ 0 −φxi ⎞⎟ dAidx ⎠ ⎝ ∂x i

=

and

K iγ

∂u0i ∂ 2w i −z 20 ∂x ∂x

i εxx =

2

Li

0 0 0 ⎤ 0 K γi ,25 K γi ,26 ⎥ ⎥ 0 K γi ,35 K γi ,36 ⎥ 0 0 0 ⎥ ⎥ K γi ,55 K γi ,56 ⎥ ⎥ K γi ,66 ⎦

Unlike the Timoshenko beam model, the shear effects have been neglected within the framework of the Euler-Bernoulli beam theory. Even though the rotation considered in the Euler-Bernoulli beam is only from pure bending, the applicability of such beam model is also extensive in rea-life engineering applications, especially for beams having large span-to-width ratio. Within the framework of the Euler-Bernoulli beam theory, the strain field only contains the normal component due to pure bending effect, which is:

i τxz

Li

(8)

2.3. Static analysis of FGP frame structures based on the Euler-Bernoulli beam model

where and denote the normal and shear stresses, respectively; Ksi denotes the shear correction factor (e.g., Ksi = 5/6 for a rectangular cross-section). Therefore, the potential energy of the ith FGP beam element based on the Timoshenko beam model is:

Ui =

Kiε

K εi,16 ⎤ ⎥ K εi,26 ⎥ i ⎥ K ε,36 ⎥ K εi,46 ⎥ ⎥ K εi,56 ⎥ ⎥ K εi,66 ⎦

1 2

L

∂u i

L

∫0 i ∫Ai σxxi εxxi dAidx = 12 ∫0 i ∫Ai E i (z ) ⎛ ∂x0 −z ⎝

∂2w0i ∂x 2

2

⎞ dAidx ⎠

i

1 2

Li

= ∫0 ∫Ai

(6)

Furthermore, by employing the Castigliano’s theorem [22,25,26] to 289

{

∂u0i

∂2w0i

∂x

∂x 2

}

⎧ ∂u0 ⎫ i i ⎡ E (z ) −zE (z ) ⎤ ⎪ ∂x ⎪ dAidx ⎢−zE i (z ) z 2E i (z ) ⎥ ⎨ ∂2w0i ⎬ ⎣ ⎦⎪ 2 ⎪ ⎩ ∂x ⎭

(12)

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Once again, by implementing the Castigliano’s theorem [22,25,26] to Eq. (12), the stiffness matrix of the ith Euler-Bernoulli type FGP beam element can be explicitly formulated. It is highlighted that since the Timoshenko beam theory can be degenerated into the Euler-Bernoulli beam theory, thus the explicit formulation of the stiffness matrix of the ith Euler-Bernoulli type FGP beam element (i.e., KiEB ) can be expressed as Eq. (7) with Φi = 0 . That is

KiEB = KiT (Φi = 0)

i i i M(1,2) M(1,3) ⎡ M(1,1) ⎢ i i M(2,2) M(2,3) ⎢ i ⎢ M(3,3) MiT = ⎢ ⎢ ⎢ symmetric ⎢ ⎢ ⎢ ⎣

(13)

The components of

After the establishment of the stiffness matrix for each presented FGP beam element within the structural system, the global stiffness matrix of the FGP frame structure can be calculated as:

MiT

i i i M(1,4) M(1,5) M(1,6) ⎤ ⎥ i i i M(2,4) M(2,5) M(2,6) ⎥ i i i ⎥ M(3,4) M(3,5) M(3,6) ⎥ i i i ⎥ M(4,4) M(4,5) M(4,6) ⎥ i i M(5,5) M(5,6) ⎥ ⎥ i M(6,6) ⎥ ⎦

(18)

are presented in Appendix B.

3.2. Dynamic analysis of FGP frame structures based on the Euler-Bernoulli beam model

n

K=

∑

GTi K i(•) Gi

When the Euler-Bernoulli beam model is considered, the velocity field of the ith FGP beam element can be expressed as:

(14)

i=1

where K ∈ R d × d denotes the global stiffness matrix; d denotes the total number of degrees of freedom; n denotes the total number of beam element presented in the system; and Gi denotes the transformation matrix of element i; K i(•) denotes the elemental stiffness matrix formulated by the ‘ (•) ’ beam theory. Therefore, the static analysis of FGP frame structure can be conquered by solving the governing linear system of equations, that is

i

u̇ i

1 2

Ke =

where u ∈ R d denotes the structural response vector; and F ∈ R d denotes the applied load vector that containing both mechanical and thermal contributions. It is also emphasized that the proposed approach for the static analysis of FGP frame structures are also applicable to FGP trusses. By simply considering the axial effects only, the stiffness matrices presented in Eqs. (7) and (13) can be convergent to an identical matrix which can be served as the stiffness matrix of a FGP truss element. Therefore, the proposed computational approach possesses extensive applicability by being competent to encompass multiple structural types within a single analysis framework.

L

∫0 i ∫Ai ρi (z )·(u̇ i)T ·u̇ idAidx

{

L

1

= 2 ∫0 i ∫Ai u̇0i ẇ0i

∂ẇ 0i ∂x

}

i i 0 −zρi (z ) ⎤ ⎧ u̇0 ⎫ ⎡ ρ (z ) ⎪ i⎪ ⎢ 0 0 ⎥ ẇ0 dAidx ρi (z ) ⎥ ⎢ ⎨ ⎬ ∂ẇ 0i i 0 z 2ρi (z ) ⎥ ⎢ ⎪ ⎦⎪ ⎣−zρ (z ) ⎩ ∂x ⎭

(20) Once more, by utilizing the Lagrange interpolation polynomials as shape functions for the axial effects, and the cubic Hermite interpolation polynomials for the bending effects, the mass matrix of the ith Euler-Bernoulli type FGP beam element can be explicitly derived. Since the Euler-Bernoulli beam model can be considered as a special case of the Timoshenko beam theory, consequently, the explicit formulation of the mass matrix of the ith Euler-Bernoulli type FGP beam element can be determined from Eq. (18) by setting Φi = 0 . That is,

3. Dynamic analysis of FGP frame structures through FEM 3.1. Dynamic analysis of FGP frame structures based on the Timoshenko beam model

MiEB = MiT (Φi = 0)

Based on the static analysis presented in Section 2, the velocity field, u̇ i , of a Timoshenko type FGP beam element can be explicitly formulated as:

3.3. Free and forced vibration analyses of FGP frame structures

⎧ ∂u0 −z· ∂ϕx ⎪ ⎫ ⎪ u ̇i = ⎧ i ⎫ = ∂t i ∂t ∂ w ⎨ ⎬ ⎨ ⎬ 0 ⎩ ẇ ⎭ ⎪ ⎪ ∂t ⎭ ⎩

(16)

n

M=

where u̇i denotes the velocity in the x-direction; ẇ i denotes the velocity in the z-direction; and t denotes the time variable. With the consideration of the rotary inertia and the axial inertial effects, the kinetic energy of the ith Timoshenko type FGP beam can be determined as:

Ke =

1 2

L

= ∫0 ∫Ai

{

u̇0i

ẇ0i

i ϕẋ

}

i 0 −zρi (z ) ⎤ ⎧ u̇0 ⎫ ⎡ ρ (z ) ⎪ i⎪ i ⎢ 0 0 ⎥ ẇ0 dAidx ρ (z ) ⎥ ⎢ ⎨ i⎬ i 0 z 2ρi (z ) ⎥ ⎢ ϕ̇ ⎪ ⎦⎪ ⎣−zρ (z ) ⎩ x⎭

GTi Mi(•) Gi

(22)

where M ∈ R d × d denotes the global mass matrix of the FGP frame structure which incorporates both the rotary inertia and axial inertial effects. Since the global stiffness and mass matrices of the FGP frame structure are available, the free vibration analysis of the concerned structural system can be explicitly formulated as an eigen-analysis, that is:

∫0 i ∫Ai ρi (z )·(u̇ i)T ·u̇ idAidx Li

∑ i=1

i

1 2

(21)

After the competent establishment on both the stiffness and mass matrices of the FGP beam element within the two adopted beam theories, the global mass matrix of the FGP frame structure involving d degrees of freedom and n FGP beam elements can be calculated as:

i

i

u̇ i

(19)

The kinetic energy of the ith FGP beam element, which incorporates both the rotary inertia and axial inertia effects, can be calculated as:

(15)

Ku = F

i

⎧ ∂u0 − ∂ ⎛z· ∂w0 ⎞⎫ ⎪ ∂t ∂t ⎪ u ̇i ⎝ ∂x ⎠ = ⎧ i⎫ = ⎨ ⎬ ∂w0i ⎩ ẇ ⎬ ⎭ ⎨ ⎪ ⎪ ∂t ⎩ ⎭

(K−ωj2 M) Φj = 0, for j = 1,…,d (17)

ωj2

(23)

where ∈ R + denotes the jth eigenvalue of Eq. (23) and ωj ∈ R + denotes the natural frequency of the jth mode; Φj ∈ R d denotes the jth eigenvector or mode shape that is corresponding to ωj . In addition to the free vibration analysis of the FGP frame structures, the forced vibration can also be conducted within the proposed analysis framework. According to Newton’s second law of motion, the

i

where u̇0i = ∂u0i / ∂t , ẇ0i = ∂w0i / ∂t , and ϕẋ = ∂ϕxi / ∂t . Once again, by implementing a combination of Lagrange and Hermite interpolation polynomials, the mass matrix of the ith Timoshenko type FGP beam element (i.e., MiT ) can be explicitly formulated as: 290

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Table 1 Dimensionless natural frequency of the first mode of an aluminium beam with SS boundary conditions. Dimensionless natural frequency of the first mode λ1 = L/h

ω1 L2 h

Proposed method

10 30 100

ρA EA

Ref. [31]

TBT

CBT

TBT

CBT

TBT

CBT

2.8042 2.8440 2.8486

2.8375 2.8478 2.8490

2.8023 2.8439 2.8496

2.8375 2.8478 2.8496

2.797 2.843 2.848

2.849 2.849 2.849

based approach is competent to simultaneously incorporate various material types (i.e., laminated composites, functionally graded materials, pure metallic materials etc.) within a single analysis framework. Consequently, the proposed computational method is more adequate for practical engineering applications by possessing the ability to handle both simple FGP beams as well as complex FGP truss and frame structures.

governing equation of the forced vibration of the FGP structural system with d degrees of freedom can be formulated as:

Mu ¨ (t ) + Cu̇ (t ) + Ku (t ) = F (t )

(24)

Rd×d

Rd

¨ ,u̇ ∈ denotes the damping matrix; u,u denote the where C ∈ FGP structural displacement, velocity and acceleration vectors, respectively; F(t ) denotes the applied dynamic loading functions. For the purpose of this paper, the Rayleigh damping model is adopted such that:

4. Numerical investigations

(25)

C = a0 M + a1 K

4.1. Result validation through free vibration analysis

where a 0 and a1 are the proportionality constants. Generally, there are many time-stepping algorithms that can be implemented to adequately solve the second order ordinary differential equation formulated in Eq. (24). Within the scope of this paper, the prevalently employed Newmark’s method [27–30] is adopted within the proposed analysis framework. For the purpose of achieving the completeness of this study, the algorithm of the Newmark’s scheme is briefly outlined as follows. Within the scheme of Newmark’s method, the constant average acceleration method is adopted with β = 0.25, and γ = 0.5 for all subsequent numerical investigations.

In order to demonstrate the accuracy of the proposed method, a substantial amount of verifications have been conducted. Three distinctive investigations have been performed through the free vibration analysis of beams with diverse materials. For the first result verification, a simply supported (SS) aluminium beam is investigated for the free vibration analysis by utilizing the Timoshenko beam theory (TBT) and Euler-Bernoulli or classical beam theory (CBT). Identical problems have also been investigated in [31,32]. Regarding the pure aluminium beam, it has a Young’s modulus of EA = 70 GPa , and a density of ρA = 2700 kg/m3 . For this particular exercise, a square cross-section is adopted with b = h = 0.1 m . By implementing the proposed method, the dimensionless natural frequencies of the first model of the aluminium beam are calculated and then reported in Table 1 for various L/h ratios. Also, the calculated natural frequencies from the proposed method are compared with the results provided in [31,32]. For the second result verification, a functionally graded (FG) beam manufactured with alumina (Al2O3) top surface and aluminium bottom surface is investigated. For this particular model, the power-law material constitutive relationship with k = 0.3, where k denotes the exponent of the power-law, is adopted. In particular, the top alumina has material properties as: Et = 380 GPa , ρt = 3800 kg/m3 , υt = 0.23 [31]; and the bottom aluminium has material properties as: Eb = 70 GPa , ρb = 2700 kg/m3 , υt = 0.23 [31]. The cross-section has dimensions b = 0.1 m and h = 0.125 m . Once again, a free vibration through TBT is conducted, and the corresponding dimensionless fundamental natural frequencies of the functionally graded beam with various boundary conditions are summarized in Table 2. In the second verification exercise, identical problem has also been investigated in [31–33], and the results provided in [31–33] are also summarized in Table 2. In this particular case, both analytical [31,32] and experimental [33] results are available for the purpose of result verification. In this particular validation, the considered boundary conditions are including simply supported-simply supported (SS), clamped–clamped (CC), clampedsimply supported (CS) and clamped-free (CF). For the last validation, a simply supported (SS) functionally graded beam with alumina at top and steel at the bottom is considered. The material is also governed by the power-law with k = 1. The material properties of the alumina at the top layer are Et = 390 GPa , ρt = 3960 kg/m3 , υt = 0.25 [31]; and the material properties of the steel at the bottom layer are Eb = 210 GPa , ρb = 7800 kg/m3, υt = 0.31 [31].

Algorithm. Forced Vibration through Newmark’s method Start 1. Initial calculation (step i = 0) ¨ 0 = M−1 (F0−Cu̇ 0−Ku 0) ; 1.1. u 1.2. Select time step: Δt ; γ 1 C+ 1.3. K ̂ = K + 2 M; β Δt

1.4. a =

1 M β Δt

Ref. [32]

+

β (Δt )

γ C, β

b=

1 M 2β

+ Δt

(

) C;

γ −1 2β

2. Calculate results for each time step i = 1,2,…,m 2.1. ΔFi ̂ = ΔFi + au̇ i + bu ¨ i; −1 2.2. Δui = K ̂ ΔFi ;̂ 2.3. Update and save the displacement, velocity and acceleration:

(

Δu̇ i =

γ γ Δui− β u̇ i β Δt

Δu ¨i =

1 1 1 Δui− β Δt u̇ i− 2β u ¨ i; β (Δt )2

γ

)

+ Δt 1− 2β u ¨ i;

ui + 1 = ui + Δui , u̇ i + 1 = u̇ i + Δu̇ i , u ¨ i+1 = u ¨ i + Δu ¨i 3. Replace i by i + 1 and repeat steps 2.1–2.3 for the next step. End

3.4. Remarks The finite element method for static and dynamic analyses of FGP frame structures presented in Sections 2 and 3 provides an equally accurate, but more generalized, analysis framework when comparing with the analytical approaches. Even though FGP beam element has been illustrated in Sections 2 and 3, the proposed numerical approach can also be extended to FGP trusses and grid-like structures. Moreover, another advantage associated with the proposed method is that the FEM 291

Engineering Structures 165 (2018) 287–301

D. Wu et al.

with two distinctive material models through both the Timoshenko and Euler-Bernoulli beam theories are calculated, and the results are summarized in Table 4. In addition, the results of the FGP beam with CF boundary condition are summarized in Table 5. From Tables 4 and 5, it can be observed that the FGP beam with symmetric material constitutive relationship generally possesses higher natural frequencies than the beam with monotonic material constitutive relationship for both the Timoshenko and Euler-Bernoulli beam theories. Otherwise, the structural natural frequencies determined based on the Timoshenko beam theory are generally lower than the magnitudes of the natural frequencies of the beam with identical L/h ratios obtained from the Euler-Bernoulli beam theory.

Table 2 Dimensionless natural frequency of the first mode of an FG beam with L/h = 10, k = 0.3 against various boundary conditions. Dimensionless natural frequency of the first mode

λ1 =

ω1 L2 h

∫ ρ (z ) dA ∫ E (z ) dA

Boundary conditions

Proposed method

Ref. [31]

Ref. [32]

Ref. [33]

SS CC CF CS

2.7367 5.8691 0.9701 4.1485

2.7450 5.9544 0.9858 4.2030

2.774 6.013 0.996 –

2.803 6.078 1.008 4.291

4.2. Forced vibration of FGP beam subjected to concentrated moving harmonic loads

Similar to the first two verifications, the dimensionless fundamental natural frequencies of the functionally graded beam with various L/h ratios are calculated and then presented in Table 3. Identical problem has been investigated in [31,34], whose results are also reported in Table 3. As illustrated in Tables 1–3, excellent agreements between proposed method and reported studies can be clearly observed. Within the four distinctively reported works, both analytical and experimental results were adopted. Through the rigorous cross-validation, the accuracy of the proposed method has been evidently demonstrated. After the thorough validation of the proposed method, the free vibration analysis for FGP beams with two different material models (i.e., SMCR and MMCR) subjected to various boundary conditions are also conducted. For this particular part of investigation, two types of boundary conditions, namely clamped-clamped and clamped-free, are separately examined. The material properties of the FGP beam are E1 = 200 GPa , ρt = 7850 kg/m3 , υt = 0.33, and e0 = 0.5. By utilizing the proposed method with discretization of 100 FGP beam elements, the first five modes of dimensionless natural frequencies of the FGP beam

In prior to the formal investigation on the forced vibration of a simply supported FGP beam subjected to concentrated moving harmonic loads, one additional validation of the proposed method for forced vibration involving a concentrated moving load is investigated in the first part of Section 4.2. In order to further verify the accuracy of the proposed analysis framework, one additional benchmark example, which was firstly investigated in [35], is re-visited here by utilizing the proposed FEM approach. The general structural layout is illustrated in Fig. 3. For the benchmark example, a simply supported alumina-steel FG beam, which is governed by the power-law, is analysed. The material properties of the FG beam are adopted from [35], which are Et = 390 GPa , ρt = 3960 kg/m3, Eb = 210 GPa , ρb = 7800 kg/m3. In this particular case, the forced vibration is executed within the Euler-Bernoulli beam theory. The dimensions of the simply supported FG beam are: b = 0.4 m , h = 0.9 m , and L = 20 m . The applied moving load is:

Table 3 Dimensionless natural frequency of the first mode of an FG beam with SS boundary conditions, and k = 1. Dimensionless natural frequency of the first mode λ1 = 100ω1 h Proposed method

ρb Eb

Ref. [31]

Ref. [34]

L/h

CBT

TBT

CBT

TBT

CBT

TBT

100 10 5

0.039105 3.8945 15.388

0.039100 3.8476 14.715

0.039218 3.9059 15.436

0.039213 3.8586 14.756

0.039219 3.9060 15.436

0.039215 3.8663 14.861

Table 4 Dimensionless natural frequency of FGP beam with CC boundary conditions. Dimensionless natural frequency of FGP beam (CC) λ i = L/h = 5

ωi L2 h

ρ1 E1

L/h = 20

L/h = 50

Mode

SMCR

MMCR

SMCR

MMCR

SMCR

MMCR

CBT

λ1 λ2 λ3 λ4 λ5

6.3393 14.3794 16.5216 28.7624 30.0004

5.7687 14.3658 15.1039 27.4569 28.8300

6.4716 17.7708 34.6311 56.7868 57.5187

5.8807 16.1529 31.4917 51.6680 57.5170

6.4792 17.8492 34.9578 57.7105 86.0649

5.8872 16.2191 31.7679 52.4504 78.2314

TBT

λ1 λ2 λ3 λ4 λ5

5.0185 11.2724 14.3794 18.6110 26.4276

4.7216 10.7878 14.3780 17.9640 25.6576

6.3476 17.0542 32.3755 51.5447 57.5178

5.7872 15.6088 29.7656 47.6214 57.5161

6.4588 17.7265 34.5502 56.6999 83.9717

5.8718 16.1267 31.4604 51.6863 76.6454

292

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Table 5 Dimensionless natural frequency of FGP beam with CF boundary conditions. Dimensionless natural frequency of FGP beam (CF) λ i = L/h = 5

ωi L2 h

ρ1 E1

L/h = 20

L/h = 50

Mode

SMCR

MMCR

SMCR

MMCR

SMCR

MMCR

CBT

λ1 λ2 λ3 λ4 λ5

1.0099 6.0333 7.1895 15.7814 21.5703

0.9181 5.5005 7.1947 14.4465 21.5194

1.0179 6.3590 17.7152 28.7581 34.4614

0.9249 5.7793 16.1062 28.7568 31.3503

1.0184 6.3788 17.8461 34.9293 57.6512

0.9253 5.7960 16.2168 31.7436 52.3996

TBT

λ1 λ2 λ3 λ4 λ5

0.9819 5.1639 7.1895 12.0916 19.8225

0.8969 4.8188 7.1819 11.4543 18.9968

1.0160 6.2788 17.2051 28.7580 32.7317

0.9235 5.7189 15.7196 28.7544 30.0317

1.0181 6.3656 17.7589 34.6175 56.8410

0.9251 5.7861 16.1512 31.5084 51.7874

Fig. 3. Simply supported FGP beam subjected to a moving harmonic concentrated load.

⎧Q (t ) = Q0 L L L ⎨ x Q (t ) = v·t − 2 ,− 2 ⩽ x Q ⩽ 2 ,t1 = 0 ⩽ t ⩽ t2 = ⎩

L v

MMCR) with various porosity coefficients are explicitly investigated. In the first parametric study, a moving concentrated load with constant magnitude is considered. Also, the normalized central deflection of the FGP beam (w0 (x = 0,t )/ Ds ) has been selected once again for the observation reference. Consequently, by employing the proposed method, the normalized central FGP beam deflections are calculated and the results are plotted in Fig. 4. As demonstrated in Fig. 4, when the porosity coefficient is increased, the normalized central deflection would also be increased regardless of the adopted material constitutive relationships. However, by comparing two different material constitutive relationships, the symmetric material constitutive relationship provides better performance than the FGP beam with monotonic material constitutive relationship by experiencing less structural deformation. In addition, the normalized central deflection calculated by the Timoshenko beam theory is more severe than the ones calculated from the Euler-Bernoulli beam theory in both the SMCR and MMCR cases. Since the Timoshenko beam theory is more critical than the EulerBernoulli model for the forced vibration analysis of the FGP beam, therefore only the Timoshenko beam theory is adopted for the second parametric study. In the second analysis, a moving harmonic concentrated load is applied to the FGP beam structural system as outlined

(26)

where Q (t ) denotes the applied concentrated load; Q0 denotes the magnitude of the applied load; x Q (t ) denotes the location of the moving load; v denotes the velocity of the moving load. Also, in order to compare results with the ones reported in [35], the maximum central deflection of the FG beam is selected, that is, max(w0 (x = 0,t )/ Ds ) , with

Ds =

Q0 L3 48Esteel I

(27)

where Ds denotes the central deflection of the steel beam with dimensions b = 0.4 m , h = 0.9 m , and L = 20 m subjected to static load Q0 being applied at the mid-span; I denotes the second moment of area of the steel beam. Consequently, by implementing the proposed FEM method with discretization size of 500, the maximum normalized central deflection of the FG beam can be calculated, and then reported in Table 6 for various power-law exponents and velocities. From Table 6, excellent agreement can be clearly observed for all calculated cases when comparing with the results reported in [35]. Therefore, the accuracy of the proposed analysis framework has been further verified for forced vibration analyses. Since the accuracy of the proposed method has been satisfactorily verified by the investigation of FG beam, subsequently, the proposed method is implemented to conduct the forced vibration analysis of a simply supported FGP beam subjected to a moving harmonic concentrated load. Regarding the FGP beam, the material properties are E1 = 200 GPa , ρ1 = 7850 kg/m3 , and υ = 1/3. The dimensions of the FGP beam presented in Fig. 3 are b = 0.4 m , h = 0.9 m , and L = 20 m . Both symmetric and monotonic material constitutive relationships (i.e., SMCR and

Table 6 Maximum normalized structural central deflection at specific velocity.

293

Power-law Exponent k

Velocity of the load (m/ s)

Max (w0 (x = 0,t )/ D ) [35]

Max (w0 (x = 0,t )/ D ) Proposed method

0.2 0.5 1 2

222 198 179 164

1.0344 1.1444 1.2503 1.3376

1.0348 1.1446 1.2504 1.3378

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Fig. 4. Normalized central deflection of the simply supported FGP beam subjected to moving load with constant magnitude.

Fig. 5. Parametric study on the effects of velocity and excitation frequency of the moving harmonic load acting on the structural dynamic responses of the FGP beam.

294

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Fig. 6. Parametric study on the effects of the porosity coefficients of the FGP beam acting on the structural dynamic responses of the FGP beam.

responses, a third parametric study is also conducted to investigate the relationship between the porosity coefficient of the FGP beam and the structural dynamic responses. Once again, by implementing the proposed method, additional 12 cases of analysis were conducted and the results are presented in Fig. 6. As clearly demonstrated in Fig. 6, it can be observed that by increasing the porosity coefficients of the FGP beam, the normalized structural central deflection of the beam were also fluctuating with larger amplitude. That is, more severe vibration can be anticipated when the porosity coefficient of the FGP beam increases. Also, another interesting phenomenon can be observed is that if the velocity of the moving load was increased, but the excitation frequency was maintained constant, the absolute maximum amplitude of the normalized structural central deflection would decrease. Therefore, from both Figs. 5 and 6, it can be observed that the FGP beam would experience more violent vibration when the moving harmonic load is travelling at lower velocities.

in Fig. 3. In particular, the moving harmonic concentrated load considered here can be expressed as:

⎧Q (t ) = Q0sin(Ω·t ) L L L ⎨ x Q (t ) = v·t − 2 ,− 2 ⩽ x Q ⩽ 2 ,t1 = 0 ⩽ t ⩽ t2 = ⎩

L v

(28)

where Ω ∈ R + denotes the excitation frequency of the moving harmonic concentrated load. Once again, the normalized central deflection of the FGP beam is selected and the time-displacement responses of the FGP beam against various velocities and excitation frequencies are investigated by the proposed method. The time-displacement responses are presented in Fig. 5. From the second parametric study presented in Fig. 5, it can be observed that by increasing the excitation frequency of the moving harmonic load, the frequencies of the oscillations of the normalized structural central deflections were also increased. In order to thoroughly investigate the effects of the porosity coefficients (i.e., e0 ) of the FGP beam acting on the structural dynamic 295

Engineering Structures 165 (2018) 287–301

D. Wu et al.

materials which are governed by the symmetric material constitutive relationship (SMCR); all beams of the FGP frame structure are porous materials which are governed by the monotonic material constitutive relationship (MMCR); and all bracings are pure steel members. The properties of all incorporated materials are summarized in Table 7. The reason for such selection on the porosity types for different structural members is to demonstrate the applicability of the proposed FEM analysis framework on possessing the ability to simultaneously incorporate multiple porosity types only. There are not any practical reasons for such selection of material types. In the first part of this example, two free vibration analyses were conducted based on the Timoshenko and Euler-Bernoulli beam theories. By implementing the proposed method, the natural frequencies of the first 10 modes were calculated and their results are reported in Table 8 for Timoshenko beam theory (TBT) and Table 9 for Euler-Bernoulli beam theory (CBT). By examining all reported natural frequencies of the first 10 modes from Tables 8 and 9, it can be observed that all natural frequencies determined based on the Timoshenko beam theory are all smaller than the natural frequencies of the same structure but determined by the Euler-Bernoulli beam theory. Furthermore, it can be identified that for both beam theories, a general result convergence can be noticed at discretization of 100 FGP beam elements per structural member. Therefore, for all subsequent forced vibration analyses, a meshing size of 100 FGP beam elements per structural member is adopted. Consequently, there are in total 2100 FGP beam elements have been adopted for the forced vibration of the FGP frame structure within the proposed analysis framework. Also, another point should be highlighted here is that even for the finest mesh adopted in Tables 8 and 9 (i.e., 1000 FGP beam elements per structural member), only 137 s were consumed for an eigen-analysis with size of 21,000 × 21,000 on the PC with Intel® Core™ i7-4770 [email protected] GHz. In the second part of this example, the dynamic responses of the FGP frame structure subjected to various transient loads are investigated. For the purpose of this investigation, two types of loads are applied at the location as indicated in Fig. 7. The time history of these two loads is also presented in Fig. 8. In addition, the damping effect through the Rayleigh damping model is considered. For the purpose of demonstrating the applicability of the proposed approach to incorporate the damping effects, the two proportionality constants a0 = 2.62626 and a1 = 8.6328 × 10−7 were adopted from [36]. Consequently, by implementing the proposed approach, the structural time-displacement response of the top left node of the FGP frame structure can be calculated, and the results are presented in Figs. 9 and 10 for the analysis of the FGP frame structure subjected to loads 1 and 2, respectively. From Figs. 9 and 10, it can be observed that when the porosity coefficient was increased, the FGP frame structure would experience more critical vibration. Such phenomenon actually seems logical

Fig. 7. Three-storey single bay FGP frame structure.

Table 7 Structural information of three-storey single bay braced FGP frame structure.

E1 (GPa) ρ1 (kg/m3) ν b (m) h (m) Material model

FGP column

FGP beam

Steel bracing

200 7850 0.333 0.20 0.20 SMCR

200 7850 0.333 0.08 0.10 MMCR

200 7850 0.333 0.05 0.10 –

4.3. Three-storey single-bay braced FGP frame structure Unlike the first two examples, the third example presents an investigation on a three-storey single-bay braced FGP frame structure which involves multiple material types. Through the investigation of this complex FGP frame structure, the applicability and robustness of the proposed analysis framework can be further demonstrated. The schematic diagram of the general layout of the considered FGP frame structure is depicted in Fig. 7. Regarding the FGP frame structure considered in Fig. 7, there are three types of structural member that are presented within the system, which are including column, beam, and bracing. In total, there are 21 structural members involved in the concerned structural system. In addition, all three types of structural member are having different materials. In specific, all columns of the FGP frame structure are porous Table 8 Natural frequencies of the first 10 modes of the FGP frame structure based on TBT.

Numbers of FGP beam elements per structural member (TBT, e0 = 0.5 ) Natural frequency (rad/s)

5

10

50

100

500

1000

λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10

23.7008 80.5850 167.6070 214.7703 215.0515 217.3076 233.7083 236.4046 244.6142 334.5002

23.6923 80.5664 167.5690 214.6568 214.9367 217.1895 233.5658 236.2610 244.4604 334.3656

23.6894 80.5597 167.5542 214.6396 214.9192 217.1718 233.5437 236.2384 244.4359 334.3101

23.6893 80.5595 167.5536 214.6392 214.9188 217.1714 233.5432 236.2378 244.4353 334.3077

23.6893 80.5594 167.5534 214.6391 214.9187 217.1713 233.5430 236.2377 244.4351 334.3069

23.6893 80.5594 167.5533 214.6391 214.9187 217.1713 233.5430 236.2377 244.4351 334.3069

296

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Table 9 Natural frequencies of the first 10 modes of the FGP frame structure based on CBT. Numbers of FGP beam elements per structural member (CBT, e0 = 0.5) Natural frequency (rad/s)

5

10

50

100

500

1000

λ1 λ2 λ3 λ4 λ5 λ6 λ7 λ8 λ9 λ10

23.7610 80.7929 168.3945 215.6454 215.9283 218.1908 234.8008 237.5108 245.7550 337.0549

23.7610 80.7929 168.3945 215.6454 215.9283 218.1908 234.8008 237.5108 245.7550 337.0549

23.7498 80.7706 168.3626 215.5393 215.8206 218.0800 234.6684 237.3778 245.6114 336.9890

23.7498 80.7704 168.3625 215.5391 215.8204 218.0798 234.6682 237.3776 245.6112 336.9889

23.7496 80.7704 168.3624 215.5391 215.8204 218.0797 234.6681 237.3775 245.6111 336.9888

23.7609 80.7721 168.3628 215.5390 215.8203 218.0798 234.6682 237.3776 245.6111 336.9900

Fig. 8. (a) Time history of load 1 and (b) load 2.

Fig. 9. Horizontal displacement response of the top left node of the FGP frame structure subjected to load 1.

than the Euler-Bernoulli model. However, one exception can be observed at the minimum displacement of the FGP frame structure with e0 = 0.95 subjected to load 2 (i.e., Fig. 12(d)). In this particular case, the minimum structural displacement predicted by the Euler-Bernoulli beam is more critical than the one determined by the Timoshenko model. According to the currently presented numerical studies, there is not a dominant model between Timoshenko and Euler-Bernoulli beam theories which can provide absolute critical structural responses for the forced vibration analysis of FGP frame structures.

because that relatively weaker material can be anticipated for FGP beam due to higher porosity coefficient. Consequently, higher porosity coefficient of the FGP frame structure reduces the structural capacity against dynamic vibration. Furthermore, the time-domain responses of the FGP frame structure under different beam theories are examined. The corresponding results at the top left node are presented in Figs. 11 and 12 for loads 1 and 2, respectively. From both Figs. 11 and 12, it can be observed that all time-domain displacement responses of the FGP frame structure under the Timoshenko and Euler-Bernoulli beam theories generally agree very well with each other. However, by closely examining on Figs. 11 and 12, it can be clearly realized that the amplitudes of the displacement responses predicted by the Timoshenko beam model are generally larger

5. Conclusion In this paper, a finite element analysis framework is proposed for the free and forced vibration analyses of functionally graded porous 297

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Fig. 10. Horizontal displacement response of the top left node of the FGP frame structure subjected to load 2.

Fig. 11. Horizontal displacement responses of FGP frame structure subjected to load 1.

functionally graded porous beam problems as well as practically motivated complex functionally graded porous frame structures involving multiple material types. Therefore, the proposed method offers a more comprehensive and systematic analysis framework which can be implemented for real-life engineering applications. In order to demonstrate the applicability, accuracy and robustness of the proposed approach, three numerical examples have been thoroughly investigated in this study. By rigorously verifying results with both analytical and experimental results, the accuracy of the proposed approach has been thoroughly demonstrated. In addition, forced vibration analyses for both functionally graded porous beam and frame

frame structures. Within the proposed computational approach, both Young’s modulus and material density are graded through the thickness direction of the beam element. Two distinctive material models, namely the symmetric and monotonic material constitutive relationships (i.e., SMCR and MMCR), are incorporated. Thus, explicit formulations of the stiffness and mass matrices of the functionally graded porous beam element are presented for both the Timoshenko and Euler-Bernoulli beam models in this study. Moreover, the damping effects are also considered in the time-domain analysis through the implementation of the Rayleigh damping model. Unlike previous analytical approaches, the proposed computational scheme is competent for solving 298

Engineering Structures 165 (2018) 287–301

D. Wu et al.

Fig. 12. Horizontal displacement responses of FGP frame structure subjected to load 2.

Acknowledgements

structural systems subjected to diverse dynamic loading regimes have been investigated. Since this is a relatively less developed area of research comparing with finite element analysis of isotropic materials, further research works such as static and dynamic analyses of functionally graded porous frame structures with the considerations of material and geometrically nonlinearities will be developed. Also, finite element analysis (e.g., static, dynamics, linear and nonlinear analyses) for functionally graded porous plates and shells will be thoroughly proposed in the future.

The work presented in this paper has been supported by the National Natural Science Foundation of China (No. 51578166), Project of Chief Scientist of the Yangcheng Scholar (No. 1201541551), Technology Planning Project of Guangdong Province (No. 2016B050501004), Australian Research Council projects DP160103919 and IH150100006. The authors sincerely appreciate the reviewers for their valuable and constructive comments.

Appendix A. The components of the stiffness matrix through Timoshenko beam theory

K εi,11 = K εi,44 = −K εi,14 =

Ai bi Li

(A1)

K εi,13 = −K εi,16 = −K εi,34 = K εi,46 = −

K εi,22 = K εi,55 = −K εi,25 =

(A2)

12Dibi Li3 (1 + Φi)2

K εi,23 = K εi,26 = −K εi,35 = −K εi,56 = K εi,33 = K εi,66 =

Bibi Li

(A3)

6Dibi Li2 (1 + Φi)2

(A4)

Dibi (Φi2

+ 2Φi + 4) Li (1 + Φi)2

(A5)

299

Engineering Structures 165 (2018) 287–301

D. Wu et al.

K εi,36 = −

Dibi (Φi2 + 2Φi−2) Li

(A6)

I Φi = 24ksi (1 + νi ) ⎛⎜ 2i∼ ⎞⎟ ⎝ Li Ai ⎠ K γi ,22 = −K γi ,25 = K γi ,55 =

(A7)

AGi bi ksi

Φi2 (1 + Φi)2Li

K γi ,23 = K γi ,26 = −K γi ,35 = −K γi ,56 = AGi bi ksi

K γi ,33 = K γi ,36 = K γi ,66 =

(A8)

AGi bi Φi2 i ks 2(1 + Φi)2

(A9)

Li Φi2 4(1 + Φi)2

(A10) ∼ Ai and Ii denote the cross-sectional area and second moment of area of the ith FGP beam element, respectively; Li denotes the length, bi denotes the width of the ith element; and ksi = 1/ Ksi . Moreover, the elastic constants Ai , Bi , Di , and AGi of the ith FGP beam element can be explicitly determined as:

Ai =

0.5hi

∫−0.5h

E i (z ) dz =

i

Bi =

0.5hi

∫−0.5h

2hi i ⎧ ⎪ E1 (hi− π ) + ⎨ E i h − 2hi + ⎪ 1 i π ⎩

(

z·E i (z ) dz =

i

)

=

0.5hi

∫−0.5h

z 2·E i (z ) dz

i

AGi

=

0.5hi

∫−0.5h

Gi (z ) dz

(SMCR)

π 2E0i hi π

(MMCR)

(A11)

0(SMCR) ⎧ hi2 ⎨ 2 (4−π )(E1i−E0i )(MMCR) ⎩π i 3

Di

2E0i hi

(A12)

3

E1 hi h ⎧ − i 3 (E1i−E0i )(6π 2−48)(SMCR) ⎪ 12 12π = ⎨ E1i hi3 hi3 ⎪ 12 − 12π 3 (E1i−E0i )(6π 2 + 48π −192)(MMCR) ⎩

=

i

⎧ G1i (hi− 2hi ) + π

2G0i hi

⎨ i 2h G (h − i ) + ⎩ 1 i π

2G0i hi

π π

(A13)

(SMCR) (MMCR)

(A14)

Appendix B. The components of the mass matrix through Timoshenko beam theory

i i i M(1,1) = M(4,4) = 2M(1,4) =

bi Li i AM 3

i i i i M(1,2) = −M(1,5) = M(2,4) = −M(4,5) =

i i M(1,3) = M(4,6) =−

(B1)

bi i BM 2(1 + Φi)

(B2)

bi Li (1 + 4Φi) i BM 12(1 + Φi)

(B3)

i i M(1,6) = M(3,4) =

bi Li (1−2Φi) i BM 12(1 + Φi)

i i M(2,2) = M(5,5) =

i i 6DM AM bi Li bi (70Φi2 + 147Φi + 78) + 2 5Li (1 + Φi)2 210(1 + Φi)

i i M(2,3) = −M(5,6) =

i M(2,5) =

i AM bi Li2

840(1 + Φi)2

(B4)

(35Φi2 + 77Φi + 44) +

(B5)

i DM bi (1−5Φi) 10(1 + Φi)2

(B6)

i i 6DM AM bi Li bi (35Φi2 + 63Φi + 27)− 5Li (1 + Φi)2 210(1 + Φi)2

i i M(2,6) = −M(3,5) =−

i i M(3,3) = M(6,6) =

i M(3,6) =−

(B7)

i AM bi Li2 D i bi (1−5Φi) (35Φi2 + 63Φi + 26) + M 10(1 + Φi)2 840(1 + Φi)2

i AM bi Li3

840(1 + Φi)2

i AM bi Li3

840(1 + Φi)2

(7Φi2 + 14Φi + 8) +

(7Φi2 + 14Φi + 6) +

The mass constants

i AM ,

i BM

and

i DM

i DM bi Li

210(1 + Φi)2

i DM bi Li

30(1 + Φi)2

(B8)

(70Φi2 + 35Φi + 28)

(5Φi2−5Φi−1)

(B9)

(B10)

of the ith FGP beam element can be explicitly determined as: 300

Engineering Structures 165 (2018) 287–301

D. Wu et al.

i

i AM

=

0.5hi

∫−0.5h

i

i BM =

0.5hi

∫−0.5h

ρi (z ) dz

⎧ ρ i (hi− 2hi ) + 2ρ0 hi (SMCR) ⎪ 1 π π = 2ρ0i hi ⎨ i 2hi ⎪ ρ1 (hi− π ) + π (MMCR) ⎩

(B11)

0(SMCR) ⎧ hi2 ⎨ 2 (4−π )(ρ1i −ρ0i )(MMCR) ⎩π

(B12)

z·ρi (z ) dz =

i

i 3

i DM

=

0.5hi

∫−0.5h

i

z 2·ρi (z ) dz

3

ρ1 hi h ⎧ − i 3 (ρ1i −ρ0i )(6π 2−48)(SMCR) ⎪ 12 12π = ⎨ ρ1i hi3 hi3 i i ⎪ 12 − 12π 3 (ρ1 −ρ0 )(6π 2 + 48π −192)(MMCR) ⎩

(B13)

2016;107:39–48. [18] Kitipornchai S, Chen D, Yang J. Free vibration and elastic buckling of functionally graded porous beams reinforced by graphene platelets. Mater Des 2017;116:656–65. [19] Chen D, Yang J, Kitipornchai S. Nonlinear vibration and postbuckling of functionally graded graphene reinforced porous nanocomposite beams. Compos Sci Technol 2017;142:235–45. [20] Gibson LJ, Ashby M. The mechanics of three-dimensional cellular materials. Proc Roy Soc A: Math, Phys Eng Sci 1982;382(1782):43–59. [21] Choi J, Lakes R. Analysis of elastic modulus of conventional foams and of re-entrant foam materials with a negative Poisson’s ratio. Int J Mech Sci 1995;37(1):51–9. [22] Przemieniecki JS. Theory of matrix structural analysis. Dover Publication; 1985. [23] Reddy JN. Analysis of functionally graded plates. Int J Numer Meth Eng 2000;47(1–3):663–84. [24] Wu D, Gao W, Song C, Tangaramvong S. Probabilistic interval stability assessment for structures with mixed uncertainty. Struct Saf 2016;58:105–18. [25] Wu D, Gao W, Tangaramvong S. Time-dependent buckling analysis of concretefilled steel tubular arch with interval viscoelastic effects. ASCE J Struct Eng 2017;143(7):04017055-1–04017055-13. [26] Wu D, Gao W, Feng J, Luo K. Structural behaviour evolution of composite steelconcrete curved structure with uncertain creep and shrinkage effects. Compos Part B: Eng 2015;86:261–72. [27] Newmark NM. A method of computation for structural dynamics. ASCE J Eng Mech 1959;85(3):67–94. [28] Rao SS. Mechanical vibrations. 4th ed. Prentice-Hall; 2004. [29] Clough RW, Penzien J. Dynamics of structures. 2nd ed. McGraw-Hill INC.; 1993. [30] Paz M. Structural dynamics: theory and computation. 3rd ed. Chapman & Hall; 1991. [31] Su J, Banerjee JR. Development of dynamic stiffness method for free vibration of functionally graded Timoshenko beams. Comput Struct 2015;147:107–16. [32] Sina SA, Navazi HM, Haddadpour H. An analytical method for free vibration analysis of functionally graded beams. Mater Des 2009;30:741–7. [33] Wattanasakulpong N, Prusty BG, Kelly DW, Hoffman M. Free vibration analysis of layered functionally graded beams with experimental validation. Mater Des 2012:182–90. [34] Giunta G, Crisafulli D, Belouettar S, Carrera E. Hierarchical theories for the free vibration analysis of functionally graded beams. Compos Struct 2011;94:68–74. [35] Şimşek M, Kocatürk T. Free and forced vibration of a functionally graded beam subjected to a concentrated moving harmonic load. Compos Struct 2009;90(4):465–73. [36] Amini Y, Fatehi P, Heshmati M, Parandvar H. Time domain and frequency domain analysis of functionally graded piezoelectric harvesters subjected to random vibration: finite element modelling. Compos Struct 2016;136:384–93.

References [1] Liew KM, Zhao X, Ferreira AJM. A review of meshless methods for laminated and functionally graded plates and shells. Compos Struct 2011;93(8):2031–41. [2] Swaminathan K, Naveenkumar DT, Zenkour AM, Carrera E. Stress, vibration and buckling analyses of FGM plates—a state-of-the-art review. Compos Struct 2015;120:10–31. [3] Ferreira ADBL, Nóvoa PRO, Marques AT. Multifunctional material systems: a stateof-the-art review. Compos Struct 2016;151:3–35. [4] Shafiei N, Mirjavadi SS, Afshari BM, Rabby S, Kazemi M. Vibration of two-dimensional imperfect functionally graded (2D-FG) porous nano-/micro-beams. Comput Methods Appl Mech Eng 2017;322:615–32. [5] Shafiei N, Mirjavadi SS, Afshari BM, Rabby S, Hamouda AMS. Nonlinear thermal buckling of axially functionally graded micro and nanobeams. Compos Struct 2017;168:428–39. [6] Rezaei AS, Saidi AR. Application of Carrera Unified Formulation to study the effect of porosity on natural frequencies of thick porous–cellular plates. Compos Part B: Eng 2016;91:361–70. [7] Smith BH, Szyniszewski S, Hajjar JF, Schafer BW, Arwade SR. Steel foam for structures: a review of applications, manufacturing and material properties. J Constr Steel Res 2012;71:1–10. [8] Biot MA. Theory of buckling of a porous slab and its thermoelastic analogy. ASME J Appl Mech 1964;31(2):194–8. [9] Magnucki K, Stasiewicz P. Elastic buckling of a porous beam. J Theoret Appl Mech 2004;42(4):859–68. [10] Magnucka-Blandzi E. Dynamic stability if a metal foam circular plate. J Theoret Appl Mech 2009;47(2):421–33. [11] Jasion P, Magnucka-Blandzi E, Szyc W, Magnucki K. Global and local buckling of sandwich circular and beam-rectangular plates with metal foam core. Thin-Walled Struct 2012;61:154–61. [12] Belica T, Magnucki K. Stability of a porous-cellular cylindrical shell subjected to combined loads. J Theoret Appl Mech 2013;51(4):927–36. [13] Jabbari M, Farzaneh Joubaneh E, Khorshidvand AR, Eslami MR. Buckling analysis of porous circular plate with piezoelectric actuator layers under uniform radial compression. Int J Mech Sci 2013;70:50–6. [14] Jabbari M, Mojahedin A, Haghi M. Buckling analysis of thin circular FG plates made of saturated porous-soft ferromagnetic materials in transverse magnetic field. ThinWalled Struct 2014;85:50–6. [15] Chen D, Yang J, Kitipornchai S. Elastic buckling and static bending of shear deformable functionally graded porous beam. Compos Struct 2015;133:54–61. [16] Chen D, Yang J, Kitipornchai S. Free and forced vibrations of shear deformable functionally graded porous beams. Int J Mech Sci 2016;108–109:14–22. [17] Chen D, Kitipornchai S, Yang J. Nonlinear free vibration of shear deformable sandwich beam with a functionally graded porous core. Thin-Walled Struct

301