Aeroelastic behavior of composite panels undergoing progressive damage

Accepted Manuscript Aeroelastic behavior of composite panels undergoing progressive damage Douglas Quintanilha Tsunematsu, Maurício Vicente Donadon PII: DOI: Reference:

S0263-8223(18)31450-8 https://doi.org/10.1016/j.compstruct.2018.11.065 COST 10430

To appear in:

Composite Structures

Please cite this article as: Tsunematsu, D.Q., Donadon, M.V., Aeroelastic behavior of composite panels undergoing progressive damage, Composite Structures (2018), doi: https://doi.org/10.1016/j.compstruct.2018.11.065

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.

Aeroelastic behavior of composite panels undergoing progressive damage Douglas Quintanilha Tsunematsu∗, Maur´ıcio Vicente Donadon Divis˜ ao de Engenharia Aeron´ autica, Instituto Tecnol´ ogico de Aeron´ autica, S˜ ao Jos´ e dos Campos, SP, Brazil

Abstract A finite element model for predicting the nonlinear aeroelastic behavior of composite panels undergoing intralaminar and translaminar progressive damage in supersonic flow is presented. The classical plate theory in conjunction with the von K´ arm´ an nonlinear strains is used for structural modeling, and the linear piston theory is used to model the aerodynamic loads. Progressive damage is modeled by a smeared cracking formulation in which stress-based, continuum damage mechanics and fracture mechanics approaches are combined. No modal reduction is performed and an iterative form of the Newmark method is used for the numerical direct integration in time of the nonlinear equations. Simulations considering different lay-ups are conducted, in which the influence of progressive damage on the aeroelastic behavior of the panels is investigated, and damage extent and failure mechanisms are assessed. The results obtained in the analyses consist in important insights concerning the flutter-induced damage in composite panels. Keywords: Aeroelasticity, Panel flutter, Composites, Progressive damage, Finite element method

1. Introduction Panels that constitute the external structure of aerospace vehicles which operate at supersonic speeds are prone to present an aeroelastic instability called panel flutter. This instability is characterized by a self-excited oscillation of the panel after the flow velocity has reached a critical level. 5

Many authors have studied the aeroelastic behavior of panels over the past decades, and notable surveys conducted by Dowell [1] and Mei et al. [2] give an outstanding overview of the subject. It is widely known that the occurrence of flutter may damage the panel, leading to catastrophic or fatigue failure. However, this matter has been scarcely investigated until now. Most of the studies found in the literature are focused on evaluating the influence of some kind of damage present in the panel on the critical

10

flutter speed. Chen and Lin [3] studied cracked isotropic panels in supersonic flow, and their results showed that the existence of cracks decreases the speed at which flutter occurs. Similarly, cracked composite panels were considered in the works of Lin et al. [4], Pidaparti [5] and Wang et al. [6]. Their results revealed ∗ Corresponding

author Email address: [email protected] (Douglas Quintanilha Tsunematsu)

Preprint submitted to Composite Structures

October 16, 2018

that the presence of a crack in the composite panel usually also reduces the flutter boundary. However, some exceptions were found within certain range of fiber orientation angles. Natarajan et al. [7] conducted 15

an analysis for functionally graded material panels with cracks. Recently, Abdullah et al. [8] studied the aeroelastic behavior of cracked composite plates using the doublet lattice method for the aerodynamic loads modeling. Shiau [9] investigated the effects of delaminations on the flutter boundaries of composite panels. The results of his study showed that in most cases the existence of a delamination also decreases the critical flutter speed. Yazdi and Rezaeepazhand [10] studied the flutter characteristics of composite plates with

20

delaminations using structural similitude. In a recent work by Fazilati [11], delaminated composite panels with curvilinear fibers were considered. Pidaparti [12] studied damaged composite plates in supersonic flow using two different stiffness degradation models. The results once more showed that the stiffness degradation reduces the flutter boundary. Kim et al. [13] presented a study of bimodular flutter oscillations of microcrack damaged composite panels.

25

Nonlinear analyses which incorporate geometric nonlinearities in order to determine the amplitudes of damaged fluttering panels were also conducted. AsadiGorgi et al. [14] investigated the effects of all-over part-through cracks on the aeroelastic behavior of isotropic panels. The results showed that the amplitudes of the flutter motion can increase or decrease depending on the location of the cracks. Dowell [15] and Xue and Mei [16] provided some examples of fatigue life estimation for fluttering isotropic

30

plates and discussed their implications for design. Strganac and Kim [17] studied the aeroelastic behavior of composite panels considering the natural progression of matrix microcracks. The growth of the microcrack density was modeled using an internal state variable which is a function of the number of cycles. In all the previously cited works, whether damage was already present in the panel since the beginning of the analysis and does not change with time, or gradually progresses with the number of cycles in a

35

fatigue context. To the best of the authors’ knowledge, the damage induced in a composite panel by the high stresses that arise throughout the flutter motion, and its progression and influence on the aeroelastic behavior of the structure have never been investigated. A better understanding of the effects of progressive damage on the aeroelastic behavior of composite panels is of paramount importance for conducting a damage tolerance design. In order to contribute to this matter, this work presents a finite element model for

40

predicting the nonlinear aeroelastic behavior of composite panels undergoing intralaminar and translaminar progressive damage in supersonic flow. The classical plate theory is used for structural modeling and geometric nonlinearities are considered through the use of the von K´arm´an strains. The quasi-steady firstorder piston theory aerodynamics is used to model the aerodynamic loads. Progressive damage is modeled by a smeared cracking formulation in which stress-based, continuum damage mechanics and fracture mechanics

45

approaches are combined. No modal reduction is performed and an iterative form of the Newmark method is used for the numerical direct integration in time of the nonlinear equations. Simulations considering different lay-ups are conducted, in which the influence of progressive damage on the aeroelastic behavior of 2

the panels is investigated, and damage extent and failure mechanisms are assessed.

2. Finite element formulation 50

Consider a flat rectangular laminated panel of length a, width b and thickness h, subjected to a supersonic flow with velocity V , located in a Cartesian coordinate system defined by the axes x, y, and z, as illustrated in Fig. 1. The length a is aligned with the x direction and the width b is aligned with the y direction. The airflow above the panel is in the positive x direction and the effects of the cavity under the panel are neglected. The principal material coordinate system of a ply that constitutes the laminate is defined by

55

