Nonlinear random responses and fatigue prediction of elastically restrained laminated composite panels in thermo-acoustic environments

Composite Structures 229 (2019) 111391

Contents lists available at ScienceDirect

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

Nonlinear random responses and fatigue prediction of elastically restrained laminated composite panels in thermo-acoustic environments

T

Yilong Wanga, , Dengqing Caoa, , Jiaqi Penga, Hao Chengb, Huagang Lina, Wenhu Huanga ⁎

a b

⁎

School of Astronautics, Harbin Institute of Technology, Harbin 150001, PR China Beijing Institute of Structure and Environment Engineering, Beijing 100076, PR China

ARTICLE INFO

ABSTRACT

Keywords: Laminated composite panel Elastic boundary condition Thermo-acoustic vibration Rayleigh-Ritz method

This paper presents a formulation for predicting the nonlinear random response of the elastically restrained laminated composite panel subjected to thermo-acoustic loads. Based on the laminated plate theory and Von Kármán large deflection and classical thin plate theories, the natural characteristics are obtained via RayleighRitz method and then the governing equations of the panel subjected to combined acoustic and thermal loads are formulated. The nonlinear partial differential equations of motion are transformed to a set of coupled nonlinear ordinary differential equations in truncated modal coordinates. A numerical example where the acoustic load is considered as the Gaussian band-limited white noise is given to perform the process of obtaining the mode and responses of the panel. Taking the natural frequency obtained from the finite element method as a reference value, the process of obtaining the natural frequencies is validated by comparing the frequency results. Numerical results show that the buckling, snap-through, and nonlinear random vibrations of the thermal-elastic restrained panel can be predicted accurately. Comparing stress PSD distributions with fatigue damage distributions, the first-order mode is proved to be valid for determining the most dangerous area for fatigue life prediction.

1. Introduction Skin plates are a commonly used form of structural components especially in aerospace vehicles, such as high-speed aircraft, hypersonic flights, rockets, and spacecraft, which are subjected to thermal loads due to aerodynamic and/or solar radiation heating, and random acoustic loads due to engine and/or aerodynamic transonic noise [1,2]. These complicated working environments may lead to elevated temperatures and pressure distributions over the surface of the panel. Additionally, the presence of these thermal fields and acoustic loads may also cause large thermal deflections (thermal buckling) and thermalacoustic fatigue of the skin panels. In general, thermal buckling does not indicate structural failure [3]. However, the thermal environment is still a cruel factor which results in thermal stresses internal of the structures, modifies the stiffness of the structural system and alters the dynamic characteristics of the system essentially when subjected to acoustic and/or aerodynamic loads [4], consequently resulting in a reduction in the flight performance. Accordingly, it is significant to consider the interactive effects of thermal and acoustic loads. In recent years, the thermo-acoustic random responses of the skin panels of aerospace vehicles at elevated temperatures and intensity acoustic

⁎

loads have become a matter of considerable importance. The range of the engine acoustic sound pressure levels for hypersonic vehicles has been studied to lie between 130 dB and 180 dB [5–7], and the fluctuation of the operating temperatures for typical skin panels is often from 1000 °F to 3000 °F [8]. It is obvious if the skin panel work in such severe environments, there will be quite strong nonlinear characteristics which can reduce the fatigue life of the skin panel even directly lead to failures of the structure. Extensive research work on panels under thermo-acoustic loads has been carried out with the fast development of Aeronautics and Astronautics. For the prediction of the nonlinear response of structures, there are a few major analysis methods, such as perturbation, Fokker–Plank–Kolmogorov (FPK equation), Monte Carlo, equivalent linearization (EL), finite element (FE) numerical integration, etc. [9]. For systems with strong geometric nonlinearities and complex structures, all of these methods are limited or unavailable to some extent. But, with the fast development of computer hardware, the finite element method (FEM) has been becoming the most widely used method, the application of which has been extended to structures subjected to thermo-acoustic loads since many years ago [10,11]. Assuming temperature-independent material properties, Locke [12] studied the large

Corresponding authors. E-mail addresses: [email protected] (Y. Wang), [email protected] (D. Cao).

https://doi.org/10.1016/j.compstruct.2019.111391 Received 17 June 2019; Received in revised form 12 August 2019; Accepted 9 September 2019 Available online 12 September 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 229 (2019) 111391

Y. Wang, et al.

Nomenclature

a b h k th zk 0

A D T0 T 1

2

w

U UT Us k xt 0 (y )

Restraints along the edge x = 0 and x = a Stiffness functions for the translational elastic Restraints along the edge y = 0 and y = b Stiffness functions for the rotational elastic Restraints along the edge x = 0 and x = a Stiffness functions for the translational elastic Restraints along the edge y = 0 and y = b Thermal stress function along x-axis Thermal stress function along y-axis Thermal stress function along xy-plane Density of the composite lamina Circular frequency W (x , y ) Mode function amn Unknown coefficient obtained by the Rayleigh-Ritz method Orthogonal polynomials along x-axis m (x ) Orthogonal polynomials along y-axis n (y ) Wi (x , y ) i-th mode function N Total amount of the modal coordinates p (x , y , t ) Sound pressure on the surface of the panel c0 Structural damping coefficient of the panel FT Vector of thermal loads N (q) Nonlinear stiffness matrix t k xa (y ) k yt 0 (x ) t k yb (x ) k xr 0 (y ) r k xa (y ) k yr 0 (x ) r k yb (x ) NxT NyT T Nxy

Length of the panel Width of the panel Thickness of the panel Number of the lamina Thickness of the k-th lamina Membrane nonlinear strain vector Coefficient vector of the bending strain Membrane stiffness matrix Bending stiffness matrix =The initial temperature of the panel The current temperature of the panel The thermal expansion coefficient along the 1 main direction of the material The thermal expansion coefficient along the 2 main direction of the material The transverse displacement measured from the mid-surface of the panel Deformation potential energy of the panel Thermal potential energy of the panel Potential energy resulting from the deflections of the elastic constraints Stiffness functions for the translational elastic

