Remaining useful life prediction of laminated composite materials using Thermoelastic Stress Analysis

Composite Structures 210 (2019) 381–390

Contents lists available at ScienceDirect

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

Remaining useful life prediction of laminated composite materials using Thermoelastic Stress Analysis ⁎

Ricardo Marquesa, Mustafa Unelb,c,d, Mehmet Yildizb,c,d, , Afzal Sulemane,

T

⁎

a

Instituto Superior Técnico, Universidade de Lisboa, Av. Rovisco Pais, 1049-001, Lisbon, Portugal Integrated Manufacturing Technologies Research and Application Center, Sabanci University, Tuzla, 34956 Istanbul, Turkey c Composite Technologies Center of Excellence, Sabanci University-Kordsa, Istanbul Technology Development Zone, Sanayi Mah. Teknopark Blvd. No: 1/1B, Pendik, 34906 Istanbul, Turkey d Faculty of Engineering and Natural Sciences, Sabanci University, Tuzla, 34956 Istanbul, Turkey e Department of Mechanical Engineering, University of Victoria, Victoria, BC, Canada b

A R T I C LE I N FO

A B S T R A C T

Keywords: Remaining useful life Polymer matrix composite materials Thermoelastic Stress Analysis Strain energy density

Estimation of the remaining useful life (RUL) based on the elastic strain energy density using the non-contacting and full-field capabilities of Thermoelastic Stress Analysis (TSA) was applied to E-glass fiber-reinforced polymers (GFRPs) subjected to high-cycle fatigue loading conditions. A quantitative prediction was successfully attained during the material’s intermediate and final fatigue life phases. The initial state of the specimen was used as the reference data for defining a damage parameter that enabled the detection of visually imperceptible subsurface cracks.

1. Introduction Aircraft structures made of fiber reinforced polymeric (FRP) composite materials provide higher specific strength and stiffness together with improved fatigue life and corrosion resistance when compared to metallic-based components. Hence, FRP composites have found extensive application in the latest Airbus A-350 and Boeing 787 airframes containing over 50% of weight content [1] due to the achievable cost saving objectives in fuel consumption and maintenance operations. During the mission, the aircraft fuselage is subjected to several fatigue loading scenarios, namely pressurization during climb and descent, gust loads on the wing and inertial forces due to maneuvers [2]. The aforementioned loads coupled with external environmental factors and manufacturing defects enhance damage propagation and growth before ultimately leading to catastrophic failure. To circumvent the risk of imminent failure of the structures subjected to harsh operational loading and environmental conditions, and avoid the requirement of scheduled inspection and maintenance causing prolonged down-time cost, Structural Health Monitoring (SHM) concepts have been investigated as a viable approach to ensure the structural reliability of the components through continuous data collection and processing with pertinent sensors systems, hence replacing the scheduled inspection and maintenance with condition based ones. Nevertheless, the engineering tendency of using traditional non-

⁎

destructive inspection (NDI) techniques and strict regulatory demands within the framework of periodic monitoring of aircraft structures still preclude the widespread application of SHM techniques [3]. Yet, Thermoelastic Stress Analysis (TSA) based on the physics of thermoelastic effect stands out as a valid alternative to more common NDI and SHM methods. Its full-field, non-contacting and non-destructive nature provides quantitative stress/strain data without both a thorough surface preparation and a priori knowledge of prone-to-damage locations. The thermoelastic effect in solids takes place due to the mechanicalthermal coupling of the material deformed under adiabatic and reversible conditions. The reversible and adiabatic deformation can be achieved if the material’s stress–strain state is within the linear elastic domain and high loading frequencies are imparted so that heat is, neither transferred to the environment, nor between laminae in a laminated composite material [4]. John Gough’s experiments reported in 1803 using Indian rubber (or caoutchouc) were a pioneering study in the qualitative observation of the thermoelastic phenomenon [5]. The mathematical relation between principal stresses (Δσ ) and the surface temperature difference (ΔT ) measured between the corresponding loading states in isotropic materials was originally derived by William Thomson around 1850 [5]:

ΔT = −

Corresponding author. E-mail addresses: [email protected] (M. Yildiz), [email protected] (A. Suleman).

https://doi.org/10.1016/j.compstruct.2018.10.047 Received 25 August 2018; Accepted 18 October 2018 Available online 14 November 2018 0263-8223/ © 2018 Published by Elsevier Ltd.

α T0 Δσ , ρC

(1)

Composite Structures 210 (2019) 381–390

R. Marques et al.

where ρ is the density, α is the coefficient of thermal expansion (CTE), C is the specific heat capacity at constant pressure (all assumed to be temperature independent) and T0 is the absolute surface temperature measured at the mean point of the cyclic temperature or stress, here onward, referred to as the mean absolute surface temperature. The very small change in surface temperature due to cyclic loading, ΔT , corresponding to temperatures at maximum and minimum stresses, have only become detectable with the commercial availability of highly sensitive photon detectors. In the literature, the relevant applications of TSA to FRP composite materials mainly include qualitative and quantitative damage assessment [6,7] and the interpretation of the measured thermoelastic signal source [8–11]. Moreover, the potential of using TSA’s output to directly predict fatigue life was evaluated to be an emerging research topic by Wong [12]. Emery and Dulieu-Barton [13] correlated the evolution trends of both Young’s Modulus and first strain invariant with the imminent failure scenario of glass fiber-reinforced polymers (GFRPs) under fatigue. Horn et al. [14] devised an experiment to obtain a localized and single valued stress concentration factor (SCF) for a specimen subjected to impact loading utilizing the ratio of temperature differences ΔT , measured near and away from impact-damaged regions of the specimen. Specimens with an initial impact damage were subjected to a load-controlled fatigue test. They processed the experimental results by plotting the modified stress (the applied stress on the specimen multiplied by the corresponding SCF) versus the number of fatigue cycles until failure, and the outcome revealed similitude with a SN curve. In addition, Johnson et al. [15] applied the probabilistic Markov Chain Method (MCM) to determine the upper and lower bounds of a stochastic S-N curve for a FRP single lap joint. However, to the best of the authors knowledge, none of the relevant studies in the literature made a quantitative remaining useful life (RUL) prediction exploiting the full-field nature of TSA. There is no previous study that defines the damage parameter based on the strain energy density at the mean strain level, which enables the detection of the damage evolution, most of the time not visible with conventional imaging methods or the naked eye. The main objective of this work is to develop a TSA-based SHM methodology that can predict the number of fatigue cycles until failure of a FRP with a quasi-isotropic (QI) configuration. Given that in any structure under loading, the strain is a quantity that remains constant along the laminate’s thickness and easily measurable, a promising and readily usable Remaining Useful Life (RUL) methodology should preferably be based on the strain field as an input parameter. To this end, in this study, as the RUL model, we have used strain energy release rate (SERR) of FRP composite subjected to strain-controlled fatigue test [16]. Keulen et al. [17] utilized the SERR concept in a SHM technique that relied on embedded Fiber Bragg Grating (FBG) sensors. Due to their local measurement capability, the embedded FBG sensors were reported to be capable of revealing the onset and growth of matrix cracks and other micron- or submicron-level defects undetectable by an extensometer [17]. However, the fragile nature of an FBG sensor reveals to be disadvantageous during the manufacturing process since a FBG sensor requires special handling and implementation in a highly plausible location for damage growth. Overcoming these limitations is made possible when using TSA’s full-field and non-contacting nature. In this study, the outcome of a TSA analysis is presented where the RUL of the specimen is determined at different steps each comprising a fixed number of fatigue cycles. Moreover, it is shown that a local damage parameter independent of the total number of performed cycles can be constructed. The damage parameter can be used to determine the critical locations at which the structure is likely to fail and therefore it is an important information for developing a reliable scheduled aircraft maintenance scenario. Both the RUL and the damage parameter information are valuable tools to enhance condition based maintenance.