the axes x1 and x2 , which are parallel and normal to the fiber direction, respectively. The fiber orientation angle θ, which is the angle between the axes x and x1 , is measured counterclockwise from the x direction. Due to a constraint of minimum weight, the panels that constitute the aerospace structures are usually thin. In this case, the effects of transverse shear deformations can be neglected, and the classical plate theory can be employed for the structural modeling. Additionally, in order to determine the amplitudes of

60

the flutter motion, geometric nonlinearities must be incorporated in the aeroelastic model. For this purpose, the von K´ arm´ an nonlinear strains can be used. With this, the von K´arm´an nonlinear strain-displacement relations expressed in terms of the classical plate theory displacement field are given by [18]       
 ∂u 1 ∂w 2  ∂2w         −        2 x  ∂x    ∂x   2 ∂x                   
2 2 m ∂v 1 ∂w ∂ w = + + z = {m y − L } + {N L } + z {κ} 2 ∂y 2 ∂y ∂y                                 ∂u ∂v     ∂w ∂w    ∂2w  γxy −2 ∂x∂y ∂y + ∂x ∂x ∂y

(1)

where x , y and γxy represent the state of strain in the laminate coordinate system and u, v and w are the midplane displacements in the x, y and z directions, respectively. 65

The laminate constitutive relations are given by [18]    {N }  {M }



[A]

= [B]

     {m } {m } L NL   = [C] {} +  {κ}   {0}  [D] [B]

(2)

where {N } and {M } are the in-plane force and moment resultants, and the matrices [A], [B] and [D], that form the laminate constitutive matrix [C], represent the extensional, bending-extensional coupling and bending stiffnesses, respectively. These quantities are expressed in terms of the lamina constitutive matrix in the principal material coordinate system [Q] and the fiber orientation angle θ. The terms of the matrix 70

[Q] are calculated as a function of E1 , E2 , G12 and ν12 , which are the lamina elastic moduli, lamina shear modulus and Poisson’s ratio, respectively. Since the flow above the panel is in the supersonic regime, the aerodynamic load can be modeled by the

3

quasi-steady first-order piston theory [19]. In this case, the aerodynamic pressure is given by [20] 2q pa = − β



∂w M 2 − 2 1 ∂w + 2 ∂x M − 1 V ∂t



0 D0 ∂w D11 ∂w = − λ 11 + g a a3 ∂x ω0 a4 ∂t



 (3)

0 where λ is the nondimensional aerodynamic pressure, ga the nondimensional aerodynamic damping, D11 75

the reference stiffness, ω0 the reference frequency, q the dynamic pressure, M the Mach number, ρm the √ material density and β = M 2 − 1. The Hamilton’s Principle states that Z

t2

(δU − δT − δWN C ) dt = 0

(4)

t1

where δU is the virtual strain energy, δT the virtual kinetic energy, δWN C the virtual work done by nonconservative forces, t1 the initial time and t2 the final time. For a finite element, these quantities 80

are expressed as Z

T

{δ}

δU = A

   {N }  {M }

 T  ρ h 0 0    δu     m       ρm h 0  δv    0    Z      0 ρm h  0 δT = − δw   A     δw    0 0 0  ,x            δw  ,y 0 0 0 Z

0

0 ρm h3 12 0

A



Z δw A

T

{δ} [C] {} dA

(5)

A

   u    ¨        0    v¨     Z    T 0  w {δum } [Im ] {¨ um } dA   ¨  dA = − A     w 0   ¨     ,x    3   ρm h  w ¨,y  12 0

0

δwpa dA = −λ

δWN C =

Z dA =

  0  Z 0 D11 D11 w,x dA − ga δw w˙ dA 3 a ω0 a4 A

(6)

(7)

where A is the element area and comma and dot represent spatial and time derivative, respectively. Since a considerable number of elements are usually used in a progressive damage analysis, an element with as few degrees of freedom as possible was searched for the finite element discretization in order to reduce the computational cost. For this purpose, the Bogner-Fox-Schmit conforming quadrangular element 85

[21] shown in Fig. 2 is employed, since it is one of the most simple conforming elements for finite element models that employ the classical plate theory. It has 4 nodes and 6 degrees of freedom per node. The 24

4

nodal displacements can be represented as   e    {w1 }        {we }  2 e {w } = {we }    3       {we }  4

(8)

with   i    u          i     v         wi   e {wi } = i     w,x        i      w   ,y      w i  

(i = 1, 2, 3, 4)

(9)

,xy

where {wie } represents the nodal displacements of a single node. 90

The BFS element uses bicubic Hermite and bilinear Lagrange polynomials to interpolate the bending and membrane degrees of freedom, respectively [22]. Thus, the following approximations can be written

{} = [B] {we }

(10)

  ¯ {δwe } {δ} = B

(11)

{um } = [Hm ] {we }

(12)

{δum } = [Hm ] {δwe }

(13)

w = [Ha ] {we }

(14)

δw = [Ha ] {δwe }

(15)

w,x = [Ba ] {we }

(16)

5

  ¯ , [Hm ], [Ha ] and [Ba ] are the interpolation matrices. where [B], B Substituting Eqs. (10)-(16) in Eqs. (5)-(7), and the resultant equations in Eq. (4), we can write Z

t2

T

{δwe }

 e  e   e  e e [m ] w ¨ + ga [cea ] w˙ e + ([kN dt = 0 L ] + λ [ka ]) w

(17)

t1

with

e [kN L] =

Z

 T ¯ [C] [B] dA B

(18)

A

[me ] =

Z

T

[H] [Im ] [H] dA

(19)

A

0 D11 a3

Z

D0 = 114 ω0 a

Z

[kae ] =

[cea ] 95

T

(20)

T

(21)

[Ha ] [Ba ] dA A

[Ha ] [Ha ] dA A

e e e where [kN L ] is the element nonlinear stiffness matrix, [m ] the element mass matrix, [ka ] the element

aerodynamic stiffness matrix and [cea ] the element aerodynamic damping matrix. Due to the arbitrariness of {δwe }, the term multiplied by it in Eq. (17) must vanish. Hence, the equations of motion for a single element can be written as

 e   e e e [me ] w ¨ + ga [cea ] w˙ e + ([kN = {0} L ] + λ [ka ]) w

(22)

Considering the boundary conditions, the system equations of motion are assembled from the element 100

equations of motion. Hence, we can write    ¨ + ga [Ca ] W ˙ + ([KN L ] + λ [Ka ]) W = {0} [M ] W

(23)

where [KN L ] is the global nonlinear stiffness matrix, [M ] the global mass matrix, [Ka ] the global aerodynamic stiffness matrix, [Ca ] the global aerodynamic damping matrix and {W } the global nodal displacement vector.