deflection random vibration of a thermally buckled thin isotropic plate under a stationary Gaussian pressure with zero mean and uniform magnitude and phase over the plate surface, using the method of equivalent linearization, while Dhainaut et al. [9] presented a finite element formulation for the prediction of nonlinear random response of isotropic and composite panels simultaneously subjected to high acoustic loads and elevated temperatures. Przekop et al. [13] and Radu et al. [14] investigated the random response and fatigue life of panels subjected to thermo-acoustic loads by numerical simulation, also using the finite element method. Przekop and Rizzi [15] studied the nonlinear snap-through characteristics of thin-wall structures subjected to combined thermo-acoustic loads. Ibrahim et al. [3] provided a time-domain finite element formulation to investigate the nonlinear random response of composite plates impregnated with pre-strained SMA fibers subject to thermal and random acoustic loads. Shukla et al. [16] and Sha et al. [17–20] pointed out that thin-walled structures will exhibit a complex nonlinear response under the combined thermo-acoustic loadings, including nonlinear random vibration around the initial equilibrium position, snap-through between post-buckled equilibrium positions and nonlinear vibration around one post-buckled equilibrium position. The influences of thermal-acoustic excitations on nonlinear dynamics response for panels and curved panels were then analyzed in Sha et al. [21,22], and the corresponding multi-axial fatigue life estimation was also presented in their study. There are also experiments performed to study thermally loaded panels under random excitation by Raymond Istenes et al. [23], Ng and Clevenson [24] and Murphy et al. [25], where the snap-through phenomenon and frequency shifting due to nonlinear large-amplitude vibration were observed. The probabilistic solutions of a nonlinear plate were investigated in Er and Iu [26] by employing the state space-split method and exponential polynomial closure method, where the plate is excited by uniformly distributed Gaussian white noise which is fully correlated in space. A relatively rigor verification of the proposed method was also given in [26] via the comparison of the computational effort and numerical results with those obtained by Monte Carlo simulation and EL, respectively. Although considerable attention has devoted to researches on the nonlinear random response of panels subjected to thermo-acoustic loads, the boundary conditions of these panels were mainly considered classical, such as simply supported boundary or clamped boundary, etc.

However, few researchers study the elastically restrained panels subjected to the thermo-acoustic loads. When it comes to the background of the aerospace vehicle, where the panel made by composite materials, such as the functionally graded [1,27] and the curvilinear fiber laminated composites [28–30], etc., are widely used, the related investigation is much less. Therefore, the investigation on the response of the LCP under combined thermal and random acoustic loads is of great significance. The purpose of this paper is to explore the method of establishing the model of the elastically restrained LCP and study its nonlinear random response when subjected to acoustic loads and elevated temperatures. A LCP with the elastic restraint boundary is presented, and the Rayleigh-Ritz method using the boundary characteristic orthogonal polynomials is employed to obtain the modes of the system under thermal loads. The governing equations of motion are derived using the Hamilton principle, which are then reduced to a set of coupled nonlinear modal equations. The band-limited Gaussian white noise is chosen as the random acoustic load and is considered to act on the panel with an angle, meaning that the sound pressure is not uniformly distributed over the surface. The numerical simulation is performed to obtain the modes for validating the formulation, then the responses under different combinations of sound pressure level (SPL) and temperature environment, and the PSD of the response, where three types of motion can be predicted. The stress/strain PSD distributions and the damage distributions are also presented to determine the most dangerous locations on the panel for fatigue life prediction and strength evaluation. 2. Formulation The LCP elastically restrained along edges (x = 0, a and y = 0, b) are shown in Fig. 1. The linear and rotational springs are evenly distributed along the edges of the panel (for clarity, only several springs are shown here, but actually there are a number of springs along the edges.). The total kinetic energy is simply calculated from

T=

2

1 h 2 Area

w t

2

dxdy ,

(1)

Composite Structures 229 (2019) 111391

Y. Wang, et al.

U=U +U +

1 2

+

1 2

0

1 2

0

1 2

0

+ +

+

Fig. 1. The schematic of an LCP elastically restrained along edges.

T

=

y

=

T,

y

T=T

¯k =Q

y

T0,

xy

xy

x

y

y

xy

k

¯ k( =Q

k

+

T T [NxT , NyT , Nxy ] =

H /2

¯k

Q

T

=

1 2

Area

NxT

w x

2

(3)

T dz .

+ NyT

2

w y

2

w y

2

w w x y

dxdy

dy x=0

dy x=a

dx y=0

dx y=b

(6)

(7)

T + Nxy

w w x y

n

y . b

(8)

U ,max

Us,max aij

U

T ,max )

= 0.

(K s + K + K

2 M) a

T

(10)

= 0,

where Ks , K , K T and M are all (um ·un) × (um ·un) matrices and a is a y x vector of (um ·un) dimensions. Letting = a and = b , the entries of these matrices and vector can be expressed as (only containing the linear part)

k = (i

1) × un + j,

l = (m

1) × um + n,

(i, m = 1, 2, 3 um ; j, n = 1, 2, 3 un)

dxdy.

(9)

The minimizing condition along with energy function yields the following generalized eigenvalue problem:

(4)

w y

2

T + Nxy

) dxdy

By substituting the Expressions (7) and (8) into the expressions of kinetic energy and potential energy (1) and (6), meanwhile neglecting the nonlinear components in U , the Rayleigh-Ritz quotient to the corresponding factor is [31]

The potential energy resulting from the thermal loads can be obtained as

U

w x

x a

m

m=1 n=1

¯ k is the transformed reduced stiffness matrix of the kth lamina. where Q ¯ k can be found in Appendix I. The expression of Q Integrating the thermal stresses over the panel thickness H, the internal loads due to the thermal strain can be written as H /2

amn

T ),

xy

k

un

W (x , y ) =