2. Theoretical background 2.1. Stress-based thermoelastic equations The thermoelastic effect detected on the surface of multidirectional laminates under in-plane stress conditions can be expressed in relation to the principal stress ( x , y ) in the global coordinate and principal material directions (1, 2 ) in the local coordinate in matrix form:

ΔT = −

T0 T [α ]Tx , y [Δσ ]x , y = − 0 [α ]T1,2 [Δσ ]1,2 ρCp ρCp

(2)

where [σ ]x , y and [α ]x , y are respectively the stress and CTE vectors in the global coordinate with the components of σx , σy, τs and α x , α y, αs while [σ ]1,2 and [α ]1,2 are respectively the stress and CTE vectors in the local coordinate with the components of σ1, σ2, τ6 , and α1, α2, α6 . When loading the composite material along the principal stress and principal material directions, τs and α6 are null, respectively. Therefore, the expansion of Eq. (2) yields the well-known TSA surface stress-based formulation valid for orthotropic materials:

ΔT = −

T0 ⎛ T0 ⎛ ⎞ ⎞ ⎜α1 Δσ1 + α2 Δσ2 ⎟ = − ⎜α x Δσx + α y Δσy ⎟ ρC ⎝ ρC ⎠ ⎠ ⎝

(3)

For isotropic materials, CTE is independent of the considered direction, namely, α x = α y = αs = α . Therefore, for the principal stress directions, Eq. (3) is simplified as:

ΔT = −

T0 ⎛ α ⎜Δσx + Δσy ⎞⎟ ρC ⎝ ⎠

(4)

2.2. Strain-based thermoelastic equations In uniaxially loaded multidirectional laminates, a different stress field is generated in each layer depending upon the fiber orientation, thus locally distinct stress-induced temperature changes are attained. However, if perfectly bonded layers are considered, the strain field along (x , y ) direction will be constant throughout the laminate’s thickness. Therefore, writing the thermoelastic equation as a function of strain is a more useful approach. A uniaxially loaded symmetric and balanced laminated composite material deforms with a null shear strain, γs = 0 . If the Poisson’s ratio νxy - defined as the ratio of strains along y and x directions - is used, the transverse strain component can be eliminated in the forthcoming equations. Therefore, the in-plane strain vector for the principal stress directions, [ε ]x , y is represented as:

[ε ]x , y = ( εx − νxy εx 0 )T

(5)

In general, the composite manufacturing process leads to the accumulation of epoxy resin in the outer surface of the laminate, immediately after the first ply. During the transmission phenomena of the heat wave generated in the laminate’s first ply through the epoxy layer, the wavelength and amplitude can be reduced with the increase in loading frequency and epoxy thickness, thus leading to the thermal drag-down phenomena [18], which can be detected from an increasing phase shift of the measured signal i.e. a phase angle different from zero degrees [19]. The phase shift is described as the angle difference between the strain and temperature signals where the strain signal is provided by the fatigue testing machine. The temperature changes being measured become a sole result of the epoxy’s thermoelastic effect that, thus acting as a strain-witness of the laminate, i.e. Δεxr = Δεxc and Δεyr = Δεyc = −νxy Δεxc , where the subscripts r and c stand for resin and composite, respectively. References [8,20] respectively reported that the minimum resin thickness thresholds of 25 μ m and 30 μ m are required to achieve the desired signal attenuation in glass and carbon fiber-reinforced composites. Composite materials with low-conductive reinforcements (in 382

Composite Structures 210 (2019) 381–390

R. Marques et al.

particular, glass fibers) or with the surface ply oriented transversely to the applied load proved to better fit the resin-rich layer (RRL) model (or strain-witness model) [19,21,22]. The stress–strain relation for isotropic materials reads:

Er ⎛ ⎞ Δσx + Δσy = ⎜1 − νxy ⎟ Δεxc 1 − νr ⎝ ⎠

incorporating the fatigue life within phases I and III in the RUL prediction. Thus, the total number of cycles until failure, Nt , can be written as:

Nt = (6)

(11)

Another valid strain energy density calculation arises from the definition of specific heat capacity at constant strain (or constant volume), Cε = ∂U / ∂T . Due to the fact that the epoxy resin is considered to be responsible for the measured thermoelastic temperature and incompressible, Cε is equal the specific heat capacity at constant pressure, Cr , namely, Cε = Cr [23]. Considering a linear dependency of Cr with the local mean absolute temperature, T0, p , the energy density for discrete locations of the specimen’s surface, u 0, p , can be alternatively written as:

where Er is the resin’s Young’s Modulus. Given that the resin is an isotropic material, the strain-based thermoelastic equation for the RRL model can be obtained through substituting Eq. (6) into Eq. (4):

ΔT α ⎡ Er ⎛ ⎤ ⎞ =− r ⎢ ⎜1 − νxy ⎟ Δεxc ⎥ T0 ρr Cr 1 − νr ⎝ ⎠ ⎦ ⎣