3. Progressive damage model The progressive damage model proposed in this work is based on the efforts of Donadon et al. [23, 105

24, 25, 26]. The formulation combines stress-based, continuum damage mechanics and fracture mechanics approaches. Initially, the damage onset is detected by a set of stress-based criteria for the different failure modes of a unidirectional ply under plane stress state. After that, irreversible damage variables which are 6

based on a smeared cracking formulation are calculated. These damage variables gradually degrade the stiffness in the lamina level of each element of the structural model. The intralaminar and translaminar 110

fracture toughnesses are considered in the damage variables formulation, enabling the control of the energy dissipation associated with each failure mode. 3.1. Damage onset detection The damage onset is detected by the Hashin-Rotem criteria [27], which are given by Fiber tension (σ1 > 0):

115

σ1 =1 Xt

(24)

|σ1 | =1 Xc

(25)

Fiber compression (σ1 < 0):

Matrix tension (σ2 > 0): 

σ2 Yt

2

σ2 Yc

2

 +

τ12 S

2

τ12 S

2

=1

(26)

=1

(27)

Matrix compression (σ2 < 0): 

 +

where σ1 , σ2 and τ12 represent the state of stress in the principal material coordinate system, Xt is the tensile strength in the fiber direction, Xc the compressive strength in the fiber direction, Yt the tensile 120

strength in the matrix direction, Yc the compressive strength in the matrix direction and S the in-plane shear strength. 3.2. Matrix damage evolution law The matrix damage evolution in composite laminates is primarily related to matrix cracking generated by transverse and in-plane shear loading. However, when fiber breakage occurs, it usually gives rise to cracks

125

in the adjacent matrix portions. Hence, when any of the failure criteria given by Eqs. (24)-(27) are met, matrix damage commences and grows according to the following damage evolution law [24]:

dm = dtm + dcm − dtm dcm

7

(28)

3.2.1. Tensile matrix damage evolution law The tensile matrix damage evolution (2 > 0) is described by the bilinear law [24] dtm =

¯tm,u t ¯m,u − ¯tm,0

¯tm,0 ¯tm

 1−

 (29)

where ¯tm is the resultant strain associated with tensile matrix damage [25], and is given by ¯tm = 130

q 2 22 + γ12

(30)

where 2 and γ12 are the transverse and in-plane shear strains in the principal material coordinate system, respectively. Using polar coordinates, we can write

2 = ¯tm cos φtm




γ12 = ¯tm sin φtm

φtm

 = arccos

(31)




2 ¯tm

(32)

 (33)

The quantity ¯tm,0 represents the damage onset resultant strain associated with tensile matrix damage. Considering that

t2,0 = ¯tm,0 cos φtm




γ12,0 = ¯tm,0 sin φtm 135




(34)

(35)

and also that

t σ2,0 = E2 t2,0

(36)

τ12,0 = G12 γ12,0

(37)

t where t2,0 and γ12,0 are the damage onset strains and σ2,0 and τ12,0 are the damage onset stresses associated

with transverse tension and in-plane shear, respectively, and substituting Eqs. (34) and (35) in Eqs. (36)

8

and (37), and the resultant equations in Eq. (26), we can write

¯tm,0 =

"

E2 cos (φtm ) Yt

2

 +

G12 sin (φtm ) S

2 #−1/2 (38)

The ultimate resultant strain ¯tm,u associated with tensile matrix damage is derived considering the power 140

law [25] 

λ¯

Gtm Gtmc

 +

Gsm Gsmc

λ¯ =1

(39)

¯ = 1 for a unidirectional lamina, Gt and Gs are the intralaminar fracture toughnesses and Gt where λ mc mc m and Gsm are the energy release rates associated with transverse tension and in-plane shear, respectively. In the ultimate state, we can write

Gtm

=l

∗

Z

t2,u

σ2 d2 = l∗

0

Gsm

=l

∗

γ12,u

Z

t σ2,0 t2,u 2

(40)

τ12,0 γ12,u 2

τ12 dγ12 = l∗

0

(41)

where t2,u and γ12,u are the ultimate strains associated with transverse tension and in-plane shear, respec145

tively, and l∗ is the characteristic length, which is, for rectangular elements, expressed as [28] l∗ = lx sin (|θ|) + ly cos (|θ|)

(42)

t2,u = ¯tm,u cos φtm

(43)

where lx and ly are the element sizes. Considering also that




γ12,u = ¯tm,u sin φtm




(44)

and substituting Eqs. (36), (37), (43) and (44) in Eqs. (40) and (41) and the resulting equations in Eq. (39) we have that

¯tm,u = 150

2 ¯tm,0 l∗



E2 cos2 (φtm ) G12 sin2 (φtm ) + Gtmc Gsmc

All the other damage evolution laws are derived in the same way.

9

−1 (45)

3.2.2. Compressive matrix damage evolution law The compressive matrix damage evolution law (2 < 0) is given by ¯cm,u c ¯m,u − ¯cm,0

dcm =

¯cm,0 ¯cm

 1−

 (46)

with

¯cm =

¯cm,0

" =

¯cm,u =

E2 cos (φcm ) Yc 2



¯cm,0 l∗

q 2 22 + γ12

2

 +

(47)

G12 sin (φcm ) S

2 #−1/2

E2 cos2 (φcm ) G12 sin2 (φcm ) + Gcmc Gsmc

φcm

 = arccos

2 ¯cm

(48) −1 (49)

 (50)

where ¯cm , ¯cm,0 , ¯cm,u are the resultant strain, damage onset resultant strain and the ultimate resultant strain 155

associated with compressive matrix damage, respectively, and Gcmc is the intralaminar fracture toughness associated with transverse compression. 3.3. Fiber damage evolution law The fiber damage evolution in composite laminates is directly related to fiber breakage caused by longitudinal loading. Hence, once the failure criteria given by Eqs. (24) and (25) are met, fiber damage commences

160

and grows according to the following evolution law:

df = dtf + dcf − dtf dcf

(51)

3.3.1. Tensile fiber damage evolution law The tensile fiber damage evolution law (1 > 0) is given by dtf =

t1,u t1,u − t1,0

 1−

t1,0 1

 (52)

with

t1,0 =

10

Xt E1

(53)

t1,u =

2Gtf c Xt l∗

(54)

where 1 is the longitudinal strain in the principal material coordinate system and t1,0 , t1,u and Gtf c are the 165

damage onset strain, the ultimate strain and the translaminar fracture toughness associated with longitudinal tension, respectively. 3.3.2. Compressive fiber damage evolution law The longitudinal compressive damage evolution law (1 < 0) is given by dcf =