(Tmax T

0

um

(2)

x

t r k yb w 2 + k yb

2

TD

where W (x , y ) can be constructed by Gram-Schmidt orthogonal polynomials, i.e.,

where x , y and xy donate respectively the thermal expansion coefficients of the material constituting the panel along the x and y directions, and in the x-y plane, whose expressions are given in Appendix I. For the kth composite lamina, according to Hook’s law, the stressstrain relation can be expressed as x

a

2

w y

w x

+

w (x , y , t ) = W (x , y ) sin( t ),

x

T xy

k yt 0 w 2 + k yr0

0

In this study, the solution will be sought in weak form using the Rayleigh-Ritz procedure, by which the frequencies and mode shapes of the panel can be obtained. The Gram-Schmidt orthogonal polynomials are employed to construct the mode function. The transverse displacement of the panel can be expressed as

T T

t r k xa w 2 + k xa

a

A

+ NyT

k xt 0 w 2 + k xr 0

b

T 0

2.1. The Rayleigh-Ritz method

where w = w (x , y, t ) is the transverse displacement of the LCP. Since the panel is assumed to work at elevated temperatures, there will be thermal stress caused by thermal expansion mismatch deformation within the panel. Herein, neglecting the temperature difference between upper and lower surfaces, the thermal strain of the panel at elevated temperatures can be written as x

b

2

w x

NxT

Area

1 2

1 ( 2 Area

+ Us =

T

(5)

a = [a11, a12, a13,

, a1un , a21,

= [a1, a2, a3,

Accordingly, for a purely bending LCP, based on the laminated plate theory and Von Kármán large deflection and classical thin plate theories, the total potential energy can be given by [31] (see Nomenclature for notations and Appendix I for the formulation of U ).

M kl = abh

3

1 0

m id

, al 1 0

(11)

aum un ]T (12)

]T , n jd

,

(13)

Composite Structures 229 (2019) 111391

Y. Wang, et al.

bD11 a3

K kl =

1 0

m i

j

d +

1 n

0

1

+

1

+

b T Nx a

K klT =

Kskl

(

1

(0)

n(

0

kr + x21 a 1

1

(1)

m(

0

a T Ny b

b2

d

d

m id

0

1

d 1

n

0 1 0

,

1 0

n

Substituting Eq. (23) into Eq. (6), the expressions of potential and kinetic energies in the form of the matrix are respectively

jd

n jd

j

) j ( )d

( 0T A 0 + T D ) Area 1 A = 2 [qT K1q + 411 (qT K1n1q) 2 +

d +2

n (0) j

(0)

1

1 0

) j ( )d

) i ( )d + m(

k yt1 n (1) j

) i ( )d ,

) = 1,

k+1 (

2(

)=(

)=(

Bk )

B1 )

k(

)

Ck

k 1(

),

Bk =

a

[

b

Ck =

a

k

k(

)]2 d

/

b

(16)

(k = 2, 3, 4, …)

a

[

2 k( ) d /

1( )

k( b a

[

k

, 2 1 ( )] d .

Wwork =

)=

k(

1

)/

0

[

k(

)]2 d .

F p1 =

(18)

(19)

0

m

0, id = 1,

if m i , if m = i

Since the restraints at each edge are the same, the function taken as the same form of φ.

(20)

w t

2

dxdy =

(27) (28)

1 T q Mq, 2

Ksy0,

Ksy1

p + c0

w t

(29)

wdxdy .

p (x , y, t ) Wdxdy = ab Area

Area

t1

(21)

t0

(T

(30)

c0 WT Wdxdy = ab

1 0

1 0

1

0 1 0

p (x , y , t ) Wd d ,

c 0 W T Wd d .

(31) (32)

p1

+ qT F

(33)

p2 q.

U + W ) dt = 0.

(34)

Substituting Eqs. (26)–(29), (33) into Eq. (34), finally, the governing equation in the form of the matrix is written as

Mq¨ + F p2q + (K1 + FT + Ks) q + N (q) q =

F

p1,

(35)

where FT and N (q) are respectively the thermal stiffness matrix and nonlinear items of the equation, whose expressions are given in Appendix I.

(22) can be

2.3. Acoustic load simulation The acoustic excitation is assumed as a stationary White-Gaussian random pressure with zero mean and uniform magnitude over the panel surface. The relation between the time history of sound pressure p (t ) and the cross-spectral density Gp ( ) can be expressed as

p (t ) = FFT

N i

q,

The governing equation of the panel is established based on Hamilton’s equation

With the results from Section 2.1, the formula deduction process of handling the condition of combined thermal and random acoustic loads will be innovatively presented in this section. The transverse displacement of the panel can be expressed as follows

Wi (x , y ) qi (t ) = WT q,

dd

W WT

Then, the external virtual work Wwork can be rewritten as

2.2. Governing differential equation

w (x , y , t ) =

W WT

Area

Wwork = qT F

It is easy to verify that all polynomials satisfy the geometry boundary conditions if the first item 1 (x ) does. Accordingly, the relationship among polynomials can be obtained as follow 1

b T W WT a Nx + NyT a b

Let

The polynomials resulting from Eqs. (17)–(20) can be normalized as k(

(26)

are modal mass matrix and stiffness mawhere M and , trices, respectively. The expressions are given in Appendix I. The external virtual work can be expressed as

F p2 =

)]2 d

1 0

Ksx1,

where b

0

1 h 2 Area

T=

(17)

1,

A22 T 3 (q Kn1q)2 4

1 T x 1 q (Ks0 + Ksx1 + Ksy0 + Ksy1) q = qT Ksq, 2 2

Us =

In order to satisfy the geometry boundary conditions of the elastically restrained panel, the free boundary condition is here chosen to construct the Gram-Schmidt orthogonal polynomials. As explained in Refs. [31–34], each item of the Gram-Schmidt orthogonal polynomials can be obtained as 1(

1

= qT

T

T + abNxy

k yt0 n (0) j (0)

m(

0

(1)

(

n(

+

+ qT Kn2qqT Kn3q + qT Kn4qqT Kn5q],

m (0) i

0

A12 T 2 (q Kn1q)2 2

+ A66 (qT Kn41q)2

U

1

+a

1 2

U=

(14)

(15)

k xr0 a2

(25)

, qN ]T .

q = [q1, q2 , q3,

n jd

0

.

k yr1 + 2 n (1) j b

) i ( )d

m i

m i 1

k xt1 m (1) i (1)

k yr0

m i

1

) j ( )d +

n(

0

) i ( )d +

1

(1)

) j ( )d +

m (1) i m(

0

n(

0

0

0

+

0

1

26 2

n jd

1

m (0) i (0)

1

16 2

(24)

, WN ]T ,

W = [W1, W2, W3,

n j

jd

n jd

where Wi (x , y ) can be obtained from Eq. (7). The mode function vector W and modal coordinate vector q are

m id

1 0

66

0 0

1 0

) + 2Da ( 2D d )+ b ( 4D d )+ ab

1 1

m id

0

n j

d

m i

0 1

= b k xt 0

1

aD22 b3

+

m id

n

n j

0

1

T Nxy

1 0

0

m id

0

1 0

1

d

m i

0

(

d

m i

0

n jd

0

D12 ab

1

d +

1

d

G( ) =

(23) 4

1[

Gp ( 0

G ( ) ei

( ) ], 1)

0

(others )

,

(36) (37)

Composite Structures 229 (2019) 111391

Y. Wang, et al.

traveling in the air). Assuming the sound pressure is uniformly distributed over the plane vertical to the direction of the sound pressure. Accordingly, the lag time t x, y can be calculated based on the position (x, y). The sound pressure at any position of the panel can be rewritten as

p (x , y , t ) = p (x , y, t¯

tx , y ) = p (x , y , t (x , y )),

(41)

where t¯ is the time of the earliest position subjected to the sound pressure along the x-axis. Substituting Eq. (41) into Eq. (31), Eq. (31) can be rewritten as

F p1 = ab = ab

where ϕ(ω) and i respectively represent the random phase spectrum within 0 to 2 and the imaginary number. Based on the Signal Processing, the relation among the power W, bandwidth f , variance of the sound pressure prms , the cross-spectral density Gp ( ) and the SPL can be presented as follows

W = Gp × f = SPL = 10 log

prms p0

F p1 = ab

(38)

SPL 10

0

(39)

).

10

1 0

1 0 n i=1

p (x , y , t ) Wd d p ( , , t ( , )) Wd d .

(42)

1 0

p ( , , t ( , )) Wd d

i/n (i 1)/ n

p¯ (t

ti )

1 0

Wd d ,

(43)

a (i vn

1) cos( ) and p¯ = p (t ti ) sin( ) which is the where ti = vertical component of the sound pressure acting on the surface of the panel.

where p0 represents the standard sound pressure in the air and is equal to 2 × 10 5Pa. Then, by substituting Eq. (38) into Eq. (39), the relation among the SPL, bandwidth, the cross-spectral density can be obtained as [35]

Gp = 4/ f × 10 (

1

= ab

2

,

1 0

Generally, since the sound pressure p(t) cannot be expressed as a continuous function but a discrete sequence used for the integral of Eq. (42), an approximate method to calculate the integral is here proposed to solve this problem, which is based on the infinitesimal calculus. The panel is evenly divided into n strips along the x-axis, as shown in Fig. 4, and the sound pressure at any position of each area is considered to be the same. Therefore, Eq. (42) can be rewritten as

Fig. 2. Band-limited White-Gaussian acoustic excitation (SPL = 160 dB).

2 prms ,

1 0

3. Numerical results and discussions Several natural frequencies of the LCP involving various boundary conditions are calculated first in this section, in order to check the convergence of the solution. The nonlinear vibration behavior of the LCP with all edges elastically restrained is then investigated with two parameters in the study: temperature rise T , sound pressure level (SPL). The example considered here is an eight layered [0/90]4s panel. The dimensions of the LCP are chosen to be a = 0.4 m, b = 0.15 m, h = 0.0021 m, whose material properties are a1 = 0.6 × 10 6 K 1, G12 = 5.403 GPa, = 1616.5 kg/m3, a2 = 29.1 × 10 6 K 1, E1 = 140.35 GPa, E2 = 9.44 GPa, ν12 = 0.253, ν 21 = 0.253 [32].

(40)

As an example, a typical simulated White-Gaussian acoustic excitation time history and its power spectral density (PSD) at SPL = 160 dB are shown in Fig. 2. Considering the acoustic loads from all directions cannot always act vertically on the surface of the panel, the sound pressure p (x , y , t ) should vary randomly with time and position of the panel. Unlike most researches where the noise was considered to act on the surface of the panel vertically, in this paper, the noise is assumed to travel along the xaxis, which is not always vertical to the surface of the panel. And, only one noise source is considered, as shown in Fig. 3. As a result, there is an angle between the normal direction of the panel and the direction of the sound pressure. In Fig. 3, Sx, y is the lag distance of the position (x, y), which is equal to the product of the lag time Δtx,y, and v (the speed of sound

3.1. Validation of the formulation For the elastic restraint boundary, all the classical homogeneous boundary conditions can be directly obtained by accordingly setting the spring constants equal to an extremely large or small number. Thus, simple support (or clamped) edge can be viewed as a special case when the stiffness for the translational (and rotational) boundary springs

Fig. 3. The sound pressure acting on the surface of the panel.

Fig. 4. The panel divided into n strips along the x-axis. 5

Composite Structures 229 (2019) 111391

Y. Wang, et al.

3.2. Nonlinear response under thermo-acoustic loads

Table 1 Natural frequencies (Hz) of the panel with the increase of the orthogonal polynomial’s order. i*

Order of the orthogonal polynomial

This section presents the thermo-acoustic response of an elastic restraint LCP having temperature rises T = 0 , 50, 100, 150, 200 °C, sound pressure level SPL = 140, 150, 160, 170 dB, and then carries out the PSD analysis. Herein, both of um and un are taken as 12, so that the first 6 order modes (N = 6) of the system can be accurate. It can be seen from Table 2 that the natural frequencies of the elastic restraint system are approaching those of the clamped system when k xr and k yr are more than 106 N·m/rad . This means that the condition of the elastic restraint boundary is equivalent to the clamped boundary when the values of the rotational spring’s stiffness are over 106 N·m/rad . Thus, to simulate the situation when the restraint of the boundary is not strong enough, the translational and rotational springs are considered to be a bit weaker than those of simple support and clamped boundaries. The stiffness coefficients k xr and k yr for the rotational boundary springs are taken as 10 4 N·m/rad here. Similarly, the stiffness coefficients for the translational boundary springs are set to be 107 N·m/rad . The first 4 order natural frequencies and mode shapes of the elastic restraint panel are shown in Fig. 5. It can be observed from these mode shapes that the large amplitude is located in the middle point of the panel for the first mode and the 3rd mode. Since the transverse displacement of the panel is dominated the first mode in general, the middle point is selected as a sample point to study the displacement response of the panel thereafter. To investigate the influence of the acoustic loads, we assume that the panel is subjected to the acoustic loads alone and the responses of the system for different SPLs are given for illustration. Let W = w/ h , the responses of the middle point and their PSD for different SPLs are shown in Fig. 6. It can be observed from Fig. 6 that the magnitude of the maximum amplitude is over the thickness of the panel when the SPL is 150 dB. Once the maximum amplitude reaches this level, the complicated nonlinear vibration behaviors will occur. The broadening of bandwidth and shifting of the peaks in the PSD plots are further observed with the increase of SPL. The resonant peaks corresponding to the first and third natural frequencies are shifted clearly when the SPL increase to 170 dB. In addition, observed from the figure for PSD when the SPL = 170 dB, there are contributions in the frequency range over 1500 Hz. Certainly, for a linear system, since the external loads do not contain the signal of such a range, there should not be responses. But, for a nonlinear system, the nonlinear components may gradually play leading roles in the system with the increase of maximum amplitude, resulting in many phenomena such as frequency shifting, internal resonance, and hyper harmonic resonance, etc. For a panel under thermo-acoustic loads, the response indicates three distinct types of motion when ΔT is over the critical buckling temperature: (1) small-deflection random vibration about one of the two thermally buckled equilibrium positions, (2) snap-through or oilcanning phenomenon between the two buckled positions and (3) largeamplitude nonlinear random vibration covering both thermally buckled positions [5]. Herein, to determine the relationship among SPL, ΔT and the three type motions, several examples of different SPL and ΔT are considered and the responses are shown in Fig. 7. The critical buckling temperature of the panel is 138.6 °C, written as ΔTcr. For clarifying the three type motions better, the responses of a 1/4 point (x = 0.25a, y = 0.25b) instead of the middle point are taken as the example, as shown in Fig. 7. This is because the value of the mode function W1 (x , y ) at the middle point is zero when the thermal buckling occurs (Fig. 8 when SPL = 160 dB), and the response of the corresponding mode coordinate q1 dominates the three type motions (as shown in Fig. 9 when SPL = 160 dB, T = 200 °C ) while the other mode coordinates (e.g. q2 q 4 , Fig. 9) have no contribution to those motions

FEM result

6

8

10

12

14

1 2

613.6 678.8

611.7 675.8

611.7 675.8

611.7 675.8

611.7 675.8

609.9 673.6

C-C

3 4 5 6 1 2

1692.0 1729.3 35262.2 35297.7 275.2 342.3

852.7 1187.9 1659.4 1698.2 275.2 341.8

835.5 1114.8 1658.9 1697.6 275.2 341.8

834.8 1110.5 1517.7 1658.9 275.2 341.8

834.8 1110.4 1500.2 1658.9 275.2 341.8

831.8 1105.6 1504.7 1664.3 275.0 341.3

S-S

3 4 5 6

674.1 1068.9 1105.0 1279.3

517.4 818.7 1064.5 1100.6

512.0 792.4 1064.5 1100.6

511.9 791.6 1064.5 1100.6

511.9 791.6 1064.5 1100.6

510.9 789.9 1065.2 1101.4

* i is the order number of the natural frequency.

become infinitely large, which is actually represented by a very large number, such as 10 8 to 1010. In order to check the convergence of the solution, the convergence analysis for the influence of orthogonal polynomial’s order on the natural frequencies is performed, in which the simple support (S-S) and clamped (C-C) boundary conditions are selected as examples. Taking the natural frequency (listed in Table 1) obtained from the ANSYS simulation software as a reference value, the process of obtaining the natural frequencies is validated by comparing the frequency results. It is noteworthy that, as reported in the literature (e.g. in Refs. [30,36]), there is a significant shortcoming in the thin plate theory based on the FEM, which is due to the shear locking phenomenon caused by the inclusion of both the bending and the shear stiffness in a unique rotational degree of freedom [36]. The locking of the finite element solution will result in inaccurate numerical results. Thus, in order to alleviate this problem and make sure the obtained solution from the ANSYS software is reliable, the shell element (Shell 181) is used to establish the FEM model of the LCP, where the corresponding approach is also adopted to increase the reliability of the result. It is obvious with the increase in the order of the orthogonal polynomials listed in Table 1 that the natural frequency of the LCP is gradually convergent to the results of the FEM. Additionally, the natural frequencies of the panel involving different stiffness for the rotational boundary springs are calculated when the stiffness coefficients k xt and k yt for the translational boundary springs are taken as 1010 N/m , which are compared with FEM results of S-S and C-C boundary conditions, as shown in Table 2. It can be observed from Table 2 that the natural frequencies obtained by the Rayleigh-Ritz method match very well with those from FEM when k xr and k yr are taken as 0, 1010 N·m/rad or larger. This shows that the proposed method is valid. Table 2 Natural frequencies (Hz) of the panel under different k xr , k yr . i*

1 2 3 4

FEM for S-S

275.0 341.3 510.9 789.9

FEM for C-C

Rayleigh-Ritz [k xr = k yr (N·m/rad) ] 0

102

10 4

106

108

1010

275.2 341.8 511.9 791.6

290.9 356.5 524.0 801.4

549.3 612.2 767.8 1035

610.9 675.0 834.0 1110

611.7 675.8 834.8 1111

611.7 675.8 834.8 1111

609.9 673.6 831.8 1106

* i is the order number of the natural frequency. k xr , k yr are the stiffness

coefficients for the rotational boundary springs.

6

Composite Structures 229 (2019) 111391

Y. Wang, et al.

Fig. 5. The natural modes of the panel.

Fig. 6. Time histories and PSD of the panel with different SPLs.

7

Composite Structures 229 (2019) 111391

Y. Wang, et al.

Fig. 7. Time histories for the panel with different SPLs and at different temperature rises.

Fig. 8. The mode of the panel corresponding to the mode function W1(x,y) and the time history of the middle point with SPL = 160 dB for T = 200 °C respectively. 8

T = 150 °C and

Composite Structures 229 (2019) 111391

Y. Wang, et al.

Fig. 9. Time histories of different coordinates (q1 q

after buckling. It means that if the middle point was taken as the example, there would be no three type motions seen in the responses of the panel, as shown in Fig. 8. When ΔT is less than the critical buckling temperature, the panel exhibits the basically small deflection random vibration as is shown in Fig. 7. With the increase of ΔT, the panel shows a gradually weaker stiffness as the root mean square (R.M.S.) of the response in Fig. 7 increases under different SPL acoustic loads. Additionally, it is worth noting that the increase of R.M.S. with the rise of ΔT will enlarge as the acoustic loads become stronger. When the buckling of the panel occurs, the responses under different conditions indeed show the three types of motion. It is found from Fig. 7a–c that the thermal post-buckling deflections dominated the responses, and the nonlinear stiffness that comes from thermal deflection is seen to counteract the acoustic pressure, leading to a small-amplitude random vibration about one of the

4)

when SPL = 160 dB and T = 200 °C .

equilibrium buckling positions at T = 150 °C or T = 200 °C and SPL = 140 dB or SPL = 150 dB, as shown in Fig. 7a, b. With the increase of SPL, the response will show the snap-through phenomenon at the same condition of ΔT, such as when T = 200 °C, SPL = 160 dB. If the acoustic load continues becoming strong, the snap-through motion will be more frequent as shown in Fig. 7c T = 150 °C and Fig. 7d T = 200 °C, or the thermal buckling will fail to counteract the acoustic pressure, leading to a large-amplitude nonlinear random vibration covering both thermally buckled positions, as shown in Fig. 7d T = 150 °C. 3.3. Analysis of stress/strain and fatigue life prediction In the engineering application, both the fatigue life prediction and strength evaluation are problems of the main concerns. Before

Fig. 10. The 1st order stress/strain PSD distributions of the surface of the panel (SPL = 160 dB and T = 0 ). 9

Composite Structures 229 (2019) 111391

Y. Wang, et al.

evaluating the strength and predicting the fatigue life of a structure, the most dangerous location of stress/strain in the structure should be found first. Since the distribution of stress/strain is often not in correspondence with that of displacement, the most dangerous location of stress/strain cannot be predicted directly from the mode shapes, but it can be predicted by the PSD distribution of stress/strain. Generally, the 1st order resonance contains enormous energy. Accordingly, the energy distribution of the first order resonance can provide good reference resources for the dangerous location prediction. Considering that the acoustic load may lead to the failure of the panel when SPL = 170 dB, therefore, taking SPL = 160 dB and T = 0 as the example, the frequency is set to be 515 Hz, where the first order resonance of the panel occurs. The 1st order PSD distribution of stress/ strain can be calculated via the constitutive equations [31] after obtaining the displacement time histories, as shown in Fig. 10. The strain PSD distributions of the surface along x-axis and y-axis as shown in Fig. 10a, b shows that the most strains are located at the 1/2 points of the panel’s edge, while the figure (Fig. 10c) for that of shear (xy-direction) strain shows that the most shear strains lie at the 1/4 points of the panel. It can be observed from Fig. 10d for the Von Mises stress PSD distribution that the most stress is located at the middle point. Comparing Fig. 10a with b, it is also obvious that the stress along x-axis dominates the distribution of the Von Mises stress. Although the dangerous locations of stress can be obtained from the PSD distribution of the Von Mises stress, however, the Von Mises stress is not a valid reference resource for the fatigue life prediction and even strength evaluation, since the LCP is not isotropic. It means that the fatigue strength, S-N curves and ultimate strength of the panel vary from different directions, which are extremely time-consuming and costly to obtain by experiments. At the meantime, the direction of the Von Mises stress also varies with time. Thus, it is quite hard to use Von Mises stress to predict the fatigue life and evaluate the strength of the LCP. Additionally, the fatigue failure mechanism of the laminated composite is very complicated and has not been clarified even though extensive research work on the fatigue of LCPs has been carried out

since the 1960s. Based on the experimental investigation, Reifsnider [37] concluded that fatigue damage evolution is nonlinear in composite materials. During the initial period of fatigue life, many non-interactive cracks occur in the matrix. When the matrix crack density reaches saturation, the fiber failure, interfacial deboning and delamination occur in the composites. The stress of the cross section will be redistributed and increased, leading to the rapid development of damage and causing ‘‘sudden death” of the material in the end period of fatigue life [38]. Based on plenty of experimental investigations, many damage models [37–42], which have been, respectively, defined by strength degradation, stiffness degradation and energy dissipation of composites, have been employed to macroscopically describe the damage development of materials in the recent decades. But, the sonic fatigue mechanism of the panel is different since the deformation of the panel’s cross-section is not uniform. This means the stress distribution of the cross-section under the thermo-acoustic load is quite different from that under the axial tension or compression force, and accordingly, the damage caused by the time-varying cycling load varies with the location of the cross section even if the panel consists of the same material. When it comes to the LCP, the mechanism is more complicated since the fatigue properties change with the ply angle and materials of every lamina. Therefore, the existing theories are less employed directly to predict the fatigue life of the LCP under the thermo-acoustic loads. But, considering that 1) the failure of a LCP often occurs layer by layer, 2) the failure initiates by fiber breakage at their weakest cross sections and phases [43], the rough fatigue life prediction of the LCP can be carried out by determining when and where the initial fatigue failure occurs in a single lamina. As known from the constitutive equations [31], the most stresses along x-axis and y-axis often lie on the surface (ply-0° ) of the panel in the simulation. But, from Fig. 10b, it can be observed that the most strain of the surface along y-axis is much larger than that along x-axis. This means that 1) the fatigue failure along y-axis may occur first, and 2) the strain along y-axis in the second lamina (ply-90° ) from the surface may be also large, resulting in the strong stress along its fiber. Thus,