u 0 (r − 1) 0.75[a (εmax / εult )b ]

(7)

,where ΔT / T0 will be subsequently referred to as normalized temperature, T . From here onwards, the term Δεxc will be replaced by Δεx .

u 0, p = ρ . Cr (T0, p). ΔTp

2.3. Remaining useful life (RUL) calculations

where ΔTp is the temperature difference at pixel p determined during the TSA stage.

Natarajan et al. [16] showed that the expended SERR, dU / dN , of FRP’s subjected to constant strain amplitude conditions under fatigue is characterized by three distinct phases where U and N respectively represent strain energy and fatigue cycle number. In the first fatigue phase comprising the initial 15–20% of fatigue life, the rapid formation and interconnection of matrix cracking causes a sharp, non-linear increase in the expended strain energy. The second phase accounts for 15–20% to 90% of the fatigue life where there is a gradual and linear increase in the expended strain energy, which is attributed to propagation of several cracks, fiber debonding and delamination. The final phase is distinguished by a sharp nonlinear increase in the expended strain energy due to the plurality of fiber breakages. The same authors developed a promising method based on the linear evolution of the SERR during the second phase and found a power law relation between the slopes (dU / dN ) and applied strain ratios (εmax / εult ) such that:

2.4. Temperature correction methodology As the fatigue process continues, matrix cracking and some other types of defects such as fiber debonding and delamination form and evolve, hence leading to additional temperature increases in the specimen due to the frictional effects occurring at the defect surfaces. To compensate this non-stress related temperature increase in the specimen, the temperature correction methodology reported in [24] is employed. A calibration parameter Rp is determined for each pixel in the structure as the ratio between the harmonic average of the mean absolute surface temperature of the specimen at the undamaged state, < T0(0) >, and T0, p :

Rp = (< T0(0) > / T0, p )n ,

⎜

⎟

(8)

where εult and εmax are respectively the ultimate strain of the specimen and the maximum strain applied to the specimen under the fatigue loading, a and b are material coefficients which depend on the stacking sequence, fiber volume fraction, specimen thickness and loading type (e.g. tension–tension, tension–compression, bending). Integrating the previous equation between the moment before fatigue initiation and the end of the second phase, the following RUL equation can be obtained:

Nf =

U0 (r − 1) , a (εmax / εult )b

Tp, c = Rp Tp

1 Ex Δεx2, 2

(14)

Any further reference to a normalised temperature must always undergo this correction procedure, therefore Tp, c will be simply referred to as T . 3. Experimental procedure and implementation

(9)

3.1. GFRP’s fabrication and processing

where r = Uf / U0, U0 and Uf are respectively the strain energy at the initial mean load level and at the end of phase II, and Nf is the number of cycles up to the end of phase II. To avoid dependency of the RUL calculations on the specimen’s volume selected for strain averaging, hereinafter the strain energy density, u is to be considered. From the reversible condition inherent to TSA, the fatigue mean strain level is limited by the stress–strain linear elastic condition. Therefore, the elastic strain energy density at the mean strain level, u 0 , can be calculated as:

u0 =

(13)

where the coefficient n depends upon the infra-red wavelength range detectable by the camera. The corrected normalised temperature for a given pixel p in the structure, Tp, c , is obtained from the measured normalized temperature at each pixel p , Tp , as it follows:

b

dU ε = a ⎛ max ⎞ , dN ⎝ εult ⎠

(12)

A quasi-isotropic (QI) configuration ([0/45/ −45/90]s ) was selected for further RUL prediction tests (hereinafter referred to as QI(0)). Here, the fiber direction 0 refers to the axial direction. Composite laminates were produced through employing vacuum bagging of glass fiber prepregs on an aluminum mold with the dimensions of 400 × 450 mm. The prepreg was fabricated by Kordsa Company (with the code KOM10 EGF UD330) composed of stitched unidirectional E-glass fabric with a total areal weight of 330 g/m2 from Metyx Composites and the epoxy resin with the code OM10. The complete vacuum bagging assembly consisted of peel-ply and breather layers enclosed by a vacuum bag sealed along the mold edges and connected to an air suction pump through a vacuum valve with quick disconnect coupling and hose. The aluminum mold surface was prepared with sealer and release agents before stacking. The vacuum bag curing process included a first dwell period of 1h30 at 70 °C followed by a second dwell period of 2 h at 120 °C inside a temperature-controlled oven. Finally, the assembly cooled down to room temperature at a slow rate to avoid generating

(10)

where Ex is the specimen’s longitudinal Young’s Modulus and Δεx is the axial global strain of the test specimen obtained using Eq. (7). Given that the fatigue coefficients in Eq. (9) are obtained using the information from phase II of the expended strain energy versus cycle number curve - corresponding to 75% of the specimen total real fatigue life - it is proposed to correct Eq. (9) with a factor of 0.75, thereby 383

Composite Structures 210 (2019) 381–390

R. Marques et al.

thermal residual stresses.

Table 1 GFRP’s and epoxy’s material properties.

3.2. Preparation and curing of epoxy resin samples

Property

Value

Py Pult (KN) εy εult (% ) ρ ρr (kg/m3) Cr (J/(kg.°C)) v v (% )

The matrix material for Kordsa’s prepreg is the hot melt resin system composed of structural epoxy resin, dicyandiamide hardener and urone-based accelerator. The ratio between epoxy and hardener was kept at 1 whereas 3–5 parts per hundred (phr) of accelerator for 100 phr of resin were mixed. To produce mechanical test specimens for the matrix material, the epoxy resin was initially heated up to 50 °C inside an oven and then mixed with room temperature hardener and accelerator under vacuum. Subsequently, the mixture was cast inside a preheated metal mold at 80 °C with a dog bone shape cavity. Curing cycle consisted of a first step at 120 °C for 1 h followed by a second step at 140 °C for 2 h.

v f (% ) v m (% )

Property

Value

3.72 17.86

Er (GPa)

3.08

0.39 2.59

νr

0.35

1834 1142 1290 [4.71,6.97] [46.77,49.03]

ab αr (/°C) α x (/°C) αy (/°C)

4619.5 7.4298 6.60 × 10−5 1.20 × 10−5 1.20 × 10−5

46.26