c1,u c1,u − c1,0

 1−

c1,0 |1 |

 (55)

with Xc E1

(56)

2Gcf c Xc l∗

(57)

c1,0 =

c1,u = 170

where c1,0 , c1,u and Gcf c are the damage onset strain, the ultimate strain and the translaminar fracture toughness associated with longitudinal compression, respectively. Once damage onset is detected, the respective damage variables must be always calculated throughout the analysis. However, in order to account for damage irreversibility, they must be updated only if the new calculated values are greater than the previously calculated ones.

175

3.4. Constitutive matrix degradation The lamina constitutive matrix in the principal material coordinate system [Q] is degraded in each layer of each element using the degradation approach proposed by Van Paepegem and Degrieck [29]. The damaged lamina constitutive matrix in the principal material coordinate system is expressed as



(1 − df ) Q11 p p  [Qd ] =  (1 − df ) (1 − df ) (1 − dm )Q12  0

p p (1 − df ) (1 − df ) (1 − dm )Q12

0

(1 − df ) (1 − dm ) Q22

0

0

    

(1 − df ) (1 − dm ) Q66 (58)

The degradation also takes into account the coupling between fiber failure and matrix cracking [26].

11

180

4. Solution method Applying a small initial perturbation to the panel, the system equations of motion are fully solved in the time domain by employing an iterative form of the Newmark method [30]. No modal reduction is performed in the present formulation. The solution method is illustrated in the flowchart in Fig. 3, and consists of the following steps:

185

1. At the initial time t = 0, prescribe the quantities λ and ca and the time step ∆t, and calculate the matrices [M ], [Ka ], [Ca ]. Set the initial global displacement vector {W }0 as the first in-vacuo natural mode of the panel with w/h = 0.01 at x/a = 0.75 and y/b = 0.5, and the initial global velocity and   ˙ ¨ acceleration vectors W and W to zero; 0 0 2. In all elements, evaluate the Hashin-Rotem criteria in all layers;

190

3. For all elements, compute the damage variables df and dm in the layers in which the Hashin-Rotem criteria were satisfied; 4. Set the iteration i to zero; 5. Use the global displacement vector at time t as a first estimate of the global displacement vector at time t + ∆t;

195

6. Compute the global nonlinear stiffness matrix [KN L ] using the estimated global displacement vector at time t + ∆t and the previously calculated damage variables in conjunction with Eq. (58). 7. Using the Newmark method expressions given in Ref. [30], compute the new estimate of the displacement vector at time t + ∆t and the error; 8. If the error is greater than the tolerance, increment the iteration i and go back to step 6.

200

9. Compute the global velocity and acceleration vectors at time t + ∆t using the Newmark method expressions given in Ref. [30]; 10. If the final time t = tf is reached, end the simulation; 11. Update the global vectors at time t, increment the time t and go back to step 2. This solution method has the advantage of not needing the computation of a tangent stiffness matrix

205

throughout the iterative procedure.

5. Results and discussion In this section, the verification of the aeroelastic and progressive damage models, and simulations in order to investigate the effects of progressive damage on the aeroelastic behavior of different composite panels are presented.

12

210

5.1. Aeroelastic model verification The verification of the nonlinear aeroelastic model was made comparing the limit-cycle oscillation amplitudes obtained by Dixon and Mei [31] for a simply supported square antisymmetric cross-ply graphite-epoxy panel of six layers [0◦ /90◦ ]3 with the ones obtained by the present model. Their results were calculated by employing the linearized updated-mode with a nonlinear time function approximation. The panel dimen-

215

sions are a = b = 305 mm and h = 3.05 mm. The material properties are E1 = 207 GPa, E2 = 5.17 GPa, G12 = 2.59 GPa, ν12 = 0.25 and ρm = 2766 kg/m3 . It was employed ca = 0.01. The comparison for 16 × 16 and 20 × 20 meshes is shown in Fig. 4, where WA represents the limit-cycle amplitude divided by the panel thickness h at the point x/a = 0.75 and y/b = 0.5. The results obtained are in very good agreement with the literature, and clearly demonstrate the accuracy of the present nonlinear

220

finite element model. 5.2. Progressive damage model verification Since experiments in order to validate the predictions of the damage induced by a supersonic flow in composite panels were not possible to be conducted, results for comparison are not available in the literature and the existing commercial codes are unable to simulate the event in study, a quasi-static analysis of the

225

progressive failure of a composite panel subjected to a uniform pressure was conducted in Abaqus/Explicit using Hashin progressive damage model [32] and the results were compared with the ones generated by the present model. The panel used for the verification is a simply supported square symmetric AS4/PEEK panel of twelve layers [0◦ /45◦ /0◦ / − 45◦ /0◦ /90◦ ]S , which represents a typical lay-up used in the aerospace industry. This

230

panel will be called IND panel. The panel dimensions are a = b = 200 mm and h = 1.52 mm. The material properties are E1 = 127.6 GPa, E2 = 10.3 GPa, G12 = 6 GPa, ν12 = 0.32, ρm = 1600 kg/m3 , Xt = 2023 MPa, Xc = 1234 MPa, Yt = 80 MPa, Yc = 160 MPa, S = 82.6 MPa, Gtf c = 128 kJ/m2 , Gcf c = 128 kJ/m2 , Gtmc = 5.6 kJ/m2 , Gcmc = 9.31 kJ/m2 and Gsmc = 4.93 kJ/m2 [33, 34, 35, 36]. In Fig. 5, it is shown the comparison of the load-displacement curves until failure obtained by the two

235

models considering 40 × 40 and 44 × 44 meshes, where Wc = w/h at the center of the panel. The progression of damage until the failure of the structure is dominated by matrix damage generated by tension. With this, the comparison of the tensile matrix damage maps on the verge of failure of the layers 1, 2 and 6, which have the orientations 0◦ , 45◦ and 90◦ , respectively, for a 40 × 40 mesh is given in Fig. 6. The layers are numbered from the bottom to the top layer. It can be seen that both models predict

240

that the failure of the structure occurs due to the emergence of cracks in two parallel edges of the panel. Considering the results, it can be concluded that the present model can accurately predict the progressive failure of a composite panel subjected to bending loads, as in the case of a fluttering panel.

13

5.3. Influence of progressive damage on the aeroelastic behavior The influence of progressive damage on the nonlinear aeroelastic behavior of simply supported square 245