Fig. 11. The 1st order stress PSD distributions of the surface and the second laminate of the panel (SPL = 160 dB and T = 0 ). 10

Composite Structures 229 (2019) 111391

Y. Wang, et al.

Fig. 12. The dimensionless fatigue damage distributions of the surface and the second laminate of the panel (SPL = 160 dB and T = 0 ).

there will be both equally important things that have to be done in fatigue life prediction. That is to determine in which direction and lamina the initial failure will occur respectively. As shown in Fig. 11, taking SPL = 160 dB and T = 0 as the example, the most stresses along x-axis and y-axis on the surface and the second lamina from the surface can be found in the 1st order stress PSD distributions.

Comparing Fig. 11a with Fig. 11b, it is obvious that the most stress along the fiber is located at the 1/2 points of edge y = 0 or y = b on the second lamina (ply-90° ) from the surface instead of on the surface. Additionally, the most stress along the direction vertical to the fiber still lies at the 1/2 points of edge y = 0 or y = b on the surface. Referring the experiment data from [44], the fatigue life of a single

Fig. 13. The 2nd order stress PSD distributions of the surface and the second laminate of the panel (SPL = 160 dB and T = 0 ). 11

Composite Structures 229 (2019) 111391

Y. Wang, et al.

lamina can be estimated by the usage of both rain flow cycle counting method and conventional S-N curves, together with the Miner’s theory when the location of the most possible failure is found [21]. In order to validate whether the 1st order stress PSD distribution shown in Fig. 11 can provide a good reference resource for fatigue life prediction, the dimensionless fatigue damage distributions of the corresponding laminas during 5 s under the same condition are shown in Fig. 12. It can be seen from the Fig. 12 that the initial fatifue failure is estimated to occur around the middle points of edge y = 0 or y = b on the surface, which is caused by the stress along y-axis, i.e. the direction vertical to the fiber. Comparing Fig. 11 with Fig. 12, there are some differences, although the results of the damage distribution match those of the 1st order stress PSD distribution well. Actually, these differences are mainly caused by the higher order modes (e.g. The 2nd order mode, whose stress PSD distributions are shown in Fig. 13) which often have a small contribution to the dynamic stress of the panel at most time. In this simulation, the 2nd order mode obviously has a noticeable effect on the dynamic stress of the panel. Thus, with due consideration for the higher order modes, the dangerous locations for fatigue life prediction can be found more accurately. However, the 1st order stress PSD distribution is proved to be valid for providing a good reference resource for fatigue life prediction. This may help to determine the most dangerous areas for fatigue life prediction, even strength evaluation preliminarily, which can narrow the search scope for localization and improve efficiency.