resin as Cr (T ) = 4.6T + 1176.4 . The Cr value presented in Table 1 was obtained for T = < T0(0) >. Density measurements were based on the hydrostatic weighing method using water as the immersion liquid. Five GFRP samples and one epoxy resin sample with the dimensions of 20 × 25 mm were firstly weighed in air before suspending in water. The reduction in weight after immersion – equal to the sample’s volume - was determined. Subsequently, the matrix burnoff test (ASTM D3171) was performed to determine the void content, v v as well as fiber and matrix volume fractions, v f and v m . Since upper and lower bounds for the reinforcement’s density were known, a range of v f and v v values were therefore determined. To remove any water content absorbed by the specimens during the density measurement, samples were conditioned inside a temperature-controlled oven for 24 h at 70 °C. Having recorded their initial mass, samples were placed inside different crucibles and held in a nitrogen-purging furnace for 1 h at 560 °C and the crucible’s content was remeasured afterwards. Relevant material properties obtained through the characterization methods described in the preceding parts for the GFRP and epoxy resin are presented in Table 1. Three tension–tension strain-controlled cyclic tests were ran on QI (0) specimens at each considered strain ratio (0.4, 0.5, 0.6, 0.7 and 0.8) until failure. The frequency of all fatigue tests was chosen to be 5 Hz to avoid the adverse influence of autogeneous heating in the fatigue life of GFRP’s [16,17]. A minimum strain level of 8% of εult assured positive load values until failure as the specimen’s axial stiffness degraded during the test’s course. Using the data collected, the fatigue coefficients, a and b, in the Eq. (11) were obtained by fitting a power law curve to the strain energy density variation, du/ dN (obtained between 15% and 90% of the fatigue life) versus strain ratio data, εmax / εult , as shown in Fig. 1a. A similar power law trend for the ratio between the strain energy densities at the end of the fatigue phase II and initial mean strain level (r = uf / u 0) was observed in Fig. 1b. The r value for the strain ratio imparted during the fatigue stage was obtained from the fitting curve’s equation. A summary of the adopted strategy for the implementation of RUL prediction is now presented and schematically shown in Fig. 2.

3.3. Material characterization tests The mechanical properties of laminate’s, namely, the longitudinal strains and loads at the yield (εy and Py , respectively) and at the breakage (εult and Pult , respectively) points were determined in accordance with the ASTM D3039 standard using an Instron 8802 servo hydraulic testing system equipped with a ±250 kN load cell. A total of 5 GFRP tensile specimens were cut into the dimensions of 25 × 250 × 2 mm and bonded with 25 × 50 × 1.5 mm aluminum tabs using epoxy Araldite 2011. A quasi-static load under constant upper grip head displacement of 2 mm/min was continuously applied until failure while acquiring axial strain εx employing an Instron 2620-603 dynamic extensometer with 50 mm gauge length. For determining tensile properties of epoxy resin system, at least two dog-bone specimens were prepared without tabbing and tested following the ASTM D638 standard using an Instron 5900 electromechanical testing system with the test conditions described above for the testing of the composite laminate. Load and strain data collected by an Instron 2650-563 biaxial extensometer with 50 mm axial and 25 mm transverse gauge lengths for εx values between 0.1% and 0.3% were employed to calculate Young’s modulus (Er ) and the Poisson’s ratio (νr ). The tests were terminated for a εx equal to 0.31%. Experiments for determining the CTE of specimens were performed inside a temperature-controlled oven using a strain gage method described in [25,26]. 20 × 25 mm samples made of isotropic materials (epoxy resin, and aluminum as the reference material) were bonded with a Kyowa KFG 120Ω uniaxial strain gage each whereas one equally dimensioned QI(0) GFRP sample had two mutually perpendicular strain gages bonded. A K-type thermocouple was also attached to the reference aluminum sample whose CTE value had been determined to be 2.40 × 10−5 °C−1 through performing a Thermomechanical Analysis (TMA) using a Mettler Toledo TMA/SDTA 1 STARe System. To enable free dilatation and contraction of samples, they were placed on top of a low-friction medium. Temperature increased at the rate of 1 °C/min from 20 °C to 50 °C and strain data was acquired for the temperature range between 25 °C and 45 °C. Tests were considered as valid after attaining five consecutively equal linear regressions of the difference between the measured strain in the material of interest (either the GFRP or epoxy sample) and the reference aluminum specimen. The specific heat capacity of the matrix material at constant pressure, Cr was measured through following the Differential Scanning Calorimetry (DSC) method described in ASTM E1269. The Mettler Toledo DSC 3+ device underwent an initial heat flow calibration using a synthetic sapphire disk whose C is known over a range of temperatures. An epoxy sample with the weight of 8.50 mg was heated up under a constant heat flow rate. For every 5 °C increase in temperature, the difference between the heat flow transmitted to the sample and reference sapphire was continuously recorded. A linear curve was fitted to the Cr data collected between 10 °C and 50 °C, thereby obtaining a temperature dependent specific heat capacity relation for the epoxy

• Determination of the global Young’s modulus (E ) and Poisson’s ratio (ν ); • ‘TSA’ stage: cyclic strain-controlled test under tensile-tensile condix

xy

• •

tions with a maximum applicable strain level limited by the yield point threshold defined in Table 1. ΔTp and T0, p data were collected; ‘Fatigue’ stage: cyclic strain-controlled test under tensile-tensile conditions where damage was allowed to initiate and evolve in the laminate; The procedure restarted from the 1st stage of Fig. 2 if failure had not occurred. Nt was continuously updated in a step-by-step strategy.

4. Results and discussion 4.1. RUL prediction The RUL procedure can be generalized in global steps consisting of 384

Composite Structures 210 (2019) 381–390

R. Marques et al.

Table 2 Specimen initial dimensions and strain controlled-test parameters for the TSA and Fatigue stages. Dimensions h (mm) w (mm)

2.08 25.09

TSA stage f (Hz) N εmin (%) εmax (%)

Fatigue stage 12 1200 0.052 0.388

f (Hz) N εmin (%) εmax (%)

5 20000 0.052 0.725