symmetric AS4/PEEK panels of twelve layers was investigated. The panels used in the analyses consist in the IND panel, a cross-ply panel [(0◦ /90◦ )3 ]S and angle-ply panels [(+θ/ − θ)3 ]S . The values of θ that were investigated in the case of the angle-ply panels were 0◦ , 30◦ , 45◦ , 60◦ and 90◦ . The cross-ply panel was called CP panel and the angle-plies, APθ panels. The panel dimensions and properties are the same of the IND panel. In all simulations, λ = 800 and ca = 0.01 were employed.

250

The AP60 panel is the one that presents the most complex aeroelastic response amongst the analyzed panels. With this, a convergence study for the undamaged and damaged conditions was conducted for the panel AP60 and the same meshes were adopted in the analysis of the other panels. The undamaged condition represents the simulation that does not take into account the progressive damage model, and the damaged condition represents the simulation in which the progressive damage model is included. In Fig.

255

7, it is shown the time histories comparison for the undamaged condition considering 28 × 28 and 32 × 32 meshes. The comparison of the time histories for the damaged condition is illustrated in Fig. 8, where 40 × 40, 44 × 44 and 48 × 48 meshes were employed. The stabilized matrix damage maps comparison of the top and bottom layers for the three meshes is depicted in Fig. 9. The top and bottom layers are the ones in which the damage is more critical due to the higher bending strains. Considering the results, the panels

260

were modeled with a 28 × 28 mesh in the undamaged condition, while in the damaged condition, a 40 × 40 mesh was employed. The aeroelastic response of the IND panel is illustrated in the time history and phase plane plot in Figs. ˙ = w/h 10 and 11, respectively, where W = w/h and W ˙ at the point x/a = 0.75 and y/b = 0.5. As can be seen, the panel undergoes a limit cycle oscillation. The CP panel presents a similar behavior with slightly

265

higher amplitudes due to the absence of 45◦ and −45◦ layers. The time history and phase plane plot of the CP panel are shown in Figs. 12 and 13, respectively. For the IND and CP panels, damage initiation is not detected. The panel AP0 presents a unstable behavior characterized by a flutter condition which gives rise to a small amount of matrix damage in the layers. The comparisons between the time histories and phase plane

270

plots of the undamaged and damaged conditions are given in Figs. 14 and 15. As can be seen, the aeroelastic behavior remains almost unaltered after the damage initiation. The matrix damage progression until its stabilization is illustrated in the damage maps of the top and bottom layers at different times in Fig. 16. No fiber damage is detected in this scenario. The aeroelastic behavior of the panels AP30 and AP45 also consists in a limit cycle oscillation. However,

275

damage onset is not present in these cases. The time histories and phase plane plots of these panels are shown in Figs. 17-20. 14

The panel AP60 undergoes a pronounced nonperiodic oscillation. Matrix damage takes place and propagates through the layers, changing the aeroelastic behavior of the panel, but still without the occurrence of catastrophic failure. The comparisons between the time histories and phase plane plots of the undamaged 280

and damaged conditions are depicted in Figs. 21 and 22. It is interesting to note that the stiffness reduction in the matrix direction caused by the damage does not change the amplitudes of the panel substantially. However, it completely changes the form of the orbits in the phase plane plot, giving rise to an even more nonperiodic behavior. The matrix damage maps are shown in Fig. 23. Fiber damage is not present in this case.

285

Finally, the panel AP90 suffers a very rapid catastrophic failure. The comparisons between the time histories and phase plane plots of the undamaged and damaged conditions presented in Figs. 24 and 25 show a substantial increase in the panel amplitude prior to the failure. The damage maps showing the damage progression until the failure of the panel are given in Fig. 26, and reveal the formation of an extensive matrix damaged region parallel to the fibers which divides the panel in two parts. Despite the

290

damage extent, the stiffness in the fiber direction also remains intact in all layers. The maximum and minimum stresses achieved throughout the simulations were evaluated and are given in Table 1. As can be seen, for all panels, the magnitudes of the maximum and minimum achieved stresses in the fiber direction are far below the strengths. The panels IND, CP, AP30 and AP45 are far from presenting damage initiation of any kind. For the angle-ply panels, except for the panel AP0, the stresses in the matrix

295

direction and in-plane shear stresses increase with the misalignment between the fibers and the airflow. Given that analyses of composite laminates are normally very configuration dependent, general considerations are difficult to be drawn, and further investigations are usually necessary. However, some expected tendencies were present in the results and some interesting aspects of the analysis should be highlighted: • All the presented damage maps show that damage is concentrated around the point x/a = 0.75 and

300

y/b = 0.5, which is roughly the point of maximum amplitude of a fluttering panel [37, 38, 39]; • As it would be expected, the results indicate that matrix is the critical component of fluttering composite panels made of unidirectional layers, since only matrix damage has been present and the fibers have remained intact in all the simulations. However, concerning matrix damage, it is interesting to note that great part of the panel has also remained intact;

305

• In the case of the angle-ply panels, except for the panel AP0, the results indicate that the increase in the misalignment between the fibers and the airflow is directly related to a more unstable aeroelastic behavior, with higher amplitudes, velocities and induced stresses in the matrix direction and in-plane shear stresses. The more unstable motion tends to damage or cause the catastrophic failure of the panel, giving rise to changes in the aeroelastic behavior. The more unstable aeroelastic behavior and 15

310

the occurrence of matrix damage onset in the panel AP0 is probably due to the absence of fibers in the transversal direction; • Results like the ones obtained by the present analysis can be used for the design of damage tolerant panels. The dimensions and materials of the panel can be chosen in a rational way, aiming a certain maximum level of damage in the structure. Additionally, since the localization of damage is known, a

315

convenient positioning of stringers and the application of variable thickness is possible. Some further considerations concerning limitations of the present model and future works should also be made: • In the present work, the value of the nondimensional dynamic pressure was fixed at 800, which is a relatively high magnitude for this quantity. However, most of the panels analyzed in this study did

320

not present a damage onset. The occurrence of damage initiation was restricted to angle-ply panels, which are not so usual in the aerospace industry. In order to investigate the progressive failure of a wider range of composite fluttering panels, analyses at higher dynamic pressures must be conducted. In this case, the inclusion of nonlinear aerodynamic effects through the use of the third-order piston theory is mandatory, since these analyses will be carried out at the hypersonic regime [40, 41];

325

• The analyses in this study were restricted to square panels of 0.2 m length and 1.52 mm thickness. In the nondimensionalized panel flutter problem, for a given pair of λ and ca , the relation w/h in the coordinates x/a and y/b is kept constant independently of the panel length [37]. With this, since most of the panels of 0.2 m length do not present damage initiation, larger square panels with the same thickness and in the same conditions will also not present damage onset due to their lower bending