combined acoustic and thermal loads have been carried out by numerical simulation, where the angled sound pressure is innovatively considered. The results of the nonlinear random response under different combinations of sound pressure level and temperature environments demonstrate that all three types of panel motion: (1) small-deflection random vibration about one of the two thermally buckled equilibrium positions, (2) snap-through or oil-canning phenomenon between the two buckled positions and (3) large-amplitude nonlinear random vibration covering both thermally buckled positions—can be accurately predicted. It is also indicated from the R.M.S. of the response that the elevated temperature will contribute to the large amplitude of vibration and this contribution will enlarge as the acoustic load becomes stronger. The stress/strain PSD distributions are carried out to determine the most probable areas of failure in the panel. In Addition to the 1st order, the higher order stress/strain PSD distributions are found to be noticeable for determining the most dangerous locations for fatigue life prediction, and the results show that the most dangerous locations mainly lie at 1/2 points of the edges and their nearby on the surface or the second lamina from the surface. But, for different elastic boundary conditions, ply angles, structure sizes and particularly composite materials, etc., it is worth noting that the dangerous locations may change due to the change of the dynamical characteristics and the fatigue properties. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

4. Conclusions The dynamic model of the LCP elastically restrained along the edges has been established based on Von Kármán large deflection and classical thin plate theories, together with the Hamilton principle. The Rayleigh-Ritz method using orthogonal polynomials has been described for the time-domain modal formulation of the panel. The modal analysis and nonlinear random responses of a sample LCP subjected to

