Laminated composite plates in contact with a bounded fluid: Free vibration analysis via unified formulation

Accepted Manuscript Laminated composite plates in contact with a bounded fluid: Free vibration analysis via unified formulation F.G. Canales, J.L. Mantari PII: DOI: Reference:

S0263-8223(16)31864-5 http://dx.doi.org/10.1016/j.compstruct.2016.11.079 COST 8040

To appear in:

Composite Structures

Received Date: Revised Date: Accepted Date:

15 September 2016 22 November 2016 25 November 2016

Please cite this article as: Canales, F.G., Mantari, J.L., Laminated composite plates in contact with a bounded fluid: Free vibration analysis via unified formulation, Composite Structures (2016), doi: http://dx.doi.org/10.1016/ j.compstruct.2016.11.079

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Laminated composite plates in contact with a bounded fluid: Free vibration analysis via unified formulation FG Canalesφ, JL Mantariφ1

φ

Faculty of Mechanical Engineering, Universidad de Ingeniería y Tecnología (UTEC)

Engineering and Technology, Jr. Medrano Silva 165, Barranco, Lima, Perú

Abstract The effects of the hydroelastic coupling are pronounced on composite structures due to the low composite density. This paper presents an analytical solution for the free vibration analysis of thick rectangular composite plates in contact with a bounded fluid. The plate displacement field is analyzed using Carrera Unified Formulation in order to consider displacement theories of arbitrary order. The displacement variables are approximated using the Ritz method. Analysis considering any of the classical boundary conditions is performed using suitable functions in the Ritz series. The fluid kinetic energy is obtained by considering an inviscid and incompressible fluid, and thus the velocity potential is applicable. The kinetic and potential energy of the plate are also obtained, and the energy functional is used to derive the eigenvalue problem. Validation is performed with results in the open literature and by using 3D finite element software. The influence of the fiber orientations on the natural frequencies is discussed, as well as the influence of the plate and fluid domain geometry. Certain natural frequencies are seen to be independent of the fluid domain height, depending on the symmetric or antisymmetric nature of the normal mode.

Keywords: Plate; Composite; Vibration; Fluid-structure interaction; Hydroelastic; Carrera Unified formulation (CUF).

1Corresponding Author email: [email protected] Tel: +00511 3540070; Cell:

+0051 96224551.

Nomenclature

a ,b

Plate length and width

c,d

Length and width of fluid domain

Cij

Constitutive matrix coefficients

D

Flexural rigidity of the plate

D p , Dnp , D nz

Linear differential operators

e

Depth of fluid domain

E

Young’s moduli

F

Fluid mass matrix

Fτ sij

Fluid mass nucleus

Fτ

Thickness expansion function

h

Plate thickness

K

Stiffness matrix

Kτ sij

Stiffness nucleus

j

Imaginary unit

J

Jacobian matrix

M

Ritz expansion order

M

Solid mass matrix

Mτ sij

Solid mass nucleus

N

CUF Expansion order

p,q

Indexes of trigonometric terms in x and y directions

P

Polynomial degree of Ritz expansion

t

Time

T

Kinetic energy of plate

TW

Kinetic energy of fluid

u, v , w

Plate displacements in x , y , z coordinates

u

Plate displacement vector

u

Plate amplitude displacement vector

U

Potential energy of plate

W

Amplitude of transverse plate deflection

x, y , z

Coordinates of plate

x
, y
, z


Coordinates of fluid domain

X , Y , Z , H

Assumed solutions of the velocity potential in x
, y
, z


axes and in time

β = e /d

Depth to width ratio of fluid domain

εn , ε p

Vector of in-plane and out-of-plane strain components

φ

Velocity potential

Γ P , ΓW

Plate and fluid area in the bottom

λ = c /d

Length to width ratio of fluid domain

ν

Poisson’s ratio

ρ , ρW

Density of structure and fluid

σn , σ p

Vector of in-plane and out-of-plane stress components

ω

Frequency of vibration

ω

Non-dimensional frequency of vibration

Ω

Fluid domain

ξ ,η

Non-dimensional x and y coordinates of plate

ξ
,η
, ζ


Non-dimensional x
, y
, z
coordinates of fluid domain

ψ u ,ψ v ,ψ w

Shape functions of the plate displacements u , v , w

Ψ

Shape functions matrix

∇

Del operator

1. Introduction Composites are now widely used in shipbuilding, aerospace and civil structures due to their superior strength to weight ratio compared to traditional materials. Many applications require modeling the dynamic response of the structures in contact with water in order to obtain the natural frequencies and avoid resonance. The effects of the hydroelastic coupling, as measured by comparing the natural frequencies of the plate in contact with water and in vacuum, are more pronounced on composite materials compared to other metallic materials such as steel since composites have lower density. Accurate modeling of the interaction between the fluid and the structure requires a complex analysis due to the hydroelastic coupling. Finite element and boundary element models can obtain a solution to the problem but at a high computational cost. Analytical methods with appropriate simplifications can obtain accurate results with a very low computational cost, and are useful to obtain a qualitative understanding of the interaction. Analytical solutions for the vibration of a plate in contact with a fluid domain have been obtained for a variety of fluid and plate arrangements. Analysis of circular plates assuming that the dry modal functions can approximate the wet dynamic displacements has been developed for incompressible [1] and compressible [2] fluid domains. Vibration of circular plates in asymmetric conditions has been investigated by Tariverdilo et al. [3]. Viscosity and damping has been introduced in the work by Phan et al. [4] and Kozlovsky [5]. Finite element modeling using 2D plate elements has been developed by Kerboua et al. [6] and Bermudez et al. [7]. Closed-form solutions using Mindlin plate theory and certain boundary conditions have been obtained by Hashemi et al. [8]. Other developments have analyzed the modal energy associated with the fluid and the plate [9], and the behavior of magneto-electro-elastic plates [10]. Experimental

investigation for the vibrational behavior of structures in contact with fluid is given in Refs. [11] and [12]. The assumption of equal dry and wet modal displacements is seen to be less accurate for composite structures. The analysis of rectangular plates coupled with an inviscid, finite fluid domain is an important topic for structural design of fuel tanks and liquid containers. A formulation using added mass terms has been introduced by Kramer et al. [13]. A coupled hydroelastic analysis of rectangular plates in contact with a fluid domain has been presented by Cheng and Zhou [14] using the Kirchhoff plate theory and Ritz method. Further development has been done by considering Mindlin plate theory and stiffeners [15], compressibility of the fluid [16], elastic foundations [17], in-plane loads [18], geometric non-linearity and sloshing effects [19] and harmonic excitation forces [20]. Analysis of multiple plates in contact with fluid, applicable for fuel assemblies in a research reactor, has been developed by Jeong and Kang [21]. The majority of the references above-mentioned use Kirchhoff or Mindlin plate theory, since most of them consider thin isotropic plates. However, accurate shear strain modeling is required for thick laminated plates. Higher-order shear deformation theories (HSDTs) can consider nonuniform shear distributions in the plate cross-section and are used to obtain more accurate results. Analytical derivations using HSDTs of arbitrary order can be performed easily using the Carrera’s Unified Formulation (CUF). This formulation can obtain accurate results compared to 3D analysis and is computationally efficient. The formulation was presented in Ref. [22]. Analysis of thermal stresses in plates [23-25], multifield problems [26-29], shells [30-31] and functionally graded materials [32] has been developed using this formulation. Analysis of plates using nonpolynomial displacement fields has been developed in Refs. [33-36]. The core of the formulation is described in detail in Refs. [37-39]. In the hydroelastic problem, the displacement field of the plate must be approximated in some manner. The Ritz method is a common choice due to the capability to use displacement shape functions that can consider arbitrary boundary conditions. The method is described in detail in Refs. [40-42]. The choice of the shape functions has a great influence on the accuracy and convergence of the results. Analysis of simply supported plates using trigonometric shape functions has been presented by Fazzolari and Carrera [43-46]. Analysis of shell elements using Ritz method has been developed in Refs. [47-49]. Free vibration analysis of plates with arbitrary boundary conditions has

been presented by Dozio [50-53], Vescovini and Dozio [54], and in the paper by Dozio and Carrera [55]. The present paper develops an analytical solution for the free vibration of thick rectangular composite plates coupled to an inviscid and incompressible fluid domain. The main novelty relies on the analysis of the hydroelastic problem of composite plates by coupling a hydrodynamic model with a structural model that can consider displacement theories of arbitrary order, using Carrera Unified Formulation. The displacement field is approximated using the Ritz method. Validation of the results is performed by using a 3D finite element solution and open literature for comparison. Parametric studies are performed to assess the influence of various geometrical parameters and the fiber orientation on the natural frequencies of the plate.

2. Analytical modeling A rectangular plate with length a , width b and thickness h is considered, as shown in Figure 1. The motion of the plates is described by using the coordinate system x − y − z . The plate is in contact with an incompressible, inviscid and irrotational fluid bounded by a rectangular domain with length c , width d and depth e , as shown in Figure 2. The fluid motion is described by the coordinate system x
− y
− z
. Only the bulging modes are analyzed and sloshing effects are neglected. Consequently, a simplified free surface condition is sufficient. The vertical walls and the bottom of the fluid domain are assumed to be rigid, except for the region in contact with the top surface of the plate at