330

stiffness. In order to study the failure of larger panels, their thickness must be increased. However, an increasing in the panel thickness can lead to the violation of the thin panel hypothesis in some situations. In this case, the first-order laminate plate theory should be implemented in the structural model [18]. The use of the first-order laminate plate theory when combined with cohesive zone interface elements [42] also allows the prediction of delaminations, since out-of-plane shear stresses are taken

335

into account in the model; • The vastly observed fiber bridging phenomenon cannot also be predicted by the present model. An additional fiber bridging fracture energy must be included in the damage evolution laws. For this purpose, the approach proposed by H¨ ower et al. [43] could be adopted; • The characteristic length in smeared crack progressive damage models must be lower than a maximum

340

value which is directly proportional to the translaminar and intralaminar fracture toughness of the material [24]. The aeroelastic analysis considering progressive failure of panels made of composite 16

materials with very low values of these properties can become prohibitive, since the element size necessary to conduct the analysis can be very small in these circumstances. For large panels, the problem is more critical, because more elements are necessary for the modeling. The panels analyzed 345

in this work are made of a thermoplastic resin with a high intralaminar fracture toughness, which has significantly reduced the computational cost; • If distorted quadrangular or triangular elements are needed for the modeling, a more general procedure for the calculation of the characteristic length of the element must be used, since the determination of the crack propagation direction is not so straightforward in these cases. For this purpose, the

350

objectivity algorithm employed in the work of Donadon et al. [24] could be used.

6. Conclusions A finite element model for predicting the nonlinear aeroelastic behavior of composite panels undergoing intralaminar and translaminar progressive damage in supersonic flow has been presented. The classical plate theory in conjunction with the von K´ arm´an nonlinear strains is used for structural modeling, and the 355

linear piston theory is used to model the aerodynamic loads. Progressive damage is modeled by a smeared cracking formulation in which stress-based, continuum damage mechanics and fracture mechanics approaches are combined. No modal reduction is performed and an iterative form of the Newmark method is used for the numerical direct integration in time of the nonlinear equations. Simulations considering different layups are conducted, in which the influence of progressive damage on the aeroelastic behavior of the panels is

360

investigated, and damage extent and failure mechanisms are assessed. Limitations of the model and future works are also discussed. The results obtained in the analyses consist in important insights concerning the flutter-induced damage in composite panels and in relevant information for conducting a damage tolerance design.

7. Acknowledgments 365

This work was supported by the Coordination for the Improvement of Higher Education Personnel (CAPES), grant number 88882.180833/2018-01, the Funding Authority for Studies and Projects (FINEP), grant number 0114018300 and the National Council for Scientific and Technological Development (CNPq), grant number 301053/2016-2.

References 370

[1] E. H. Dowell, Panel flutter - a review of the aeroelastic stability of plates and shells, AIAA Journal 8 (3) (1970) 385–399. doi:10.2514/3.5680.

17

[2] C. Mei, K. Abdel-Motagaly, R. Chen, Review of nonlinear panel flutter at supersonic and hypersonic speeds, Applied Mechanics Reviews 52 (10) (1999) 321–332. doi:10.1115/1.3098919. [3] W.-H. Chen, H.-C. Lin, Flutter analysis of thin cracked panels using the finite element method, AIAA Journal 23 (5) 375

(1985) 795–801. doi:10.2514/3.8986. [4] K.-J. Lin, P.-J. Lu, J.-Q. Tarn, Flutter analysis of anisotropic panels with patched cracks, Journal of Aircraft 28 (12) (1991) 899–907. doi:10.2514/3.46115. [5] R. M. V. Pidaparti, Free vibration and flutter of damaged composite panels, Composite Structures 38 (1) (1997) 477–481. doi:10.1016/S0263-8223(97)00082-2.

380

[6] K. Wang, D. J. Inman, C. R. Farrar, Crack-induced changes in divergence and flutter of cantilevered composite panels, Structural Health Monitoring 4 (4) (2005) 377–392. doi:10.1177/1475921705057977. [7] S. Natarajan, G. Manickam, S. Bordas, Supersonic flutter analysis of functionally graded material plates with cracks, Frontiers in Aerospace Engineering 2 (2) (2013) 91–97. [8] N. A. Abdullah, J. L. Curiel-Sosa, M. Akbar, Aeroelastic assessment of cracked composite plate by means of fully coupled

385

finite element and doublet lattice method, Composite Structures 202 (2018) 151–161. doi:10.1016/j.compstruct.2018. 01.015. [9] L.-C. Shiau, Flutter of composite laminated beam plates with delamination, AIAA Journal 30 (10) (1992) 2504–2511. doi:10.2514/3.11253. [10] A. A. Yazdi, J. Rezaeepazhand, Structural similitude for flutter of delaminated composite beam-plates, Composite Struc-

390

tures 93 (7) (2011) 1918–1922. doi:10.1016/j.compstruct.2011.02.004. [11] J. Fazilati, Panel flutter of curvilinear composite laminated plates in the presence of delamination, Journal of Composite Materials 52 (20) (2018) 2789–2801. doi:10.1177/0021998318754641. [12] R. M. V. Pidaparti, C. C. Chang, Finite element supersonic flutter analysis of skewed and cracked composite panels, Computers & Structures 69 (2) (1998) 265–270. doi:10.1016/S0045-7949(98)00003-0.

395

[13] T. Kim, S. N. Atluri, R. G. Loewy, Modeling of microcrack damaged composite plates undergoing nonlinear bimodular flutter oscillations, AIAA Journal 36 (4) (1998) 598–606. doi:10.2514/2.411. [14] H. AsadiGorgi, M. Dardel, M. H. Pashaei, Effects of all-over part-through cracks on the aeroelastic characteristics of rectangular panels, Applied Mathematical Modelling 39 (23) (2015) 7513–7536. doi:10.1016/j.apm.2015.03.017. [15] E. H. Dowell, Fatigue life estimation of fluttering panels, AIAA Journal 8 (10) (1970) 1879–1881. doi:10.2514/3.6008.

400

[16] D. Y. Xue, C. Mei, Finite element nonlinear flutter and fatigue life of two-dimensional panels with temperature effects, Journal of Aircraft 30 (6) (1993) 993–1000. doi:10.2514/3.46444. [17] T. W. Strganac, Y. I. Kim, Aeroelastic behavior of composite plates subject to damage growth, Journal of Aircraft 33 (1) (1996) 68–73. doi:10.2514/3.46904. [18] J. N. Reddy, Mechanics of Laminated Composite Plates and Shells: Theory and Analysis, CRC Press, 2004.