Acknowledgment This work is supported by the National Natural Science Foundation of China under Grant Nos: 11472089 and 11732005.

Appendix I The thermal expansion coefficients of the material constituting the panel along the x and y directions, and in x-y plane can be calculated from

= =

x y xy

2 1cos 2 1sin

= 2(

+ + 2)

1

2 2 sin 2 2 cos

, , (44)

cos sin .

The transformed reduced stiffness matrix of the k-th lamina Q can be obtained by

¯k

¯k

(45)

Q = T k 1Q [T k 1]T , cos2 k sin2 k sin k cos

Tk =

k

cos2 k cos2 k sin k cos

k

2 sin k cos k 2 sin k cos k , cos2 k sin2 k

(46)

Q11 Q12 0 Q = Q21 Q22 0 , 0 0 Q66

(47) th

where Tk and k are the transformation matrix and the ply angle of k The entries of Q are given by

Q11 = E1/(1

v12 v21), Q22 = E2 /(1

v12 v21), Q12 = v12 E2 /(1

lamina respectively, and Q is the panel stiffness matrix.

v12 v21) = v21 E1/(1

v12 v21), Q66 = G12.

(48)

where E1, E2 and G12 donate respectively the elastic and shear modulus, and v12 , v21 donate the Poisson’s ratios of the laminated composite panel. For the LCP, the membrane strain 0 is determined as