z
= 0 or z = h /2 . The plate is assumed to be linearly elastic and the amplitude of motion is assumed to be small. The energy functional is used to obtain the eigenvalue problem. In section 2.1, a structural model is used to develop the plate potential energy

U and plate kinetic energy T with resemblance to the formulation given by Dozio and Carrera [55], but adapted to composites. In section 2.2, a hydrodynamic model develops the kinetic energy of the fluid TW in a similar manner as presented by Cheung and Zhou [14]. In section 2.3, the use of the Ritz method and the choice of the shape functions is described. In section 2.4, the stiffness and mass nuclei, and the assembly of the matrix is detailed. Finally, in section 2.5 the fluid and mass matrix are used to deduce the eigenvalue problem.

2.1 Structural model For convenience, the derivation of the structural model is given in non-dimensional form. The x − y domain (0 ≤ x ≤ a , 0 ≤ y ≤ b ) of the rectangular plate is mapped into a computational ξ − η domain (-1 ≤ ξ ,η ≤ 1) by using a transformation of coordinates given by:

ξ=

2x −1 a

η=

,

2y −1 b

(1)

The derivatives in the two coordinate systems are related by:

∂ 2 ∂ = ∂ x a ∂ξ

∂ 2 ∂ = ∂ y b ∂η

,

(2)

A general displacement vector is introduced: T

u (ξ ,η , z, t ) = {u (ξ ,η , z, t ) v (ξ ,η , z , t ) w (ξ ,η , z, t )}

(3)

where u , v and w are the displacement components in ξ , η and z axes respectively, and t is time. The stress and strain components are grouped as follows:

{ } T , ε p = {ε xx ε yy ε xy } T σ n = {σ xz σ yz σ zz } T , ε n = {ε xz ε yz ε zz } T σ p = σ xx σ yy σ xy

(4)

Considering small displacements, the strain-displacement relations are:

ε p = D pu

(

(5)

)

ε n = Dnp + Dnz u where the linear differential operators D p , Dnp and Dnz are given by:

∂  ∂x   Dp =  0  ∂   ∂y

0 ∂ ∂y ∂ ∂x

  0 0     0  , D np =  0     0 0  

0

0 0

∂ ∂ ∂z ∂x   ∂ , Dnz =  0  ∂y     0 0  

0 ∂ ∂z 0

 0   0    ∂  ∂ z 

(6)

and the derivatives must be evaluated using Eq. (2). The stress components are given by the constitutive law:

σ = Cε

(7)

Eq. (7) can be split by using Eq. (4):

σ p = C pp ε p + C pn ε n

(8)

σ n = Cnp ε p + Cnn ε n where the matrices C pp , C pn , C np and Cnn are given by:

C pp

 C11 = C12 C16

T C pn = C np

C12 C22 C26  0 =  0  0

C16   C55  C26  , C nn =  C45  0 C66  0 0 0

C45 C44 0

0  0  C33 

C13  C23  C36 

(9)

The stiffness coefficients Cij depend on the material properties and the fiber orientation. For the sake of brevity, the dependence is not described here. It can be found on standard texts of composite materials, as in Daniel and Ishai [56] or Reddy [57]. The displacement field is expressed within the framework of CUF:

u (ξ ,η, z, t ) = Fτ ( z ) uτ (ξ ,η, t )

τ = 0,1,..., N

(10)

where Fτ are functions of the coordinate z on the cross-section, N stands for the number of terms used in the expansion, uτ is the vector of the generalized displacements, and the repeated subscript “ τ ” indicates summation. A simple polynomial expansion is used to determine the functions Fτ . For example, the displacement field of the second-order ( N = 2 ) expansion model can be expressed as:

u = u0 + zu1 + z 2u2 v = v0 + zv1 + z 2v2 2

w = w0 + zw1 + z w2

(11)

Reduced order plate theories can be obtained by using a suitable expansion order and eliminating certain displacement variables. For example, the Mindlin plate theory is obtained by imposing kinematic constraints in a first-order expansion ( N = 1 ), resulting in the following displacement field:

u = u0 + zu1 v = v0 + zv1

(12)

w = w0 The potential energy U of the plate in the computational domain can be written as follows:

U=

1 1 1 h/2 T ε pσ p + εTnσ n J dz dξ dη ∫ ∫ ∫ − 1 − 1 − h /2 2

(

)

(13)

where J stands for the jacobian of the coordinate transformation given by Eq. (1). Substituting Eqs. (8) and (5) in Eq. (13) the strain energy is given by:

U=

1 1 1 h /2  D pu 2 ∫−1 ∫−1 ∫− h /2 

(

(

T

) ( C pp D pu + C pn Dnpu + C pnDnzu ) + T

) ( Cnp D pu + Cnn Dnpu + Cnn Dnz u ) T + ( Dnz u ) ( Cnp D p u + Cnn Dnp u + Cnn Dnz u )  J dz dξ dη  + Dnp u

(14)

The kinetic energy T of the plate is given by the following expression: T

1 1 1 h /2  ∂u   ∂u  ρ     J dz dξ dη T= ∫ ∫ ∫ − 1 − 1 − h / 2 2  ∂t   ∂t  where

(15)

ρ is the mass density of the plate.

2.2 Hydrodynamic model The fluid is considered to be incompressible, inviscid and irrotational. A simplified surface condition is considered appropriate since sloshing effects are neglected. The velocity potential

φ ( x
, y
, z
) is used to describe the motion in the fluid domain Ω . The

Laplace equation associated to the velocity potential is given by:

∇ 2φ =

∂ 2φ ∂ 2φ ∂ 2φ + + =0 ∂ x
2 ∂ y
2 ∂ z
2

in Ω

(16)

The boundary conditions in the fluid domain are given by a zero velocity normal to the rigid walls and the rigid bottom, except for the plate surface, and no surface waves on the top of the fluid domain. It is assumed that the fluid particles are perfectly attached to the plate surface, with no diffusion or cavitation:

∂φ = 0, ∂x


on x
= 0, c

( rigid wall )

(17)

∂φ = 0, ∂y


on y
= 0, d

( rigid wall)

(18)

φ = 0,

on z
= e

( no surface waves )

(19)

−

∂φ ∂z


z
=0

∂  w / ∂t = 0

( ΓP ) on the other part ( ΓW ) on plate surface

(20)

where the plate displacement w couples the hydrodynamic model to the structural model. The amplitude W of the transverse plate displacement is defined as follows:

w(ξ ,η, z, t ) = W (ξ ,η, z ) T ( t )

(21)

The solution of the Laplace equation (16) is assumed in the following form:

φ ( x
, y
, z
, t ) = X ( x
) Y ( y
) Z ( z
) H ( t )

(22)

Substituting Eq. (22) in Eq. (16), the following three ordinary differential equations are obtained:

d2 X ± px2 X = 0 2 dx


(23)

d 2Y ± p 2yY = 0 dy
2

(24)

d2 Z ∓ px2 + p y2 Z = 0 dz
2

(

)

(25)

Solving Eqs. (23) – (25) and considering the boundary conditions given in Eqs. (17) – (19), the solution of the velocity potential of the fluid is given by:

∞

φ ( x
, y
, z
, t ) = eH ( t ) ∑ p=0

∞

Apq cos ( pπξ
) cos ( qπη
) Fpq (ζ
) ∑ q=0

(26)

where

Fpq

1 − ζ

ζ = f pq ( 2 −ζ
) f ζ
e pq − e

if p = q = 0

( )

f pq = πβ

( p /λ )

2

(27) elsewhere

+ q2

(28)

λ = c /d , β = e /d and non-dimensional coordinates have been introduced for the fluid domain:

ξ
= x
/ c , η
= y
/ d , ζ
= z
/ e

(29)

Substituting Eq. (26) and Eq. (21) in Eq. (20), using the orthogonality of the trigonometric series, and since the fluid domain bottom surface is located at z = h /2 in the plate coordinates, the coefficient Apq is derived as follows:

Apq =

ε pq

Q pq ∫∫Γ

P

(

)

W (ξ ,η , z ) z =h /2 cos pπξ
cos ( qπη
) dξ
dη


(30)

where

ε pq

1  = 2 4 

if p = q = 0 if p = 0, q ≠ 0 or p ≠ 0, q = 0

(31)

elsewhere

1 Q pq =  2 f pq  − f pq 1 + e

(

if p = q = 0

)

elsewhere

(32)

The kinetic energy of the fluid is given by:

1 2 TW = ρW ∫∫∫ ( ∇φ ) dΩ Ω 2

(33)

where ρW is the fluid density and ∇ is the del operator. By applying Green’s first identity to Eq. (33) in conjunction with Eq. (16), a surface integral involving the fluid

domain boundary is obtained. By using the boundary conditions given by Eqs. (17) – (20), the integral is simplified to a single surface integral in the plate region Γ P :

1  ∂φ  TW = − ρW ∫∫  φ  ΓP 2  ∂z


(34)

dx
dy
z
=0

Substituting Eqs. (20), (21), (26), (29) and (30) in Eq. (34), the fluid kinetic energy is given by: 2 1 ρW cde  H ( t )  2 ∞ ∞ ε  pq × ∑ ∑ Fpq p =0 q =0  Qpq

TW =

(

2 



∫∫ΓP W (ξ ,η , z ) z =h/2 cos pπξ cos ( qπη ) dξ dη  

(

)

)

(35)

where

1 Fpq ζ
=  2 f pq 1 − e