a maximum εx of 0.31% to update the νxy and Ex where strains were measured using the biaxial extensometer. The forthcoming TSA stage was performed under strain-controlled conditions by using the dynamic extensometer attached to the back side of the specimen previously in contact with the peel-ply while curing. To maintain the position of extensometer, two elastics were positioned in the upper and bottom corners of the extensometer’s gauge section. Cyclic strain limits for the TSA stage were set to the minimum and mean of the fatigue strain levels (2–28% of εult ), namely 2–15% of εult . After performing several signal-frequency dependency tests, it was determined that a 12 Hz TSA frequency guaranteed an almost stabilized normalized temperature at the surface. Thermal images during the TSA stage were captured using a thermal camera (FLIR X6580sc model) with an Indium antimonide (InSb) photon detector. This thermal camera can detect the radiation with the wavelength ranging between 1.5 and 5.5 μ m, thus yielding a coefficient n = 9.74 given in Eq. (13). The thermal camera installed with a 25 mm lens was positioned at a stand-off distance of 570 mm with a sidewise offset in relation to the specimen’s surface under observation to avoid the Narcissus effect. The integration time (i.e. the inverse of the frame rate) of 1900 μ s was employed so that the thermal camera could measure the surface temperature of the specimen up to approximately 50 °C if required. During the TSA stage, the Edevis DisplayImg 6 software was used for the data acquisition and processing as well as the visualization of the phase angle and normalized temperature results in a pixelwise format. Using the dynamic extensometer’s sinusoidal strain signal as a reference, all the thermoelastic data collected corresponding to the maximum and minimum values of the entire strain cycles (i.e. 1200 TSA cycles) was integrated together with the Edevis DisplayImg 6 software. The specimen was tensioned until reaching the mean strain level before initiating the TSA cycles, while data collection was

Fig. 1. Experimental results and fitting plot for (a) du/ dN = f (εmax / εult ) = 4619.5(εmax / εult )7.4298 (b) εmax /εult = f (r ) = 0.0845. r 2.2292 .

three main stages that were sequentially repeated for every 20,000 fatigue cycles as can be recalled from Fig. 2 until notable crack/damage in the structure was observed. The notable damage was regarded as the one which caused the fatigue strain limits from Table 2 to be exceeded although a complete separation of the specimen did not occur. The specimen was mounted on an Instron 8802 servohydraulic machine with a load cell capability of ±250 kN. An initial quasi-static tensile test was performed in accordance with the recommendations from standard ASTM D3039. The specimen was loaded within its elastic domain up to

Fig. 2. Schematic of the experimental procedure for determining the remaining useful life using the strain-based methodology. 385

Composite Structures 210 (2019) 381–390

R. Marques et al.

the maximum radius of influence centered in the pixel of interest R = 3h , where Δx is the physical distance between pixels defined as 1/ np , being np the total number of pixels along the specimen’s width direction. The schematic representation of the smoothing process is shown in Fig. 4 where i refers to the pixel of interest and j identifies its neighbouring pixels located inside the domain of influence. Assuming that step 0 is the material’s undamaged state, this methodology has only been applied from step 1 onwards. The effect of the smoothing process on the normalised temperature at step 14, Ti(14) is presented in Fig. 3c. In comparison with Fig. 3b, regions with higher Ti(s) - representing the crack fronts - had their amplitudes reduced as a consequence of the influence from the surrounding crack surfaces. For each pixel, the normalized temperature Ti(s) after smoothing was rejected whenever it was higher than the original data. After completing the smoothing operation, a new global < T g(s) > value at step s was obtained using Eq. (15). A total of 17 steps - equivalent to 327,630 fatigue cycles - were performed. From the total number of fatigue cycles performed, the different fatigue life phases could be identified. Fatigue phase-I extended up to step 2, fatigue phase-II was confined between steps 3 and 14 while the remaining two steps were attributed to the fatigue phaseIII of the material’s fatigue life. For determining the correction factor, the reference temperature < T0(0) > was calculated to be 297.9 K. In Table 3 are presented the quasi-static tensile test results (namely, νxy and Ex ), the harmonic averaged normalized temperature < T g(s) > obtained before and after smoothing (simplified as BS and AS, respectively) as well as the longitudinal strain change calculated using the RRL model in Eq. (7) and the elastic strain energy density at the mean strain level, u 0 in Eq. (10), both dependent on < T g(s) >. In comparison with the extensometer measured strain change during the TSA cycle given in Table 2, the calculated strain value using the normalized temperature information at the undamaged state yield a difference of 4.5%, computed using the relation, ∣ΔεxE − ΔεxT ∣ × 100/ΔεxT , thus stating the reliability of this model’s predictions where ΔεxT and ΔεxE are the strain differences based on the extensometer measurement and normalized temperature information respectively. (see Table 4). A single Nt for each step had been defined based on u 0 . Its step-bystep evolution calculated through using the BS and AS data given in Table 3 is shown in Fig. 5. It is worthy to note that the evolution of Nt clearly resembles to the fatigue phases of a FRP material based on the plot of elastic modulus versus number of cycles to failure [16]. This is quite understandable owing to the fact that the strain energy density is linearly dependent on the elastic modulus. To validate the proposed RUL methodology, the calculated results were compared against a reference curve. The reference curve was drawn by subtracting the performed number of cycles from the total fatigue life. RUL predictions within the initial steps up to the step 4 show some notable deviations from the reference curve, which can be attributed to the random and nongradual damage formation and failure mechanisms including matrix cracking, delamination and to a certain extent fiber failures occurring during the composite material’s initial life [29]. The damage formation leads to a nonlinear RUL decay witnessed from the decrease in Young’s modulus from Table 3, uncompensated by the proposed RUL prediction model. In general, it is rather challenging to compensate the decay of stiffness in the RUL prediction models due to the requirement of stress information, which is ordinarily unavailable. As the step number increases towards a more progressive damage evolution stage of the fatigue life, the RUL predictions based on BS and AS data yield the most accurate results between steps 4 and 9 where a maximum relative error with respect to the reference result (in absolute value) of 22% is calculated. Here, it should be noted that at the initial steps (step 0 and 1), the smoothing process does not notably affect the Nt value since the cracks are rather localized and their effect cannot be extended to broader areas with smoothing process. As the test evolves onto the fatigue phase-II and approaches the last phase, and more specifically from step 9 onwards, the theoretical

manually triggered 60 s afterwards and ended after counting 1200 TSA cycles. The emissivity value of 0.95 was used in the software for TSA. The methodology for measuring the emissivity value used in this study was described clearly in the MSc thesis of the first author [27]. A first global normalized average temperature < T g(s) > at step s was obtained using a weighed harmonic averaging method:

1 1 = . p < T g(s) >

p

∑ i=1

1 Ti(s)

, (15)