0

=

0 x 0 y 0 xy

=

1 2

w x

2

,

w y

2

,

w w x y

T

. (49)

12

Composite Structures 229 (2019) 111391

Y. Wang, et al.

The bending strain κ is given as x

=

2w

=

y xy

2w

,

x2

y2

T

2w x y

,

.

(50)

Based on Eq. (49) and Eq. (50), the total strain of the LCP can be expressed as

=

+z

x

+z

y

0 xy

+z

xy

=

y xy

T

0 x 0 y

x

=

+z .

0

(51)

It is assumed that the angle of the LCP is symmetrically laid relative to the middle plane. Thus, the potential energy of the panel is expressed as

U =

1 ( T0 A 2 Area

+

0

TD

) dxdy.

(52)

The modal stiffness matrices of Eq. (26) are calculated as 1

K1 = ab

K1n1 = K3n1 =

1

0

0

2W

D11 a4

1 0

1 b W WT d 0 a

1 0

1 a W WT 0

b

2W T

2

2

+

1 A16 W WT d 0 a

Kn 4 = a

1 0

1 A26 W WT d 0 b

+

2

W WT

1 0

1 0

W WT

d , Kn3 =

1 0

1 0

W WT

d , Kn5 =

1 0

1 0

W WT

d d , Kn41 =

1 0

2W T

2

1 0

1 0

d , Kn21 =

Kn2 = b

2W

2D12 a2b2

4D16 a3b

2W

2WT

2

+

D22 b4

2W 2

2W T 2

+

4D26 b3a

2W 2

2W T

+

4D66 a2b2

2W

2W T

dd ,

(53)

dd , dd ,

(54)

dd , dd ,

(55)

The modal stiffness matrices of Eq. (28) are calculated as

Ksx0 = b

1 0

k xt 0 (W)

T = 0 (W) = 0

+

Ksx1 = b

1 0

k xt1 (W)

T = 1 (W) = 1

+

Ksy0 = a

1 0

k yt 0 (W)

T = 0 (W) = 0

+

Ksy1 = a

1 0

k yt1 (W)

T = 1 (W) = 1

+

( ( ( (

k xr 0 a2

k xr1 a2

W

k yr0 b2 k yr1 b2

) ) ) )

W

W

W

=0

=1

) ) ) )

W

=0

=1

W