405

[19] S. G. P. Castro, T. A. M. Guimar˜ aes, D. A. Rade, M. V. Donadon, Flutter of stiffened composite panels considering the stiffener’s base as a structural element, Composite Structures 140 (2016) 36–43. doi:10.1016/j.compstruct.2015.12.056. [20] X. Guo, C. Mei, Using aeroelastic modes for nonlinear panel flutter at arbitrary supersonic yawed angle, AIAA Journal 41 (2) (2003) 272–279. doi:10.2514/2.1940. [21] F. K. Bogner, R. L. Fox, L. A. Schmit, Jr., The generation of inter-element-compatible stiffness and mass matrices by the

410

use of interpolation formulas, in: Proceedings of the Conference on Matrix Methods in Structural Mechanics, 1966, pp. 397–443. [22] O. O. Ochoa, J. N. Reddy, Finite element analysis of composite laminates, in: Solid Mechanics and its Applications, Springer, 1992. doi:10.1007/978-94-015-7995-7. [23] M. V. Donadon, The structural behavior of composite laminates manufactured using resin infusion under flexible tooling,

18

415

Ph.D. thesis, Imperial College London (2005). [24] M. V. Donadon, L. Iannucci, B. G. Falzon, J. M. Hodgkinson, S. F. M. de Almeida, A progressive failure model for composite laminates subjected to low velocity impact damage, Computers & Structures 86 (11) (2008) 1232–1252. doi: 10.1016/j.compstruc.2007.11.004. [25] M. V. Donadon, S. F. M. de Almeida, M. A. Arbelo, A. R. de Faria, A three-dimensional ply failure model for composite

420

structures, International Journal of Aerospace Engineering 2009. doi:10.1155/2009/486063. [26] N. O. Yokoyama, M. V. Donadon, S. F. M. de Almeida, A numerical study on the impact resistance of composite shells using an energy based failure model, Composite Structures 93 (1) (2010) 142–152. doi:10.1016/j.compstruct.2010.06.006. [27] Z. Hashin, A. Rotem, A fatigue failure criterion for fiber reinforced materials, Journal of Composite Materials 7 (4) (1973) 448–464. doi:10.1177/002199837300700404.

425

[28] Z. P. Baˇ zant, Mechanics of Fracture and Progressive Cracking in Concrete Structures, Springer, 1984, pp. 1–94. doi: 10.1007/978-94-009-6152-4_1. [29] W. Van Paepegem, J. Degrieck, Modelling damage and permanent strain in fibre-reinforced composites under in-plane fatigue loading, Composites Science and Technology 63 (5) (2003) 677–694. doi:10.1016/S0266-3538(02)00257-9. [30] K.-J. Bathe, Finite Element Procedures, Klaus-J¨ urgen Bathe, 2016.

430

[31] I. R. Dixon, C. Mei, Finite element analysis of large-amplitude panel flutter of thin laminates, AIAA Journal 31 (4) (1993) 701–707. doi:10.2514/3.11606. [32] Abaqus 6.14-1 Documentation. [33] J. F. Chen, E. V. Morozov, K. Shankar, A combined elastoplastic damage model for progressive failure analysis of composite materials and structures, Composite Structures 94 (12) (2012) 3478–3489. doi:10.1016/j.compstruct.2012.04.021.

435

[34] D. Carlile, D. Leach, D. Moore, N. Zahlan, Mechanical Properties of the Carbon Fiber/PEEK Composite APC-2/AS-4 for Structural Applications, ASTM, 1989, pp. 199–212. doi:10.1520/STP24603S. [35] S. T. Pinho, P. P. Camanho, M. F. de Moura, Numerical simulation of the crushing process of composite materials, International Journal of Crashworthiness 9 (3) (2004) 263–276. doi:10.1533/ijcr.2004.0287. [36] Z. Zhang, G. Hartwig, Low-temperature viscoelastic behavior of unidirectional carbon composites, Cryogenics 38 (4)

440

(1998) 401–405. doi:10.1016/S0011-2275(98)00022-8. [37] E. H. Dowell, Nonlinear oscillations of a fluttering plate., AIAA Journal 4 (7) (1966) 1267–1275. doi:10.2514/3.3658. [38] E. H. Dowell, Nonlinear oscillations of a fluttering plate. II., AIAA Journal 5 (10) (1967) 1856–1862. doi:10.2514/3.4316. [39] M. A. Kouchakzadeh, M. Rasekh, H. Haddadpour, Panel flutter analysis of general laminated composite plates, Composite Structures 92 (12) (2010) 2906–2915. doi:10.1016/j.compstruct.2010.05.001.

445

[40] T. Bein, P. Friedmann, X. Zhong, I. Nydick, Hypersonic flutter of a curved shallow panel with aerodynamic heating, in: 34th Structures, Structural Dynamics and Materials Conference, 1993. doi:10.2514/6.1993-1318. [41] G. Cheng, C. Mei, Finite element modal formulation for hypersonic panel flutter analysis with thermal effects, AIAA Journal 42 (4) (2004) 687–695. doi:10.2514/1.9553. [42] P. A. A. E. Mendes, M. V. Donadon, Numerical prediction of compression after impact behavior of woven composite

450

laminates, Composite Structures 113 (2014) 476–491. doi:10.1016/j.compstruct.2014.03.051. [43] D. H¨ ower, B. A. Lerch, B. A. Bednarcyk, E. J. Pineda, S. Reese, J.-W. Simon, Cohesive zone modeling for mode I facesheet to core delamination of sandwich panels accounting for fiber bridging, Composite Structures 183 (2018) 568–581. doi:10.1016/j.compstruct.2017.07.005.

19

Figure 1: Panel model.

Figure 2: BFS element.

20

Figure 3: Solution method flowchart.

21

1 0.9 0.8 0.7

WA

0.6 0.5 0.4 0.3 0.2 Present work - 16 # 16 Mesh Present work - 20 # 20 Mesh Dixon and Mei

0.1 0 165

180

195

210

225

240

255

270

285

6 Figure 4: Comparison of limit-cycle amplitudes for a simply supported square [0◦ /90◦ ]3 graphite-epoxy panel.

4.5

4

3.5

p (M P a)

3

2.5

2

1.5

1 Present Present Abaqus Abaqus

0.5

work - 40 # 40 Mesh work - 44 # 44 Mesh - 40 # 40 Mesh - 44 # 44 Mesh

0 0

1