where the subscript g stands for global value, and p is the total number of pixels. The weighed harmonic averaging augments the influence of locations with lower normalized temperature in the RUL calculation unlike the arithmetic averaging. Having completed the first set of TSA cycles, referred to as the step 0, the specimen was subjected to the first strain-controlled fatigue stage with a maximum strain value of 28% of the εult for 20,000 cycles with a frequency of 5 Hz. The corresponding r value for the strain ratio of 0.28 can be determined to be 1.71 from Fig. 1.Table 2 summarises the TSA and fatigue test conditions as well as specimen’s dimensions. The elastic strips used to fix the position of the dynamic extensometer led to a reduction in the available area for RUL monitoring. Therefore, temperature data extraction was split into three rectangular areas as shown in Fig. 4. u 0 calculations (Eq. (10)) for each TSA stage relied on updated single νxy and Ex values taken as representative of the whole structure. Therefore, it is a prudent approach to use a single global Δεx based on < T g(s) > for the calculation of u 0 . To better analyze the experimental data for predicting the RUL of the tested material, we performed data processing on T such as smoothing with weighted harmonic averaging and incorporated the effect of crack formation into calculation steps for predicting the Nt . One may speak of two sources of temperature change within the structures under the external cyclic loading. Namely, one is due to the deformation or thermoelastic effect while the second one is associated with friction on the crack surfaces upon the formation of damage. These two sources can be separated from each other referring to the phase shift data since any friction between crack surfaces will give rise to a phase angle significantly bigger than zero degrees. Stated otherwise, if the temperature change is related with the thermoelastic effect, it should follow the cyclic nature of the deformation. Moreover, one can separate a damage or crack into two regions according to their thermoelastic response, namely the crack surface and the immediate vicinity of the crack front. The material discontinuity between crack surfaces leads to a reduction in the elongation level of the crack surfaces in response to a cyclic external load or deformation, yielding lower T and Nt , while the frictional effects between the crack surfaces result in a phase shift. On the other hand, at the crack fronts, the enhanced tendency for the material to separate leads to higher deformation levels represented by locally higher T values that consequently lead to an unrealistic increase in the fatigue lives. As the specimen remains intact and continuous at the crack fronts, deformations on these locations occur at the same testing frequency, hence the phase angle remains in-phase. The pixelwise comparison of the phase angle and T given in Figs. 3a and 3b clearly indicates that, the higher the phase shift, the lower the deformation, which is likely to happen at crack surfaces. The yellowish spots on the phase shift image directly correspond to lower T on the temperature image, implying that these regions are where damage and crack exist in the specimen. Furthermore, regions with low T are immediately accompanied by those with high T , bespeaking that these juxtaposed low and high T regions are cracks and crack fronts. To circumvent and minimize overprediction of the fatigue lives at the crack fronts, the thermoelastic response at the crack surfaces was extended into the crack fronts. To this end, we utilized a smoothing or interpolation methodology using a two-dimensional quintic spline kernel function [28] considering the smoothing length as h = 2Δx and 386

Composite Structures 210 (2019) 381–390

R. Marques et al.

Fig. 3. Datasets at step 14 (280,000 cycles performed) (a) Phase angle (b) T before applying the smoothing process (c) T after applying the smoothing process.

Table 3 Static test measurements together with the TSA-derived weighed harmonic averaged results before and after smoothing (respectively denoted as BS and AS).

Fig. 4. Schematical representation of the data collection and smoothing process.

results start to deviate from the reference curve by over predicting the RUL. This can be attributed to the fact that as the number of fatigue cycles increases, more localized damages form and grow while still maintaining a significant part of the structure able to withstand the applied deformations. Although their influence on the averaged or global strain calculation based on the TSA formulation is limited since the area they occupy is very small in comparison to the intact region, their effect in the final life of the structure is notable. In addition, the over-predicted results could be also contributed by the non-updated matrix material properties (i.e., αr , ρr , Cr , νr , Er ) in the strain-based

s

Ex (GPa)

νxy

BS⁎

AS⁎

BS⁎⁎

AS⁎⁎

BS⁎⁎⁎

AS⁎⁎⁎

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

17.92 15.93 15.53 15.39 15.17 15.05 14.99 14.66 14.61 14.59 14.65 14.46 14.45 14.25 14.12 13.89 13.19

0.31 0.29 0.26 0.28 0.27 0.28 0.27 0.26 0.25 0.25 0.27 0.25 0.27 0.26 0.26 0.24 0.22

5.14 5.07 5.08 5.05 5.10 5.18 5.04 4.91 5.03 4.88 4.87 5.17 4.97 5.25 5.05 4.97 4.87

5.14 4.80 4.80 4.78 4.81 4.89 4.75 4.61 4.74 4.60 4.58 4.89 4.66 4.92 4.73 4.64 4.53

3.51 3.36 3.23 3.31 3.29 3.39 3.25 3.12 3.16 3.06 3.14 3.25 3.21 3.34 3.21 3.08 2.94

3.51 3.18 3.06 3.13 3.10 3.20 3.06 2.94 2.97 2.89 2.96 3.07 3.01 3.13 3.01 2.88 2.74

11.0 9.02 8.12 8.41 8.21 8.66 7.94 7.15 7.29 6.85 7.24 7.63 7.44 7.96 7.30 6.60 5.71

11.0 8.07 7.26 7.52 7.31 7.72 7.03 6.33 6.46 6.09 6.41 6.82 6.54 6.99 6.41 5.75 4.93

⁎

(×10−4) . (×10−3) . ⁎⁎⁎ (×10 4) . ⁎⁎

387

u 0 [J/m3]

Δεx



Composite Structures 210 (2019) 381–390

R. Marques et al.

Table 4 CIF as a function of fatigue step. s

0

1

2

3

4

5

Γs Ns xs CIF

1.00 0 0.00 1.00

0.93 20000 0.07 0.96

0.93 40000 0.14 0.92

0.91 60000 0.21 0.87

0.90 80000 0.28 0.81

0.86 100000 0.34 0.68

s Γs Ns xs CIF

6 0.84 120000 0.41 0.59

7 0.84 140000 0.48 0.54

8 0.84 160000 0.55 0.48

9 0.85 180000 0.62 0.46

10 0.85 200000 0.69 0.42

11 0.86 220000 0.76 0.43

s Γs Ns xs CIF

12 0.81 240000 0.83 0.28

13 0.77 260000 0.90 0.18

14 0.77 280000 0.97 0.15

15 0.76 300000 1.03 0.12

16 0.75 320000 1.10 0.09 (s ) (s ) Fig. 6. Tmin and Tmax before smoothing as a function of the step number.