( )

if p = q = 0

(36)

elsewhere

2.3 Ritz series A harmonic motion is assumed for the generalized displacements given in Eq. (10):

uτ (ξ ,η , t ) = uτ ( ξ ,η ) e jωt

τ = 0,1,..., N

u ( ξ ,η , z, t ) = Fτ ( z ) uτ ( ξ ,η ) e jωt

(37) (38)

where uτ are the amplitudes of the generalized displacements and j is the imaginary unit. The amplitudes are approximated as follows: M

uτ (ξ ,η ) = ∑ cuτ i ψ uτ i (ξ ,η )

τ = 0,1,..., N

i =1 M

vτ ( ξ ,η ) = ∑ cvτ i ψ vτ i (ξ ,η )

τ = 0,1,..., N

i =1

M

wτ (ξ ,η ) = ∑ cwτ i ψ wτ i (ξ ,η ) i=1

τ = 0,1,..., N

(39)

where M is the order of the Ritz expansion, cuτ i , cvτ i , cwτ i are unknown coefficients and ψ uτ i , ψ vτ i , ψ wτ i are assumed shape functions. These shape functions are assumed as follows:

ψ θτ i (ξ ,η ) = pi (ξ ,η )ψ bθτ (ξ ,η )

θ = u , v, w

(40)

b where pi is a generic element of a complete 2D polynomial space and ψ θτ (ξ ) is a

function used to satisfy the boundary conditions. A generic pi term is given by:

pi = ξ q −rη r

q = 0,1,..., P

and

r = 0,1,..., q

(41)

where P is the degree of the polynomial space. The indices are related as follows:

i=

( q + 1)( q + 2 ) −

(q − r )

2

(42)

The order M of the Ritz expansion is related to the polynomial degree P as follows:

M=

( P + 1)( P + 2 )

(43)

2 b

The boundary-compliant function ψ θτ is given by the following expression: α1θτ

ψ bθτ (ξ ,η ) = (1 + ξ )

α 2θτ

(1 + η )

α3θτ

(1 − ξ )

α 4θτ

(1 − η )

(44)

The expressions in parentheses are the equations of the plate boundaries. The exponents

α iθτ are chosen according to the boundary conditions in the edges. If the i-edge is fully clamped then the geometric boundary conditions are u = v = w = 0 which implies that

α iθτ = 1

(θ = u , v , w) . If the i-edge is simply supported then transverse and tangential

displacements are imposed to be zero. For example, if it is a simply supported edge parallel to the x axis then the geometric boundary conditions are u = w = 0 which implies that αiuτ = αiwτ = 1 and since displacement in the y axis, i e. the displacement variable v , has no geometric restriction, then αivτ = 0 . In matrix notation, Eq. (39) can be written as follows:

uτ (ξ ,η ) = Ψτ i cτ i

τ = 0,1,..., N

(45)

where:

ψ uτ i (ξ ,η )  0 0   ψ vτ i (ξ ,η ) Ψτ i (ξ ,η ) =  0 0   0 0 ψ ξ , η )  wτ i (  cτ i = cuτ i

cvτ i

cwτ i 

T

(46)

(47)

2.4 Stiffness and mass nuclei and matrix assembly For simplicity of notation, the following generic terms are introduced: c+ d ∂ a +b ψ θτ i ∂ ψ δ sj =∫ ∫ J dξ dη −1 −1 ∂x a ∂y b ∂xc ∂y d 1

1

Eταβ =∫ ,γ s,θ

h /2

abcd Iθδ

− h /2

Eτ s = ∫

Cαβ Fτ,γ Fs,θ dz

h /2

F F dz − h /2 τ s

(48)

(49)

(50)

Substituting Eq. (38) and (45) in Eq. (14), the maximum plate potential energy is given by:

1 U MAX = cτTi Kτ sij csj 2 where K τsij is the stiffness nucleus, and is given by: 1010 0000 0101 1001 0110 Kτ sij (1,1) = Eτ11s Iuu + Eτ55,z s ,z Iuu + Eτ66s Iuu + Eτ16s Iuu + Eτ16s Iuu 1001 0110 1010 0101 0000 Kτ sij (1, 2 ) = Eτ12s I uv + Eτ66s I uv + Eτ16s I uv + Eτ26s I uv + Eτ45,z s,z I uv 1000 0010 0100 0001 Kτ sij (1,3) = Eτ13s,z I uw + Eτ55,z s I uw + Eτ36s,z I uw + Eτ45,z s I uw 0110 1001 1010 0101 0000 Kτ sij ( 2,1) = Eτ12s I vu + Eτ66s I vu + Eτ16s I vu + Eτ26s I vu + Eτ45,z s,z I vu 1010 0110 1001 Kτ sij ( 2, 2 ) = Eτ22s I vv0101 + Eτ44,z s,z I vv0000 + Eτ66s I vv + Eτ26s I vv + Eτ26s I vv 0100 0001 1000 0010 Kτ sij ( 2,3) = Eτ23s,z I vw + Eτ44,z s I vw + Eτ36s,z I vw + Eτ45,z s I vw

(51)

0010 0001 0100 36 45 Kτ sij ( 3,1) = Eτ13,z s I wu + Eτ55s,z I 1000 wu + Eτ , z s I wu + Eτ s , z I wu 0001 0100 0010 1000 Kτ sij ( 3,2 ) = Eτ23,z s I wv + Eτ44s,z I wv + Eτ36,z s I wv + Eτ45s ,z I wv 0000 0101 1010 1001 0110 Kτ sij ( 3,3) = Eτ33,z s,z I ww + Eτ44s I ww + Eτ55s I ww + Eτ45s I ww + Eτ45s I ww

(52)

The stiffness matrix is constructed by expanding the nucleus over the indices τ , s , i and j . By varying the indices τ and s of the stiffness nucleus over the range

τ , s = 0,1,..., N the following matrix is obtained:  K 00ij  K ij =  ...   K N 0ij

... Kτ sij ...

K 0 Nij   ...   K NNij 

(53)

The stiffness matrix is obtained by varying the indices i and j over the range

i, j = 1, 2,..., M :

 K11 ...  K =  ... K ij K  M 1 ...

K1M   ...  K MM 

(54)

Substituting Eq. (38) and (45) in Eq. (15), the maximum plate kinetic energy is:

1 TMAX = ω 2cτTi Mτ sij csj 2

(55)

where Mτ sij is the solid mass nucleus, and is given by:

Mτ sij = ∫

1

1

Ψ ρ Eτ s Ψ sj J dξ dη −1 ∫−1 τ i T

(56)

The components of the solid mass nucleus in Eq. (56) are given by: 0000 Mτ sij (1,1) = Eτρs I uu

Mτ sij ( 2, 2 ) = Eτρs I vv0000 0000 Mτ sij ( 3,3) = Eτρs I ww

(57)

where

Eτρs = ∫

h /2

− h /2

ρ Fτ Fs dz

(58)

The time independent part of the transverse displacement w can be expressed in a convenient form using Eqs. (21), (37) and (38) as follows:

W (ξ ,η , z ) = Fτ ( z ) wτ (ξ ,η )

(59)

Substituting Eq. (59) and Eq. (45) into Eq. (35), the maximum fluid kinetic energy is:

1 TW.MAX = ω 2cτTi Fτ sij csj 2

(60)

where Fτ sij is the fluid mass nucleus, and is given by: ∞

Fτ sij = ρW cde ∑

∞

 ε pq

∑ p=0 q=0 Q 

 Fpq ( Fτ Fs ) z =h /2 Θτ ipq ΘTsjpq  

pq

(

)

Θτ ipq = ∫∫ ψ wτ i cos pπξ
cos ( qπη
) dξ
dη
ΓP

ψ wτ i = 0

ψ wτ i 

0

T

(61)

(62)

(63)

The only non-zero component of the fluid mass nucleus in Eq. (61) is given by:

Fτ sij ( 3,3) =

∞

∞

Eτρs ∑ ∑ p =0 q=0

W

ε pq Qpq

Fpq Θτ ipq Θsjpq

(64)

where

EτρsW = ρW cde ( Fτ Fs ) z =h /2

(

(65)

)

Θτ ipq = ∫∫ ψ wτ i cos pπξ
cos ( qπη
) dξ
dη
ΓP

(66)

The solid mass matrix M and fluid mass matrix F can be obtained by expanding the indices τ , s , i and j of the mass nucleus, in a similar way as the stiffness matrix is constructed. 2.5 Eigenvalue problem and solution

The energy functional is given by:

Π = UMAX − (TMAX + TW.MAX )

(67)

Substituting Eqs. (51), (55) and (60) in Eq. (67), and minimizing the functional with respect to the undetermined coefficients cτ i , the following equation is obtained:  K − ω 2 ( M + F )  c = 0

(68)

where c is a vector of unknown coefficients. By solving Eq. (68) the natural frequencies and mode shapes of vibration can be obtained. In order to obtain the j vibrational mode it is required to perform a summation of the set of admissible functions evaluated at the point of interest times the term i of the eigenvector X j for each of the generalized displacement variables. Afterwards, the use of Eq. (10) is required in order to obtain the displacement amplitudes. For the displacement in the z axis, it is given as follows: N

Mode j ( ξ ,η , z ) = ∑

M

∑ X j,wτ i Fτ ( z )ψ wτ i (ξ ,η )

(69)

τ =0 i = 1

3. Numerical results and discussion The present formulation has been coded into a MATLAB program, and numerical results are given in the present section. A four-letter symbolic notation is used describe the boundary conditions of the plate. Letters S and C stand for simply supported and clamped edges. The boundaries are numbered in a counterclockwise manner starting from the edge at x = 0 . Graphite-epoxy and boron-epoxy laminated materials with properties given in Table 1 are considered for further analysis. The materials are denoted by MAT1 and MAT2. The fluid density is assumed to be ρW = 1000 kg/m 3 . The plate geometric center is assumed to be coincident with the geometric center of the bottom of the fluid domain in all the numerical results. Natural frequencies are given in Hz. Thickness locking is corrected for first order models ( N = 1 ) as in Carrera and Brischetto [58]. In the numerical results, PM stands for present method. The natural frequencies are given in Hz.

3.1 Convergence study In this section, the convergence of the results is analyzed in order to verify the result stability as the number of trigonometric terms p , q and the polynomial degree of the Ritz expansion P is increased. A CCCC square (0/90°) MAT1 laminated plate with geometry parameters given by a = b = 1 m and h / b = 0.1 is considered. The fluid domain geometry is assumed as c = d = e = 1 m . Table 2 presents the natural frequencies of the plate in contact with fluid as the polynomial degree P is increased, where the number of trigonometric terms used is p, q = 16 . Slower convergence is observed for a higher expansion order N . Results within 0.15% of the converged solution can be obtained by using a polynomial degree P = 8 . For the purposes of the present work, the polynomial degree P = 11 is considered sufficiently accurate and it is further used in the manuscript. Table 3 presents the natural frequencies of a CCCC square (0/90°) MAT1 laminated plate in contact with fluid as the number of trigonometric terms p, q is increased. Results are given for N = 5 and considering equal trigonometric terms in the x
and y
directions. The plate geometry is the same as in the previous case, and the fluid domain has a square bottom, i.e. c = d . Convergence is slower for larger fluid domains and for higher modes of vibration. Further numerical results are obtained considering d / b < 3 , and thus the number of trigonometric terms used in the remainder of the manuscript is

p, q = 18 in order to obtain accurate results. 3.2 Validation of the results

In this section, validation of the results is performed by developing a numerical example and comparing the results with a 3D finite element solution and open literature. ANSYS general purpose program is used to obtain the natural frequencies of a plate coupled with a fluid domain. The plate is considered to be a CCCC square (0/60°) plate with geometry given by a = b = 2 m . Two plate thickness ratios will be considered:

h / b = 0.1 (moderately thick plate) and h / b = 0.25 (thick plate). The fluid domain width and length are considered the same as the plate, and the height is assumed as

e = 1 m . The three-dimensional model is composed of three-dimensional fluid elements (FLUID30) and solid elements (SOLID185) with identical size. The number of elements is doubled until the difference of the natural frequencies with the previous

result is less than 0.15 % . For the case with h / b = 0.1 , the plate is segmented into 172800 (120 x 120 x 12) solid elements and the fluid domain is divided into 864000 (120 x 120 x 60) fluid elements, as shown in Figure 3. For the case with h / b = 0.25 , 170368 (88 x 88 x 22) solid elements and 340736 (88 x 88 x 44) fluid elements are used. Viscosity and compressibility of the fluid are neglected. The boundary conditions for the fluid domain are the same as given in Eqs. (18), (19) and (20). The fluidstructure interaction flag is activated in fluid elements adjacent to solid elements. Table 4 presents the first five natural frequencies of a plate with h / b = 0.1 in contact with a fluid domain considering MAT1 and MAT2 materials. Results for N = 1 without the correction of thickness locking given in Ref. [58] are denoted as N = 1 TL . It can be observed that arbitrarily low errors can be obtained by increasing the expansion order

N . While highly accurate results are obtained with N = 5 , from an engineering point of view using

N = 3 yields acceptable errors. Lower differences are observed for

MAT2, due to the lower orthotropicity ratio E1 / E2 compared to MAT1. In addition, the effect of the thickness locking correction is dependent on the material used. When MAT2 is used and the correction is applied, the error obtained using N = 1 is slightly lower than using N = 2 . Then the thickness locking correction is further utilized in the manuscript. However, it must be mentioned that this correction does not have a consistent theoretical proof. Figure 4 shows the first five mode shapes using MAT1 obtained by the present method and by ANSYS software, where good agreement is observed. The present formulation can also analyze thicker plates with similar precision. Table 5 presents the natural frequencies of a plate coupled with fluid considering

h / b = 0.25 for MAT1 and MAT2 materials. The capability of the present model to accurately predict natural frequencies of thick plates in contact with fluid is evident. As far as the authors are aware, there are no results in the open literature for laminated plates coupled to a fluid domain with same geometry as used in the present manuscript. In Refs. [17, 18] some studies on vertical laminated plates can be found. Comparison of results from the open literature and the ones presented in this paper is presented in Table 6, where an isotropic plate with h / b = 0.05 and ν = 0.3 coupled to a cubic fluid domain is analyzed. CPT stands for classical plate theory and MPT stands for Mindlin plate theory. Non-dimensional natural frequencies are given, defined as follows:

ω = ωa

2

12 ρ (1 − v 2 ) Eh 2

(70)

Good agreement between the results is observed. It must be mentioned that results from Ref. [17] consider a shear correction factor k = 5 / 6 . On the other hand, the results from the present method given for N = 1 do not use such shear correction factor. 3.3 Parametric studies

In this section, natural frequencies of a plate with stacking sequence (0/θ°) are obtained for various fiber orientation angles θ. The influence of various parameters is assessed by establishing a baseline case with a = b = c = d = e = 1 m , h / b = 0.1 and varying the parameters one at a time. The results are given for an expansion order N = 5 . 3.3.1 Effect of thickness ratio

Table 7 presents the first five natural frequencies of a CCCC MAT1 plate in contact with fluid for various fiber orientation angles and plate thickness ratios. As in the dry vibration case, higher natural frequencies are obtained as the thickness is increased. However, the relative influence of the fluid on the natural frequencies is dependent on the thickness ratio. Figure 5 shows the fundamental frequency of a CCCC MAT1 plate for ply orientation angles (0°), (0/90°) and various thickness ratios. It can be observed that the added mass effect of the fluid can dramatically decrease the natural frequencies of thin plates. 3.3.2 Effect of aspect ratio

Table 8 presents the first five natural frequencies of a CCCC MAT1 plate in contact with fluid for various plate aspect ratios with b = 1 m as a fixed dimension. As the aspect ratio increases, i.e. the plate length is higher, the natural frequencies becomes more dependent of the rigidity in the Y axis. For this reason, laminas oriented at 90° increase the natural frequencies at high aspect ratios. Figure 6 shows the fundamental frequency of a CCCC MAT1 plate plotted against the plate aspect ratio for various stacking sequences. At large aspect ratios the fundamental frequency approaches a constant value, which is higher as the ply orientation angle gets closer to 90°. 3.3.3 Effect of fluid domain depth

Table 9 presents the first five natural frequencies of a CCCC MAT1 plate coupled to a fluid domain for various fluid domain depths considering b = 1 m as a fixed dimension. It can be observed that certain frequencies are independent on the fluid domain depth, depending on the fiber orientation. The natural frequencies ω 2 and ω 3 are independent on the fluid depth for all the fiber orientations analyzed. The frequency ω4 for a lamination (0/90°) and the frequency ω5 for a lamination (0/0°) are also independent on the fluid depth. In order to further study this point, the mode shapes of a CCCC MAT1 plate considering e / b = 3 and fiber orientations (0/0°) and (0/90°) are shown in Figure 7. The shape is unchanged as the fluid domain depth increases. The following comments arise: •

The mode (0,0), associated to the frequency ω 1 , and the mode associated to the frequency ω 5 for a lamination (0/90°) are doubly symmetric and the natural frequencies are dependent on fluid domain depth.

•

The modes (1,0) and (0,1) associated to the frequencies ω 2 and ω 3 are symmetric – antisymmetric and the natural frequencies are independent on fluid domain depth. The mode (0,2) associated to the frequency ω4 for a lamination (0/0°) is only slightly dependent on fluid depth.

•

The mode (1,1), associated to the frequencies ω 5 and ω4 for laminations (0/0°) and (0/90°) respectively, is doubly antisymmetric and the natural frequencies are independent on fluid domain depth.

This behavior can be explained by noting that doubly symmetric modes produce a maximum of vertical fluid displacement, and thus the natural frequencies should be highly dependent on the fluid domain depth. On the other hand, antisymmetric modes produce a zero net fluid displacement and thus the effect of fluid depth on the associated natural frequencies should be much lower. 3.3.4 Effect of fluid domain width and boundary conditions

Table 10 presents the first five natural frequencies of a CCCC MAT1 plate for various fluid domain widths, considering d = c and b = 1 m . As the fluid domain size increases, the natural frequencies increase. The added mass effect of the fluid is slightly reduced since the fluid lateral motion is less restricted. The natural frequencies approach

a constant value, corresponding to an infinite fluid domain. Table 11 presents the first five natural frequencies of a MAT1 plate in contact with fluid for various plate boundary conditions. More rigid supports increase the natural frequencies. However, the relative orientation of the clamped supports is seen to play a role. 3.3.5 Effect of plate material

Table 12 presents the first five natural frequencies of a CCCC plate considering MAT1 and MAT2 materials. Results are given considering vibration of a plate in contact with fluid and vacuum vibration. It can be seen that higher natural frequencies are obtained for MAT2, since this material has a higher stiffness. On the other hand, the ratio of wetdry natural frequencies (with fluid and in vacuum) is lower for MAT1, indicating that the effects of the hydroelastic coupling on the vibrational behavior of the plate (as compared to the vacuum vibration) is higher for this material. This occurs because of the lower density of MAT1. Since the solid mass matrix M and fluid mass matrix F are proportional to the plate and fluid density respectively, from Eq. (68) it can be seen that the ratio of wet-dry natural frequencies is highly dependent on the solid-fluid density ratio. As the plate density decreases with the fluid density fixed, the influence of the fluid on the natural frequencies is higher. Since composites usually have lower density than other traditional metallic materials, the hydroelastic coupling has a larger effect on the natural frequencies. In addition, it can be observed that the lamination has negligible effect on the wet-dry frequency ratio. 3.4 Influence of higher-order kinematics

In this section, the influence of higher expansion orders in the accuracy of the numerical results for thin plates is analyzed. Table 13 presents the first five natural frequencies of a CCCC (0/60°) for various side-to-thickness ratios b / h and for materials MAT1 and MAT2. A high expansion order is used as a reference solution, being N = 5 for

b / h = 25 , N = 4 for b / h = 50 , and N = 3 for b / h = 100 . For the case with

b / h = 25 , accurate results for MAT2 are obtained with N = 1 , where the thickness locking correction has been applied. However, for MAT1 higher differences are observed due to the higher orthotropicity ratio E1 / E2 . For b / h = 50 and b / h = 100 the difference between results obtained from different expansion orders is minimal, and expansion orders higher than N = 2 are deemed unnecessary, at least for analyzing the global behavior of the plate.

4. Conclusions An analytical solution for the free vibration analysis of thick rectangular composite plates in contact with a bounded fluid has been presented. The present formulation can analyze thick plates with high accuracy, as observed from comparison with a 3D finite element solution. The influence of geometric parameters and fiber orientation on the natural frequencies has been studied. The conclusions that emerge from this paper can be summarized as follows: a) There is a marked influence of the fluid on the natural frequencies of composite plates, as compared to other more dense materials, and especially on thin plates. b) The convergence of the present formulation is highly dependent on the fluid domain size. c) As the fluid domain size increases, the natural frequencies approach a constant value, corresponding to an infinite fluid domain. d) By using an expansion order N = 4 accurate results are obtained for thick laminated plates, with a maximum difference of 0.72% with respect to the 3D FEM solution for plates with h / b = 0.25 . Even lower errors can be obtained by using a higher expansion order. e) Only certain natural frequencies are dependent on the fluid domain depth. The symmetry and antisymmetry of the mode shapes is seen to influence this dependency. Doubly symmetric mode shapes are dependent of the fluid domain depth, while doubly antisymmetric modes are independent of the fluid depth. f) The influence of the fluid on the natural frequencies is reduced when a wider fluid domain is considered since the fluid motion is less restricted. g) A 2D formulation with higher-order kinematics is necessary to accurately predict the natural frequencies of plates coupled to a fluid domain if composite laminated plates with b / h < 25 are considered. For thinner plates, first-order models yield acceptable results.