( ( ( (

T =0

T

=1

W

W

T =0

T =1

d , d , d , d .

(56)

The modal mass matrix of Eq. (29) is calculated as

M=

1 0

1 0

abph WWT d d .

(57)

The thermal stiffness matrix and nonlinear items of Eq. (35) are expressed respectively as

FT =

1 0

N (q) =

1 0

b T W WT a Nx + NyT a b

W WT

T + abNxy

W WT

dd ,

(58)

A11 1 A A [Kn1 + (K1n1)T ] qqT K1n1 + 12 [Kn2 1 + (K2n1)T ] qqT Kn21 + 22 [Kn3 1 + (K3n1)T ] qqT Kn3 1 + A66 [Kn41 + (Kn41)T ] qqT Kn4 1 4 2 4 +

1 (Kn2 + KTn2) qqT Kn3 + (Kn3 + KTn3) qqT Kn2 + (Kn4 + KTn4) qqT Kn5 + (Kn5 + KTn5) qqT Kn4 . 2

(59)

[6] Lee J, Wentz K, Clay C, Anselmo E, Crumbacher R, Vaicaitis R. Prediction of statistical dynamics of thermally buckled composite panels. In: 39th AIAA/ASME/ ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit, 1998. p. 1975. [7] McEwan M, Wright JR, Cooper JE, Leung AYT. A combined modal/finite element analysis technique for the dynamic response of a non-linear beam to harmonic excitation. J Sound Vib 2001;243:601–24. [8] Ko WI. Thermal buckling analysis of rectangular panels subjected to humped temperature profile heating. Report for Dryden flight research center report no TP2004-212041. Washington: NASA; 2004. [9] Dhainaut J-M, Guo X, Mei C, Spottswood SM, Wolfe HF. Nonlinear random response of panels in an elevated thermal-acoustic environment. J Aircraft 2003;40:683–91. [10] Locke J, Mei C. Finite element, large-deflection random response of thermally buckled beams. AIAA J. 1990;28:2125–31. [11] Mei C, Chen RR. Finite element nonlinear random response of composite panels of

References [1] Ibrahim HH, Yoo HH, Tawfik M, Lee K-S. Thermo-acoustic random response of temperature-dependent functionally graded material panels. Comput Mech 2010;46:377–86. [2] Li X, Yu K. Vibration and acoustic responses of composite and sandwich panels under thermal environment. Compos Struct 2015;131:1040–9. [3] Ibrahim H, Tawfik M, Negm H. Random response of shape memory alloy hybrid composite plates subject to thermo-acoustic loads. J Aircraft 2008. [4] Li X, Yu K, Zhao R, Han J, Song H. Sound transmission loss of composite and sandwich panels in thermal environment. Compos B Eng 2018;133:1–14. [5] Liu L, Li J-M, Kardomateas GA. Nonlinear vibration of a composite plate to harmonic excitation with initial geometric imperfection in thermal environments. Compos Struct 2019;209:401–23.

13

Composite Structures 229 (2019) 111391

Y. Wang, et al.

[12] [13] [14] [15] [16]

[17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

arbitrary shape to acoustic and thermal loads. Norfolk, VA: Old Dominion Univ.; 1997. Locke J. Finite element large deflection random response of thermally buckled plates. J Sound Vib 1993;160:301–12. Przekop A, Guo X, Azzouz S, Mei C. Nonlinear response and fatigue life of shallow shells to acoustic excitation using finite element. In: 44th AIAA/ASME/ASCE/AHS/ ASC structures, structural dynamics, and materials conference; 2003. p. 1710. Radu A, Kim K, Yang B, Mignolet M. Prediction of the dynamic response and fatigue life of panels subjected to thermo-acoustic loading. In: 45th AIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics & Materials Conference; 2004. p. 1557. Przekop A, Rizzi SA. Dynamic snap-through of thin-walled structures by a reducedorder method. AIAA J 2007;45:2510–9. Shukla A, Gordon R, Hollkamp J. Numerical investigation of the snap-through response of a curved, clamped-clamped plate with thermal and random loading. In: 49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 16th AIAA/ASME/AHS Adaptive Structures Conference, 10th AIAA Non-Deterministic Approaches Conference, 9th AIAA Gossamer Spacecraft Forum, 4th AIAA Multidisciplinary Design Optimization Specialists Conference; 2008. p. 2230. Sha YD, Li JY, Gao ZJ. Dynamic response of pre/post buckled thin-walled structure under thermo-acoustic loading. Applied mechanics and materials. Trans Tech Publ; 2011. p. 536–41. Sha YD, Gao ZJ, Xu F, Li JY. Influence of thermal loading on the dynamic response of thin-walled structure under thermo-acoustic loading. Advanced engineering forum. Trans Tech Publ; 2012. p. 876–81. Sha YD, Gao ZJ, Xu F. Influences of thermal loads on nonlinear response of thinwalled structures in thermo-acoustic environment. Applied mechanics and materials. Trans Tech Publ; 2012. p. 220–6. Sha YD, Wei J, Gao ZJ, Zhong HJ. Nonlinear response with snap-through and fatigue life prediction for panels to thermo-acoustic loadings. J Vib Control 2014;20:679–97. Sha YD, Zhu L, Zhang GZ, Luan XC. Nonlinear response and fatigue life of thin plates subjected to non-uniform temperature distributions and acoustic excitations. Applied mechanics and materials. Trans Tech Publ; 2013. p. 1659–69. Sha Y-d, Zhu L, Jie X-B, Luan X-C, Feng F-F. Nonlinear random response and fatigue life estimation of curved panels to non-uniform temperature field and acoustic loadings. J Vib Control 2016;22:896–911. Istenes DR, Jr., Raymond, Rizzi S, Wolfe H. Experimental nonlinear random vibration results of thermally buckled composite panels. In: 36th Structures, Structural Dynamics and Materials Conference; 1995. p. 1345. Ng CF, Clevenson SA. High-intensity acoustic tests of a thermally stressed plate. J Aircraft 1991;28:275–81. Murphy KD, Virgin L, Rizzi S. Characterizing the dynamic response of a thermally loaded, acoustically excited plate. J Sound Vib 1996;196:635–58. Er GK, Iu VP. Probabilistic solutions of a nonlinear plate excited by gaussian white noise fully correlated in space. Int J Struct Stab Dyn 2017;17:1750097. Kumaraian ML, Rebbagondla J, Mathew TV, Natarajan S. Stochastic vibration

[28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44]

14

analysis of functionally graded plates with material randomness using cell-based smoothed discrete shear gap method. Int J Struct Stab Dyn 2019;19:1950037. Venkatachari A, Natarajan S, Ramajeyathilagam K, Ganapathi M. Assessment of certain higher-order structural models based on global approach for bending analysis of curvilinear composite laminates. Compos Struct 2014;118:548–59. Venkatachari A, Natarajan S, Haboussi M, Ganapathi M. Environmental effects on the free vibration of curvilinear fibre composite laminates with cutouts. Compos B Eng 2016;88:131–8. Venkatachari A, Natarajan S, Ganapathi M, Haboussi M. Mechanical buckling of curvilinear fibre composite laminate with material discontinuities and environmental effects. Compos Struct 2015;131:790–8. Lin H, Cao D, Xu Y. Vibration characteristics and thermoelastic analysis of the supersonic composite laminated panel with arbitrary elastic boundaries. Proc Instit Mech Eng, Part G: J Aerospace Eng 2018. 0954410018781484. Bhat R. Natural frequencies of rectangular plates using characteristic orthogonal polynomials in Rayleigh-Ritz method. J Sound Vib 1985;102:493–9. Liu L, Cao D, Sun S. Vibration analysis for rotating ring-stiffened cylindrical shells with arbitrary boundary conditions. J Vib Acoust 2013;135:061010. Sun S, Cao D, Han Q. Vibration studies of rotating cylindrical shells with arbitrary edges using characteristic orthogonal polynomials in the Rayleigh-Ritz method. Int J Mech Sci 2013;68:180–9. He EML. Nonlinear vibration response analysis and fatigue life prediction of a thinwalled structure under thermal-acoustic loading. J Vibration Shock 2013;32:135–9. Rodrigues J, Natarajan S, Ferreira A, Carrera E, Cinefra M, Bordas S. Analysis of composite plates through cell-based smoothed finite element and 4-noded mixed interpolation of tensorial components techniques. Comput Struct 2014;135:83–7. Reifsnider K, Henneke E, Stinchcomb W, Duke J. Damage mechanics and NDE of composite laminates. Mechanics of composite materials. Elsevier; 1983. p. 399–420. Mao H, Mahadevan S. Fatigue damage modelling of composite materials. Compos Struct 2002;58:405–10. Subramanian S, Reifsnider K, Stinchcomb W. A cumulative damage model to predict the fatigue life of composite laminates including the effect of a fibre-matrix interphase. Int J Fatigue 1995;17:343–51. Epaarachchi JA, Clausen PD. A new cumulative fatigue damage model for glass fibre reinforced plastic composites under step/discrete loading. Compos A Appl Sci Manuf 2005;36:1236–45. Broutman L, Sahu S. A new theory to predict cumulative fatigue damage in fiberglass reinforced plastics. Composite materials: testing and design (second conference). ASTM Int. 1972. Hwang W, Han KS. Fatigue of composite materials—damage model and life prediction. Composite materials: fatigue and fracture, second volume. ASTM International; 1989. Jones RM. Mechanics of composite materials. CRC Press; 2014. Jen M-H, Lee C-H. Strength and life in thermoplastic composite laminates under static and fatigue loads. Part I: Experimental. Int. J. Fatigue 1998;20:605–15.