threshold temperature, are assumed to be potentially experiencing damage in the form of crack. Looking for all the pixels with normalized temperature below Tt , and subsequently taking their harmonic average, one can define a new temperature quantity referred to as damage temperature for each step s, Td(s) = < T p(s) < Tt >. The CIF can be now introduced into the RUL formulation 11 through the fraction Γs = Td(s)/ Tt . As more damage forms, Γs will consequently decrease. (s ) A power law growing trend regarding the maximum Tmax values can be also inferred from Fig. 6. This increase reflects the severity of existing crack fronts in the specimen and resembles the evolution of du/ dN as a function of εmax / εult represented in Fig. 1a. To incorporate the crack front propagation phenomena as a function of the accumulated number of fatigue cycles into the CIF, the exponent of Γs had the fatigue coefficient, b, included together with xs , where xs = Ns / Nu is the fraction of the fatigue cycles at each TSA step, Ns is the total number of cycles up to step s and Nu is equivalent to Nt defined at the undamaged state using the weighted harmonic average of normalized temperature data before smoothing. The final form of the RUL prediction equation containing the CIF yields:

Fig. 5. Comparison between the reference and calculated Nt .

thermoelastic formulation given in Eq. (7) as a function of the accumulated damage. Given that the cracks or damages are quite localized, AS data can not notable improve the predicted RUL. Therefore, a new approach is required to extend the influence of this localized damage regions into intact regions, thereby incorporating the effect of damages on the life of the structure. To this end, after smoothing, a crack intensification factor (CIF) was implemented to spread the effect of crack surfaces towards the intact regions thereby correcting the departure of TSA-based RUL results from the experimental trend. The local minimum and maximum normalized (s ) (s ) temperatures before smoothing (Tmin and Tmax respectively) are plotted as a function of fatigue step number in Fig. 6. In a perfectly homogeneous structure, the minimum and maximum temperature should be (s ) the same at the undamaged state. A nearly linear trend of Tmin is presented which drops down by 42% at the last step before failure. Since the lowest temperature can be associated with the most severe damage, this set of data quantifies the effects of permanent defects in the (0) (16) / Tmin structure’s degradation level. It is very interesting to note that Tmin is almost equal to the r ratio for the applied strain ratio during the (s ) fatigue stage, indicating that the evolution of Tmin as a function of step number is directly related to the life of the structure due to the formation and growth of damage. On the other hand, the initial maximum normalized temperature (0) (i.e. Tmax ) is representative of the most intact region of the domain. Assuming an equal ΔT decay in every point of the structure, regions (0) with a normalized temperature below Tt = 0.42 × Tmax , named as the

Nt =

(r − 1) u 0 s . Γ bx s 0.75[a (εmax / εult )b ]

(16)

The calculated Nt using Eq. (16) is plotted in Fig. 5. Having defined the CIF significantly improved the RUL predictions throughout the end of phase-II, namely starting at step 11, and phase-III where damage growth is predominant and inevitable. To incorporate the benefits from each of the aforementioned approaches, the raw data together with the smoothed data with and without CIF formed the upper and lower bounds of the RUL of the specimen. To be able to obtain a single value for each step, we performed an harmonic average of these three datasets and included our final prediction in Fig. 5. A second approach to predict the RUL of the FRP is based on incorporating the minimum value of pixelwise energy density calculated using Eq. (12) into Eq. (16) and plotted in Fig. 5. Nt calculated with this method is independent of stiffness or Poisson’s ratio measurements. The linear variation of Cr with T0 considering the resin-rich layer’s influence had been determined from the DSC tests. However, the temperature correction methodology is only applicable to the temperature ratio rather than to its separate components, hence limiting the calculations to the raw data. Although over-predicted, results near phase-III follow a decreasing linear slope trend in close agreement with the experiments. 4.2. Damage parameter For step s, a non-dimensional parameter, Dp(s) , was defined by the 388

Composite Structures 210 (2019) 381–390

R. Marques et al.

Fig. 7. Dp(s) at s = 4, 8, 12 and 16 considering the epoxy resin’s material properties.

ratio between the local number of cycles until failure, Nt(,sp) , and the harmonic averaged RUL determined at the undamaged reference state, < Nt(0) >:

Dp(s) = 1 −

Nt(,sp) < Nt(0) >

=1−

Cr(,sp) ΔT p(s) Cr ( < T0(0) > ) < ΔT (0) >

,

Acknowledgments The financial and infrastructural support provided by Sabanci University – Integrated Manufacturing Technologies Research and Application Center and Sabanci University – Kordsa Composite Technologies Center of Excellence is greatly acknowledged. Afzal Suleman acknowledges the NSERC Canada Research Chair funding.

(17)

where < ΔT (0) > is the harmonic averaged temperature change at the undamaged state. The purpose of this parameter was to provide a feasible comparison in a real maintenance scenario where the actual number of cycles until failure is unknown. Full-field plots for steps 4, 8, 12 and 16 together with the specimen’s final configuration are shown in Fig. 7. Damage visually perceived is demarcated inside a red box and detectable by an increased D value at the same location namely during the last performed step. Moreover, subsurface damage is also detected thus highlighting the TSA’s capability to assess damage through-the-thickness.