References [1] K.H. Jeong, G. Lee, T.W. Kim, Free vibration analysis of a circular plate partially in contact with a liquid, Journal of Sound and Vibration 324 (2009) 194–208. doi:10.1016/j.jsv.2009.01.061.

[2] K.H. Jeong, K.J. Kim, Hydroelastic vibration of a circular plate submerged in a bounded compressible fluid. Journal of Sound and Vibration 283 (2005) 153–172. doi:10.1016/j.jsv.2004.04.029. [3] S. Tariverdilo, M. Shahmardani, J. Mirzapour, R. Shabani, Asymmetric free vibration of circular plate in contact with incompressible fluid. Applied Mathematical Modelling 37 (2013) 228–239. doi:10.1016/j.apm.2012.02.025. [4] C.N. Phan, M. Aureli, M. Porfiri, Finite amplitude vibrations of cantilevers of rectangular cross sections in viscous fluids, Journal of Fluids and Structures 40 (2013) 52–69. doi:10.1016/j.jfluidstructs.2013.03.013. [5] Y. Kozlovsky, Vibration of plates in contact with viscous fluid: Extension of Lamb’s

model.

Journal

of

Sound

and

Vibration

326

(2009)

332–339.

doi:10.1016/j.jsv.2009.04.031. [6] Y. Kerboua, A.A. Lakis, M. Thomas, L. Marcouiller, Vibration analysis of rectangular plates coupled with fluid. Applied Mathematical Modelling 32 (2008) 2570–2586. doi:10.1016/j.apm.2007.09.004. [7] A. Bermudez, L. H. Nieto, R. Rodriguez, Finite element computation of the vibrations of a plate fluid system with interface damping. Computer Methods in Applied Mechanics and Engineering 190 (2001) 3021–3038. [8] S. H. Hashemi, M. Karimi, H. Rokni, Natural frequencies of rectangular Mindlin plates coupled with stationary fluid. Applied Mathematical Modelling 36 (2012) 764– 778. doi:10.1016/j.apm.2011.07.007. [9] D.G. Gorman, J. Horacek, Analysis of the free vibration of a coupled plate/fluid interacting system and interpretation using sub-system modal energy. Engineering Structures 29 (2007) 754–762. doi:10.1016/j.engstruct.2006.06.017. [10] T.P. Chang, On the natural frequency of transversely isotropic magneto-electroelastic plates in contact with fluid. Applied Mathematical Modelling 37 (2013) 2503– 2515. doi:10.1016/j.apm.2012.06.016. [11] S. Carra, M. Amabili, R. Garziera, Experimental study of large amplitude vibrations of a thin plate in contact with sloshing liquids. Journal of Fluids and Structures 42 (2013) 88–111. doi:10.1016/j.jfluidstructs.2013.05.013.

[12] Y.W. Kwon, E.M. Priest, J.H. Gordis, Investigation of vibrational characteristics of composite beams with fluid-structure interaction. Composite Structures 105 (2013) 269–278. doi:10.1016/j.compstruct.2013.05.032. [13] M.R. Kramer, Z. Liu, Y.L. Young, Free vibration of cantilevered composite plates in

air

and

in

water.

Composite

Structures

95

(2013)

254–263.

doi:10.1016/j.compstruct.2012.07.017. [14] Y. K. Cheung, D. Zhou, Coupled vibratory characteristics of a rectangular container bottom plate. Journal of Fluids and Structures 14 (2000) pp. 339-357. [15] D.S. Cho, B.H. Kim, N. Vladimir, T.M. Choi, Natural vibration analysis of rectangular bottom plate structures in contact with fluid. Ocean Engineering 103 (2015) 171–179. doi:10.1016/j.oceaneng.2015.04.078. [16] C.Y. Liao, C.C. Ma, Vibration characteristics of rectangular plate in compressible inviscid

fluid.

Journal

of

Sound

and

Vibration

362

(2016)

228–251.

doi:10.1016/j.jsv.2015.09.031. [17] S. Hosseini-Hashemi, M. Karimi, D.T. Hossein Rokni, Hydroelastic vibration and buckling of rectangular Mindlin plates on Pasternak foundations under linearly varying in-plane loads, Soil Dynamics and Earthquake Engineering 30 (2010) 1487–1499. doi:10.1016/j.soildyn.2010.06.019. [18] A. Shahbaztabar, A.R. Ranji, Effects of in-plane loads on free vibration of symmetrically cross-ply laminated plates resting on Pasternak foundation and coupled with

fluid.

Ocean

Engineering

115

(2016)

196–209.

doi:10.1016/j.oceaneng.2016.02.014. [19] K. Khorshid, S. Farhadi, Free vibration analysis of a laminated composite rectangular plate in contact with a bounded fluid. Composite Structures 104 (2013) 176–186. doi:10.1016/j.compstruct.2013.04.005. [20] D. Seung Cho, B. Hee Kim, J.H. Kim, N. Vladimir, T. Muk Choi, Frequency response of rectangular plate structures in contact with fluid subjected to harmonic point excitation

force.

Thin-Walled

doi:10.1016/j.tws.2015.07.013.

Structures

95

(2015)

276–286.

[21] K.H. Jeong, H.S. Kang, Free vibration of multiple rectangular plates coupled with a liquid.

International

Journal

of

Mechanical

Sciences

74

(2013)

161–172.

doi:10.1016/j.ijmecsci.2013.05.011. [22] E. Carrera, Theories and Finite Elements for Multilayered Plates and Shells: A Unified compact formulation with numerical assessment and benchmarking. Archives of

Computational

Methods

in

Engineering

10

(2003)

215–296.