2

3

4

5

6

7

8

9

Wc Figure 5: Comparison of load-displacement curves until failure, Panel IND.

22

(a) Layer 1, θ = 0◦ , Present work.

(b) Layer 1, θ = 0◦ , Abaqus.

(c) Layer 2, θ = 45◦ , Present work.

(d) Layer 2, θ = 45◦ , Abaqus.

(e) Layer 6, θ = 90◦ , Present work.

23

(f) Layer 6, θ = 90◦ , Abaqus.

Figure 6: Tensile matrix damage maps comparison, 40 × 40 Mesh, Panel IND.

5 4 3 2

W

1 0 -1 -2 -3 -4 -5 0.49

28 # 28 Mesh 32 # 32 Mesh

0.492

0.494

0.496

0.498

0.5

t (s) Figure 7: Time history convergence study, Undamaged condition, Panel AP60.

5 4 3 2

W

1 0 -1 -2 -3 -4 -5 0.09

40 # 40 Mesh 44 # 44 Mesh 48 # 48 Mesh

0.092

0.094

0.096

0.098

0.1

t (s) Figure 8: Time history convergence study, Damaged condition, Panel AP60.

24

(a) Top layer - 40 × 40 Mesh.

(b) Top layer - 44 × 44 Mesh.

(c) Top layer - 48 × 48 Mesh.

(d) Bottom layer - 40 × 40 Mesh.

(e) Bottom layer - 44 × 44 Mesh.

(f) Bottom layer - 48 × 48 Mesh.

Figure 9: Matrix damage convergence study, Panel AP60.

25

2

1.5

1

W

0.5

0

-0.5

-1

-1.5

-2 0.99

0.992

0.994

0.996

0.998

1

t (s) Figure 10: Time history, Panel IND.

12000

9000

6000

_ W

3000

0

-3000

-6000

-9000

-12000 -2

-1.5

-1

-0.5

0

0.5

1

W Figure 11: Phase plane plot, Panel IND.

26

1.5

2

2.5 2 1.5 1

W

0.5 0 -0.5 -1 -1.5 -2 -2.5 1.49

1.492

1.494

1.496

1.498

1.5

t (s) Figure 12: Time history, Panel CP.

12000

9000

6000

_ W

3000

0

-3000

-6000

-9000

-12000 -2.5

-2

-1.5

-1

-0.5

0

0.5

1

W Figure 13: Phase plane plot, Panel CP.

27

1.5

2

2.5

2.5 2 1.5 1

W

0.5 0 -0.5 -1 -1.5 -2

Undamaged Damaged

-2.5 0.49

0.492

0.494

0.496

0.498

0.5

t (s) Figure 14: Time history, Panel AP0.

15000 12000 9000 6000

_ W

3000 0 -3000 -6000 -9000 -12000 -15000 -2

Undamaged Damaged

-1.5

-1

-0.5

0

0.5

1

W Figure 15: Phase plane plot, Panel AP0.

28

1.5

2

(a) Top layer, t = 0.05 s.

(b) Top layer, t = 0.1 s.

(c) Top layer, t = 0.5 s.

(d) Bottom layer, t = 0.05 s.

(e) Bottom layer, t = 0.1 s.

(f) Bottom layer, t = 0.5 s.

Figure 16: Matrix damage maps, Panel AP0.

29

2

1.5

1

W

0.5

0

-0.5

-1

-1.5

-2 0.49

0.492

0.494

0.496

0.498

0.5

t (s) Figure 17: Time history, Panel AP30.

10000 8000 6000 4000

_ W

2000 0 -2000 -4000 -6000 -8000 -10000 -2

-1.5

-1

-0.5

0

0.5

1

W Figure 18: Phase plane plot, Panel AP30.

30

1.5

2

2.5 2 1.5 1

W

0.5 0 -0.5 -1 -1.5 -2 -2.5 0.49

0.492

0.494

0.496

0.498

0.5

t (s) Figure 19: Time history, Panel AP45.

16000

12000

8000

_ W

4000

0

-4000

-8000

-12000

-16000 -2.5

-2

-1.5

-1

-0.5

0

0.5

1

W Figure 20: Phase plane plot, Panel AP45.

31

1.5

2

2.5

5 4 3 2

W

1 0 -1 -2 -3 -4

Undamaged Damaged

-5 0.49

0.492

0.494

0.496

0.498

0.5

t (s)

40000

40000

30000

30000

20000

20000

10000

10000

_ W

_ W

Figure 21: Time history, Panel AP60.

0

0

-10000

-10000

-20000

-20000

-30000

-30000

-40000 -5

-4

-3

-2

-1

0

1

2

3

4

-40000 -5

5

-4

-3

W

-2

-1

0

W

(a) Undamaged.

(b) Damaged. Figure 22: Phase plane plot, Panel AP60.

32

1

2

3

4

5

(a) Top layer, t = 0.003 s.

(b) Top layer, t = 0.006 s.

(c) Top layer, t = 0.5 s.

(d) Bottom layer, t = 0.003 s.

(e) Bottom layer, t = 0.006 s.

(f) Bottom layer, t = 0.5 s.

Figure 23: Matrix damage maps, Panel AP60.

33

3.5 3 2.5 2

W

1.5 1 0.5 0 -0.5 -1

Undamaged Damaged

-1.5 0

0.0003

0.0006

0.0009

0.0012

0.0015

t (s) Figure 24: Time history, Panel AP90.

25000

20000

15000

10000

_ W

5000

0

-5000

-10000

-15000 Undamaged Damaged

-20000 -1

-0.5

0

0.5

1

1.5

2

2.5

W Figure 25: Phase plane plot, Panel AP90.

34

3

3.5

(a) Top layer, t = 0.0014 s.

(b) Top layer, t = 0.0015 s.

(c) Top layer, t = 0.001571 s.

(d) Bottom layer, t = 0.0014 s.

(e) Bottom layer, t = 0.0015 s.

(f) Bottom layer, t = 0.1571 s.

Figure 26: Matrix damage maps, Panel AP90.

35

Table 1: Maximum and minimum stresses in MPa.

Panel IND CP AP0 AP30 AP45 AP60 AP90 Strength

σ1max 528 572 354 413 651 1217 408 2023

σ1min -310 -253 -193 -278 -379 -600 -158 -1234

σ2max 31 45 86 34 64 96 96 80

36

σ2min -19 -16 -57 -22 -29 -50 -151 -160

max τ12 27 33 46 15 34 64 82 82.6

min τ12 -28 -33 -44 -17 -37 -62 -82 -82.6