References [1] Kassapoglou C. Design and analysis of composite structures: with applications to aerospace structures. John Wiley & Sons; 2013. [2] Broek D, Smith SH, Rice RC. Generation of spectra and stress histories for fatigue and damage tolerance analysis of fuselage repairs, Tech. rep. BATTELLE MEMORIAL INST COLUMBUS OH COLUMBUS LABS; 1991. [3] Rajic N, Galea S. Thermoelastic stress analysis and structural health monitoring: an emerging nexus. Struct Health Monit 2015;14(1):57–72. https://doi.org/10.1177/ 1475921714548936. [4] Wong A. A non-adiabatic thermoelastic theory for composite laminates. J Phys Chem Solids 1991;52(3):483–94. https://doi.org/10.1016/0022-3697(91)90180-8. [5] Stanley P. Beginnings and early development of thermoelastic stress analysis. Strain 2008;44(4):285–97. https://doi.org/10.1111/j.1475-1305.2008.00512.x. [6] Fruehmann R, Dulieu-Barton J, Quinn S. Assessment of fatigue damage evolution in woven composite materials using infra-red techniques. Compos Sci Technol 2010;70(6):937–46. https://doi.org/10.1016/j.compscitech.2010.02.009. [7] Kakei A, Epaarachchi J, Islam M, Leng J, Rajic N. Detection and characterisation of delamination damage propagation in woven glass fibre reinforced polymer composite using thermoelastic response mapping. Compos Struct 2016;153:442–50. https://doi.org/10.1016/j.compstruct.2016.06.044. [8] Emery T, Dulieu-Barton J, Earl J, Cunningham P. A generalised approach to the calibration of orthotropic materials for thermoelastic stress analysis. Compos Sci Technol 2008;68(3):743–52. https://doi.org/10.1016/j.compscitech.2007.09.002. [9] Pitarresi G, Galietti U. A quantitative analysis of the thermoelastic effect in cfrp composite materials. Strain 2010;46(5):446–59. https://doi.org/10.1111/j.14751305.2009.00660.x. [10] Pitarresi G, Found M, Patterson E. An investigation of the influence of macroscopic heterogeneity on the thermoelastic response of fibre reinforced plastics. Compos Sci Technol 2005;65(2):269–80. https://doi.org/10.1016/j.compscitech.2004.07.008. [11] Sambasivam QS, Dulieu-Barton SJ. Identification of the source of the thermoelastic response from orthotropic laminated composites. 17th International conference on composite materials (ICCM17). 2009. p. 11. [12] Wong AK. Making the invisible visible: joining the dots from lord kelvin to fighter jet fatigue. Proc. 28th congress of the international council of the aeronautical sciences. 2012. p. 1–10. [13] Emery T, Dulieu-Barton J. Thermoelastic stress analysis of damage mechanisms in composite materials. Compos Part A 2010;41(12):1729–42. https://doi.org/10. 1016/j.compositesa.2009.08.015. [14] Horn GP, Mackin TJ, Kurath P. Estimating the residual fatigue lifetimes of impactdamaged composites using thermoelastic stress analysis. Polym Compos

5. Conclusions A novel TSA-based RUL estimation using the strain energy density was developed and successfully implemented on an initially undamaged quasi-isotropic E-glass fiber-reinforced composite laminate. The initial rapid material degradation visible in the global Ex dataset, as well as the formation and growth of cracks in very small localized regions compared with the overall specimen dimension led to a deviation of the calculated Nt results from the reference linear curve within the initial and final steps. To improve the RUL predictions, the specimen’s life-limiting effect was extended caused by the existing cracks into the RUL formulation through a crack intensification factor. This parameter contains the fatigue coefficient b as a result of the observed similarity (s ) between the evolution of Tmax – in close relation with the crack fronts formation - and du/ dN as a function of εmax / εult . It is shown that this method improved the final prediction considerably. Moreover, a damage parameter based on the evolution of minimum value of pixelwise strain energy density has been proposed. Unaware of the total number of fatigue cycles, the TSA full-field capability enabled the determination of the severity and quantification of the crack responsible for final failure using a damage parameter. 389

Composite Structures 210 (2019) 381–390

R. Marques et al.

2001;22(3):420–31. https://doi.org/10.1002/pc.10548. [15] Johnson S, Wei B, Haj-Ali R. A stochastic fatigue damage model for composite single lap shear joints based on markov chains and thermoelastic stress analysis. Fatigue Fract Eng Mater Struct 2010;33(12):897–910. https://doi.org/10.1111/j. 1460-2695.2010.01531.x. [16] Natarajan V, Gangarao HV, Shekar V. Fatigue response of fabric-reinforced polymeric composites. J Compos Mater 2005;39(17):1541–59. https://doi.org/10. 1177/0021998305051084. [17] Keulen CJ, Akay E, Melemez FF, Kocaman ES, Deniz A, Yilmaz C, et al. Prediction of fatigue response of composite structures by monitoring the strain energy release rate with embedded fiber bragg gratings. J Intell Mater Syst Struct 2016;27(1):17–27. https://doi.org/10.1177/1045389X14560358. [18] Quinn S, Dulieu-Barton J. Identification of the sources of non-adiabatic behaviour for practical thermoelastic stress analysis. J Strain Anal Eng Des 2002;37(1):59–71. https://doi.org/10.1243/0309324021514835. [19] Salerno A, Costa A, Fantoni G. Calibration of the thermoelastic constants for quantitative thermoelastic stress analysis on composites. Rev Sci Instrum 2009;80(3). https://doi.org/10.1063/1.3090885. 034904. [20] Zhang D, Enke N, Sandor B. Thermographic stress analysis of composite materials. Exp Mech 1990;30(1):68–73. https://doi.org/10.1007/BF02322705. [21] Bakis CE, Reifsnider KL. The adiabatic thermoelastic effect in laminated fiber composites. J Compos Mater 1991;25(7):809–30. https://doi.org/10.1177/ 002199839102500702. [22] Pitarresi G, Conti A, Galietti U. Investigation on the influence of the surface resin

[23]

[24]

[25] [26]

[27]

[28]

[29]

390

rich layer on the thermoelastic signal from different composite laminate lay-ups. Appl Mech Mater, 3. Trans Tech Publ; 2005. p. 167–72. https://doi.org/10.4028/ www.scientific.net/AMM.3-4.167. Pitarresi G, Patterson E. A review of the general theory of thermoelastic stress analysis. J Strain Anal Eng Des 2003;38(5):405–17. https://doi.org/10.1243/ 03093240360713469. Dulieu-Barton J, Emery T, Quinn S, Cunningham P. A temperature correction methodology for quantitative thermoelastic stress analysis and damage assessment. Meas Sci Technol 2006;17(6):1627. https://doi.org/10.1088/0957-0233/17/6/ 047. Measurement of thermal expansion coefficient using strain gages, Technical Report, Vishay micro measurements. Zanjani JSM, Al-Nadhari AS, Yildiz M. Manufacturing of electroactive morphing carbon fiber/glass fiber/epoxy composite: process and structural monitoring by fbg sensors. Thin-Walled Struct 2018;130:458–66. Marques REMC. Remaining useful life prediction of laminated composite materials using thermoelastic stress analysis [Master’s thesis]. Instituto Superior Técnico, Universidade de Lisboa; 2017. Shadloo M, Rahmat A, Yildiz M. A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant newtonian fluid. Comput Mech 2013;52(3):693–707. https://doi.org/10.1007/ s00466-013-0841-z. Talreja R. Damage and fatigue in composites–a personal account. Compos Sci Technol 2008;68(13):2585–91.