doi:10.1007/BF02736224. [23] E. Carrera, Transverse normal strain effects on thermal stress analysis of homogeneous and layered plates. AIAA Journal 43 (2005). [24] A. Robaldo, E. Carrera, A. Benjeddou, Unified formulation for finite element thermoelastic analysis of multilayered anisotropic composite plates. Journal of Thermal Stresses 28 (2005) 1031-1065. [25] E. Carrera, Temperature profile influence on layered plates response considering classical and advanced theories. AIAA Journal 40 (2002). [26] E. Carrera, M. Boscolo, A. Robaldo, Hierarchic multilayered plate elements for coupled multifield problems of piezoelectric adaptive structures: Formulation and numerical assessment. Archives of Computational Methods in Engineering 14 (2007) 383–430. doi:10.1007/s11831-007-9012-8. [27] E. Carrera, S. Brischetto, P. Nali, Variational Statements and Computational Models for MultiField Problems and Multilayered Structures. Mechanics of Advanced Materials and Structures 15 (2008) 182–198. doi:10.1080/15376490801907657. [28] E. Carrera, S. Brischetto, C. Fagiano, P. Nali, Mixed multilayered plate elements for coupled magneto-electro-elastic analysis. Multidiscipline Modeling in Materials and Structures 5 (2009) pp. 251-256. [29] A. Robaldo, E. Carrera, A. Benjeddou, A unified formulation for finite element analysis of piezoelectric adaptive plates. Computers and Structures 84 (2006) 1494– 1505. doi:10.1016/j.compstruc.2006.01.029. [30] M. Cinefra, E. Carrera, L. Della Croce, C. Chinosi, Refined shell elements for the analysis of functionally graded structures. Composite Structures 94 (2012) 415–422. doi:10.1016/j.compstruct.2011.08.006.

[31] M. Cinefra, C. Chinosi, L. Della Croce, MITC9 shell elements based on refined theories for the analysis of isotropic cylindrical structures. Mechanics of Advanced Materials and Structures 20 (2013) 91–100. doi:10.1080/15376494.2011.581417. [32] E. Carrera, S. Brischetto, A. Robaldo, Variable Kinematic Model for the Analysis of Functionally Graded Material Plates, AIAA Journal 46 (2008) 194–203. doi:10.2514/1.32490. [33] E. Carrera, M. Filippi, E. Zappino, Laminated beam analysis by polynomial, trigonometric, exponential and zig-zag theories, European Journal of Mechanics A/Solids. 41 (2013) 58–69. doi:10.1016/j.euromechsol.2013.02.006. [34] M. Filippi, M. Petrolo, S. Valvano, E. Carrera, Analysis of laminated composites and sandwich structures by trigonometric, exponential and miscellaneous polynomials and a MITC9 plate element, Composite Structures, 150 (2016) 103-114. [35] J.L. Mantari, I.A. Ramos, E. Carrera, M. Petrolo, Static analysis of functionally graded plates using new non-polynomial displacement fields via Carrera Unified Formulation,

Composites

Part

B

Eng.

89

(2016)

127–142.

doi:10.1016/j.compositesb.2015.11.025. [36] I.A. Ramos, J.L. Mantari, A.M. Zenkour, Laminated composite plates subject to thermal load using trigonometrical theory based on Carrera Unified Formulation, Composite Structures 143 (2016) 324–335. doi:10.1016/j.compstruct.2016.02.020. [37] E. Carrera, G. Giunta, M. Petrolo, Beam Structures: Classical and Advanced Theories. Wiley, Chichester, United Kingdom, 2011. [38] E. Carrera, S. Brischetto, P. Nali, Plates and Shells for Smart Structures: Classical and Advanced Theories for Modeling and Analysis. Wiley, Chichester, United Kingdom, 2011. [39] E. Carrera, M. Cinefra, M. Petrolo, E. Zappino, Finite Element Analysis of Structures through Unified Formulation. Wiley, Chichester, United Kingdom, 2014. [40] A.W. Leissa, The historical bases of the Rayleigh and Ritz methods. Journal of Sound and Vibration 287 (2005) 961–978. doi:10.1016/j.jsv.2004.12.021.

[41] M.J. Gander, G. Wanner, From Euler, Ritz, and Galerkin to Modern Computing. SIAM Review 54 (2012) 627–666. doi:10.1137/100804036. [42] S. Ilanko, L. Monterrubio, Y. Mochida, The Rayleigh-Ritz Method for Structural Analysis. Wiley, London, United Kingdom, 2014. [43] F.A. Fazzolari, E. Carrera, Accurate free vibration analysis of thermo-mechanically pre/post-buckled anisotropic multilayered plates based on a refined hierarchical trigonometric

Ritz

formulation.

Composite

Structures

95

(2013)

381–402.

doi:10.1016/j.compstruct.2012.07.036. [44] F.A. Fazzolari, E. Carrera, Advanced variable kinematics Ritz and Galerkin formulations for accurate buckling and vibration analysis of anisotropic laminated composite

plates.

Composite

Structures

94

(2011)

50–67.

doi:10.1016/j.compstruct.2011.07.018. [45] F.A. Fazzolari, E. Carrera, Coupled thermoelastic effect in free vibration analysis of anisotropic multilayered plates and FGM plates by using a variable-kinematics Ritz formulation. European Journal of Mechanics A/Solids. 44 (2014) 157–174. doi:10.1016/j.euromechsol.2013.10.011. [46] F.A. Fazzolari, E. Carrera, Free vibration analysis of sandwich plates with anisotropic face sheets in thermal environment by using the hierarchical trigonometric Ritz

formulation.

Composites

Part

B:

Engineering

50

(2013)

67–81.

doi:10.1016/j.compositesb.2013.01.020. [47] F.A. Fazzolari, Reissner’s Mixed Variational Theorem and variable kinematics in the modelling of laminated composite and FGM doubly-curved shells. Composites Part B: Engineering 89 (2016) 408–423. doi:10.1016/j.compositesb.2015.11.031. [48] F.A. Fazzolari, J.R. Banerjee, Axiomatic/asymptotic PVD/RMVT-based shell theories for free vibrations of anisotropic shells using an advanced Ritz formulation and accurate

curvature

descriptions.

Composite

Structures

108

(2014)

91–110.

doi:10.1016/j.compstruct.2013.08.037. [49] F.A. Fazzolari, E. Carrera, Advances in the Ritz formulation for free vibration response of doubly-curved anisotropic laminated composite shallow and deep shells. Composite Structures 101 (2013) 111–128. doi:10.1016/j.compstruct.2013.01.018.

[50] L. Dozio, Natural frequencies of sandwich plates with FGM core via variablekinematic

2-D

Ritz

models.

Composite

Structures

96

(2013)

561–568.

doi:10.1016/j.compstruct.2012.08.016. [51] L. Dozio, On the use of the Trigonometric Ritz method for general vibration analysis of rectangular Kirchhoff plates. Thin-Walled Structures 49 (2011) 129–144. doi:10.1016/j.tws.2010.08.014. [52] L. Dozio, In-plane free vibrations of single-layer and symmetrically laminated rectangular

composite

plates.

Composite

Structures

93

(2011)

1787–1800.

doi:10.1016/j.compstruct.2011.01.021. [53] L. Dozio, Free in-plane vibration analysis of rectangular plates with arbitrary elastic boundaries. Mechanics Research Communications 37 (2010) 627–635. doi:10.1016/j.mechrescom.2010.09.003. [54] R. Vescovini, L. Dozio, A variable-kinematic model for variable stiffness plates: Vibration and buckling analysis. Composite Structures 142 (2016) 15–26. doi:10.1016/j.compstruct.2016.01.068. [55] L. Dozio, E. Carrera, A variable kinematic Ritz formulation for vibration study of quadrilateral plates with arbitrary thickness. Journal of Sound and Vibration 330 (2011) 4611–4632. doi:10.1016/j.jsv.2011.04.022. [56] I. M. Daniel, O. Ishai, Engineering mechanics of composite materials, Oxford University Press, New York, 2006. [57] J. N. Reddy, Mechanics of laminated composite plates and shells, CRC Press, Boca Raton, 2004. [58] E. Carrera, S. Brischetto, Analysis of thickness locking in classical, refined and mixed multilayered plate theories, Composite Structures 82 (2008) 549-562.

Tables caption Table 1. Material properties. Table 2: Convergence of the natural frequencies of a square CCCC (0/90°) MAT1 plate with respect to the number of terms used in the Ritz series for various CUF expansion orders, considering p = q = 16 , a = b = c = d = e = 1 m , h / b = 0.1 . Table 3: Convergence of the natural frequencies of a square CCCC (0/90°) MAT1 plate with respect to the number of terms used in the trigonometric series for various fluid domain – plate width ratios, considering N = 5 , a = b = e = 1 m , h / b = 0.1 , c = d . Table 4. Natural frequencies of a CCCC (0/60°) MAT1 and MAT2 plate in contact with fluid considering h / b = 0.1 . Table 5. Natural frequencies of a CCCC (0/60°) MAT1 and MAT2 plate in contact with fluid considering h / b = 0.25 . Table 6. Nondimensional natural frequencies of an isotropic plate in contact with a fluid domain considering h / b = 0.05 , ν = 0.3 , a / b = c / d = d / b = e / b = 1 Table 7: Natural frequencies of a CCCC (0/θ°) MAT1 plate in contact with fluid for various plate side-to-thickness ratios. Table 8: Natural frequencies of a CCCC (0/θ°) MAT1 plate in contact with fluid for various plate aspect ratios. Table 9: Natural frequencies of a CCCC (0/θ°) MAT1 plate in contact with fluid for various fluid domain depth ratios. Table 10: Natural frequencies of a CCCC (0/θ°) MAT1 plate in contact with fluid for various fluid domain width – plate width ratios. Table 11: Natural frequencies of a (0/θ°) MAT1 plate in contact with fluid for various plate boundary conditions. Table 12: Natural frequencies of a CCCC plate in contact with fluid and in vacuum for plate materials MAT1 and MAT2. Table 13: Natural frequencies of a CCCC (0/60°) plate considering MAT1 and MAT2 materials, and a = b = c = d = e = 1 m .

Figures caption Figure 1. Coordinate frame of the plate. Figure 2. Coordinate frame of the fluid domain and the plate. Note that O and o ' are located at z
= −h / 2 or z = 0 . Figure 3. Finite element model used in ANSYS software. Figure 4. Mode shapes of vibration for a square CCCC (0/60°) MAT1 plate with

h / b = 0.1 . Figure 5: Effect of plate thickness on the natural frequencies of a square CCCC MAT1 plate in contact with fluid using N = 5 . Figure 6: Effect of plate aspect ratio on the natural frequencies of a square CCCC MAT1 plate in contact with fluid using N = 5 . Figure 7. Mode shapes of vibration for a square CCCC (0/90°) MAT1 plate with

e / b = 3.

Tables Table 1.

Material property

E1 (GPa) E2 (GPa) E3 (GPa) G12 (GPa) G13 (GPa) G23 (GPa) ν 12 ν 13 ν 23 ρ (kg/m3)

Material MAT1 (Graphite-Epoxy)

MAT2 (Boron-Epoxy)

40

206.9

1

20.69

1

20.69

0.6

6.9

0.6

6.9

0.5

4.14

0.25

0.3

0.25

0.3

0.25

0.3

1000

1950

Table 2.

Mode

Theory

P

N=1

ω1

ω2

ω3

ω4

ω5

6 7 8 9 10 11 12 13

114.24 114.24 114.24 114.24 114.24 114.24 114.23 114.23

308.17 308.14 308.14 308.13 308.13 308.13 308.13 308.13

309.34 309.30 309.30 309.29 309.29 309.29 309.29 309.29

452.06 451.97 451.91 451.91 451.90 451.90 451.89 451.89

515.71 515.61 515.49 515.49 515.47 515.47 515.46 515.46

N=3

6 7 8 9 10 11 12 13

108.66 108.61 108.58 108.57 108.56 108.56 108.56 108.56

289.94 289.74 289.67 289.60 289.59 289.57 289.57 289.56

289.94 289.75 289.67 289.61 289.60 289.58 289.58 289.57

425.47 424.83 424.57 424.47 424.38 424.36 424.34 424.33

481.94 481.42 481.03 480.96 480.87 480.86 480.83 480.83

N=5

6 7 8 9 10 11 12 13

106.16 106.10 106.07 106.04 106.02 106.01 106.00 105.99

282.55 282.35 282.22 282.16 282.09 282.06 282.04 282.02

282.82 282.63 282.47 282.39 282.31 282.28 282.24 282.22

415.60 415.00 414.75 414.53 414.44 414.33 414.29 414.24

470.85 470.22 469.83 469.61 469.51 469.42 469.36 469.32

Table 3.

d /b

p,q

3

Mode

ω1

ω2

ω3

ω4

ω5

6 10 14 18 22

140.43 140.01 139.96 139.96 139.96

332.84 304.55 303.81 303.66 303.64

332.89 304.77 304.07 303.91 303.89

515.06 439.40 437.43 436.93 436.87

621.72 523.40 513.53 513.08 513.00

6

8 16 24 32 40 48

147.89 140.80 140.74 140.72 140.72 140.72

392.13 311.33 304.13 303.96 303.92 303.90

392.17 311.43 304.39 304.22 304.17 304.15

609.39 458.18 437.55 437.09 436.96 436.90

707.30 570.85 516.83 514.60 514.13 514.07

12

16 32 48 64 80 96

149.65 140.84 140.75 140.73 140.73 140.73

393.63 311.48 304.14 303.97 303.92 303.90

393.65 311.58 304.40 304.22 304.18 304.16

609.34 458.50 437.56 437.09 436.96 436.90

710.89 578.48 517.76 514.63 514.15 514.09

Table 4.

Theory

Mode

ω1

FEM 3D 67.158 PM N = 1 TL 71.641 PM N = 1 71.191 PM N = 2 69.414 PM N = 3 68.513 PM N = 4 67.622 PM N = 5 67.243

FEM 3D 197.61 PM N = 1 TL 205.67 PM N = 1 201.52 PM N = 2 202.49 PM N = 3 198.69 PM N = 4 198.46 PM N = 5 198.04

ω5

Max diff Avg diff (%) (%)

133.51 143.60 142.63 138.39 136.06 134.10 133.42

MAT1 151.76 203.28 164.07 218.73 163.46 217.39 157.89 210.82 155.10 206.89 152.90 204.20 151.82 203.08

237.31 256.51 255.25 246.89 242.02 238.41 236.97

(8.11) (7.71) (4.04) (2.20) (0.75) (0.14)

(7.61) (7.01) (3.76) (1.98) (0.56) (0.09)

390.32 407.68 400.08 402.00 392.34 391.86 390.89

MAT2 431.56 582.28 449.96 608.16 443.64 598.36 445.20 601.09 434.61 585.72 433.78 585.04 432.29 583.17

669.71 699.45 689.14 692.35 673.87 672.70 670.40

(4.45) (2.90) (3.38) (0.71) (0.51) (0.22)

(4.33) (2.59) (3.05) (0.60) (0.45) (0.16)

ω2

ω3

ω4

Table 5.

Theory

Mode

ω1

FEM 3D 127.77 PM N = 1 TL 139.40 PM N = 1 139.03 PM N = 2 133.42 PM N = 3 130.00 PM N = 4 128.62 PM N = 5 127.80

FEM 3D 366.85 PM N = 1 TL 385.09 PM N = 1 382.01 PM N = 2 382.16 PM N = 3 370.40 PM N = 4 369.68 PM N = 5 367.66

ω2

ω3

ω4

ω5

Max diff Avg diff (%) (%)

228.45 247.05 246.30 238.02 231.23 229.51 228.39

MAT1 243.42 264.71 264.34 253.77 246.80 245.16 243.50

320.26 346.80 345.88 333.53 323.55 321.71 320.16

358.86 390.68 389.84 374.79 364.75 361.37 358.94

(9.10) (8.81) (4.44) (1.74) (0.72) (0.03)

(8.63) (8.37) (4.29) (1.40) (0.60) (0.03)

630.02 658.91 653.85 655.99 636.23 635.03 631.32

MAT2 675.35 705.46 701.71 702.32 682.97 681.66 676.79

872.33 911.18 905.38 907.86 881.77 879.99 874.15

979.41 1024.6 1018.1 1019.4 990.58 988.13 981.12

(4.97) (4.13) (4.17) (1.14) (0.93) (0.22)

(4.62) (3.91) (4.09) (1.06) (0.85) (0.20)

ω1

ω2

ω3

ω4

ω5

ω6

Table 6.

Boundary condition

Theory

Mode

SSSS

CPT [14] MPT [17] PM N = 1

11.844 11.752 11.765

39.243 38.449 38.557

39.243 38.449 38.557

66.609 64.462 64.743

81.643 78.408 78.823

85.674 82.247 82.686

CCCC

CPT [17] MPT [17] PM N = 1

23.097 22.398 22.500

60.293 57.283 57.689

60.293 57.283 57.689

93.645 87.312 88.096

108.65 100.35 101.42

115.96 107.23 108.36

Table 7.

Lamination angle (0 / θ°)

Mode

ω1

ω2

ω3

ω4

ω5

0.04

0 15 30 45 60 75 90

53.038 46.871 40.225 37.137 36.238 36.368 36.553

99.851 94.941 91.132 94.655 102.61 110.51 115.22

151.98 150.34 149.42 133.28 123.55 118.38 115.34

194.04 174.21 153.22 163.06 171.63 176.94 179.13

230.48 220.21 201.72 197.99 207.95 216.60 217.30

0.1

0 15 30 45 60 75 90

121.51 116.85 109.95 105.94 104.93 105.54 106.01

235.99 235.02 237.43 247.36 261.41 273.97 282.06

361.68 349.75 329.23 310.86 296.97 288.23 282.28

369.92 371.63 379.28 390.54 401.04 410.12 414.33

444.78 440.93 438.49 447.67 464.10 470.88 469.41

0.25

0 15 30 45 60 75 90

208.47 207.45 205.64 204.63 204.86 205.73 206.20

417.22 419.36 426.21 436.92 448.83 458.46 462.57

517.77 512.01 500.56 488.28 477.78 470.91 468.27

623.16 622.90 623.52 626.38 632.23 639.93 644.53

632.57 654.78 667.21 683.30 691.56 686.09 680.73

h/b

Table 8.

Lamination angle (0 / θ°)

Mode

ω1

ω2

ω3

ω4

ω5

2

0 15 30 45 60 75 90

58.425 57.201 57.716 62.459 69.175 74.617 76.658

133.34 126.80 118.48 115.16 115.91 118.11 119.20

162.76 166.71 180.29 193.04 187.74 185.92 185.63

226.51 218.22 203.27 203.10 227.70 245.85 252.46

239.43 236.73 239.18 249.72 265.05 276.17 274.88

3

0 15 30 45 60 75 90

43.241 43.787 47.317 54.763 63.326 69.702 72.016

76.909 74.767 73.734 77.065 83.072 88.355 90.396

132.00 124.94 116.61 113.99 115.76 118.82 120.23

149.44 153.42 168.47 163.68 160.41 160.29 160.64

176.07 180.75 172.98 194.32 216.14 213.09 212.38

4

0 15 30 45 60 75 90

38.433 39.407 43.780 52.065 61.182 67.839 70.239

55.931 55.946 58.309 64.529 72.451 78.647 80.938

87.755 84.818 82.716 85.205 90.863 96.144 98.224

131.86 124.43 115.98 113.62 115.80 119.21 120.74

145.48 149.23 157.02 149.54 147.73 148.67 149.40

5

0 15 30 45 60 75 90

36.499 37.528 42.140 50.747 60.091 66.873 69.314

46.665 47.554 51.370 58.878 67.578 74.106 76.484

66.281 65.630 66.925 72.169 79.654 85.803 88.119

95.132 91.469 88.461 90.321 95.660 100.86 102.94

131.74 124.10 115.60 113.41 115.82 119.43 121.04

a/b

Table 9.

Lamination angle (0 / θ°)

Mode

ω1

ω2

ω3

ω4

ω5

2

0 15 30 45 60 75 90

90.418 87.078 82.108 79.208 78.478 78.915 79.255

235.76 234.80 237.21 247.13 261.16 273.70 281.79

361.29 349.38 328.89 310.54 296.67 287.94 282.00

366.34 367.72 374.89 386.46 398.29 409.05 414.31

444.75 440.86 438.06 445.64 458.51 461.36 458.58

3

0 15 30 45 60 75 90

75.155 72.419 68.341 65.961 65.361 65.721 66.001

235.76 234.80 237.21 247.12 261.16 273.70 281.79

361.29 349.38 328.89 310.54 296.67 287.93 282.00

365.07 366.32 373.28 384.89 397.16 408.58 414.31

444.75 440.85 437.92 444.94 456.53 458.00 454.68

4

0 15 30 45 60 75 90

65.680 63.307 59.769 57.703 57.182 57.494 57.738

235.76 234.80 237.21 247.12 261.16 273.70 281.79

361.29 349.38 328.89 310.54 296.67 287.93 282.00

364.42 365.60 372.44 384.06 396.55 408.31 414.31

444.75 440.84 437.85 444.58 455.52 456.30 452.67

5

0 15 30 45 60 75 90

59.070 56.946 53.778 51.928 51.461 51.741 51.960

235.77 234.80 237.21 247.12 261.16 273.70 281.79

361.29 349.38 328.89 310.54 296.67 287.93 282.00

364.02 365.16 371.93 383.55 396.17 408.14 414.31

444.74 440.84 437.80 444.37 454.91 455.27 451.45

e/b

Table 10.

Lamination Mode angle (0 / θ°) ω

ω2

ω3

ω4

ω5

1.4

0 15 30 45 60 75 90

143.93 138.18 129.72 124.82 123.60 124.34 124.90

247.56 246.46 248.86 259.28 274.18 287.70 296.73

381.53 369.99 347.77 328.07 313.22 303.72 296.97

383.13 383.82 393.09 406.60 418.70 428.59 433.08

465.79 461.25 457.57 466.13 484.31 494.17 492.95

1.8

0 15 30 45 60 75 90

154.21 147.94 138.74 133.41 132.10 132.91 133.52

250.92 249.79 252.19 262.75 277.91 291.71 301.02

385.84 375.85 353.15 333.08 317.96 308.25 301.26

389.35 388.31 397.89 411.33 422.74 431.84 435.90

468.92 464.30 460.66 470.07 490.50 504.16 503.95

2.2

0 15 30 45 60 75 90

158.83 152.32 142.78 137.26 135.90 136.75 137.38

252.20 251.05 253.45 264.07 279.32 293.24 302.65

387.80 378.08 355.20 334.99 319.76 309.97 302.90

391.71 390.34 399.99 413.21 424.13 432.76 436.58

469.69 465.06 461.52 471.46 493.04 508.80 509.25

2.6

0 15 30 45 60 75 90

160.93 154.31 144.61 139.01 137.64 138.50 139.14

252.75 251.59 253.99 264.63 279.93 293.89 303.34

388.74 379.04 356.07 335.81 320.53 310.71 303.60

392.73 391.30 400.97 414.06 424.72 433.12 436.83

469.97 465.35 461.87 472.09 494.22 511.01 511.84

d /b

1

Table 11.

Boundary Lamination condition angle (0 / θ°)

Mode

ω1

ω2

ω3

ω4

ω5

SSSS

0 15 30 45 60 75 90

82.001 79.166 73.598 67.071 63.165 59.788 58.372

167.39 163.52 172.83 185.96 198.60 209.69 217.08

297.03 299.79 276.51 255.92 239.56 226.46 217.74

330.70 301.38 315.54 330.91 339.81 341.25 340.29

387.30 389.56 383.35 393.61 416.46 416.22 387.30

CSCS

0 15 30 45 60 75 90

112.78 107.59 99.657 92.857 87.869 85.149 84.411

203.87 202.80 203.75 208.86 219.19 230.71 236.24

324.23 328.38 319.64 299.83 282.19 269.55 264.56

353.22 338.74 340.29 355.71 367.29 374.55 377.74

387.30 416.35 412.27 415.66 430.41 443.86 387.30

CCSC

0 15 30 45 60 75 90

106.09 102.18 96.947 94.870 94.804 94.554 94.380

218.76 218.95 224.50 237.74 252.55 260.25 260.09

350.51 336.88 313.76 294.15 279.93 272.68 272.65

357.66 359.74 368.78 380.19 389.92 395.30 396.86

434.39 429.44 426.56 438.76 457.52 462.52 459.46

SSCS

0 15 30 45 60 75 90

97.240 92.440 85.366 79.014 74.747 71.641 70.479

184.99 185.63 190.18 198.27 209.44 220.30 225.67

309.82 314.59 299.86 278.36 260.81 247.63 241.45

342.63 322.91 329.25 344.43 354.20 358.16 359.18

387.30 403.42 397.54 404.33 423.29 437.01 387.30

Table 12.

Layup

Model

Mode 1

2

3

4

5

414.33 653.98 0.634

469.41 778.35 0.603

0/90°

With fluid In vacuum Ratio

106.01 286.17 0.370

MAT1 282.06 282.28 499.14 499.14 0.565 0.566

45°/-45°

With fluid In vacuum Ratio

101.90 275.25 0.370

275.55 488.09 0.565

276.18 488.09 0.566

428.02 679.24 0.630

456.25 746.74 0.611

1168.8 1544.8 0.757

1327.2 1853.5 0.716

1211.5 1608.5 0.753

1282.3 1767.0 0.726

0/90°

With fluid In vacuum Ratio

321.74 656.88 0.490

MAT2 813.68 815.82 1172.7 1172.7 0.694 0.696

45°/-45°

With fluid In vacuum Ratio

309.40 633.00 0.489

791.72 1141.6 0.693

793.74 1141.6 0.695

Table 13.

b/h

Theory

Mode

ω1

25 PM PM PM PM PM

N=1 N=2 N=3 N=4 N=5

36.865 36.636 36.501 36.312 36.238

PM PM PM PM PM

N=1 N=2 N=3 N=4 N=5

108.95 109.71 109.08 109.04 108.98

50 PM PM PM PM

N=1 N=2 N=3 N=4

13.848 13.852 13.838 13.814

PM PM PM PM

N=1 N=2 N=3 N=4

41.225 41.540 41.470 41.464

100 PM N = 1 PM N = 2 PM N = 3

4.9966 5.0087 5.0074

PM N = 1 PM N = 2 PM N = 3

14.960 15.078 15.071

ω2

ω3

MAT1 105.21 127.36 104.17 125.77 103.64 125.03 102.86 123.99 102.61 123.55 MAT2 309.68 352.76 311.63 355.06 309.05 351.96 308.90 351.68 308.73 351.38 MAT1 41.403 50.887 41.365 50.778 41.304 50.690 41.194 50.535 MAT2 125.07 143.41 125.93 144.48 125.61 144.09 125.58 144.05 MAT1 15.258 18.829 15.292 18.857 15.286 18.849 MAT2 46.913 53.868 47.241 54.286 47.207 54.245

Avg diff (%)

ω4

ω5

176.90 174.75 173.60 172.15 171.63

215.63 212.40 210.84 208.69 207.96

(2.82) (1.67) (1.09) (0.29)

512.06 514.89 509.57 509.32 508.90

602.97 606.96 599.99 599.51 598.94

(0.41) (1.03) (0.13) (0.07)

72.334 72.160 72.014 71.782

89.487 89.200 88.998 88.636

(0.63) (0.47) (0.30)

215.80 217.05 216.31 216.27

257.94 259.83 258.83 258.75

(0.39) (0.31) (0.02)

27.116 27.151 27.137

33.721 33.772 33.751

(0.13) (0.05)

82.730 83.219 83.140

99.539 100.29 100.18

(0.64) (0.08)

Figures

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Mode no.

1

2

3

4

5

Present method

3D FEM

Figure 5.

Figure 6.

Figure 7.

Mode no.

1

2

3

4

5

Fiber orientation 0/0°

0/90°