Gigliotti-L-2016-PhD-Thesis (2)

Multiscale analysisof damage-tolerant

composite sandwich structures

by

Luigi Gigliotti

Department of AeronauticsImperial College London

South Kensington CampusLondon SW7 2AZUnited Kingdom

This thesis is submitted for the degree ofDoctor of Philosophy of Imperial College London

2016

AbstractComposite sandwich structures are widely regarded as a cost/weight-effective altern-ative to conventional composite stiffened panels and are extensively utilized for light-weight applications in various sectors, including the aeronautical, marine and trans-port industries. Nevertheless, their damage tolerance remains a critical issue.

This work aims to develop reliable analytical and numerical tools for the design ofdamage-tolerant advanced foam-cored composite sandwich structures for aerospaceapplications. It comprises of original experimental observations together with novelnumerical and analytical developments, as detailed below.

A novel analytical model for predicting the post-crushing response of crushablesandwich foam cores is presented. The calibration of the model is performed usingexperimental data obtained exclusively from standard monotonic compressive tests.Hence, the need for performing time-consuming compressive tests including multipleunloading-reloading cycles is avoided.

Subsequently, the translaminar initiation fracture toughness of a carbon-epoxyNon-Crimp Fabric (NCF) composite laminate is measured. The translaminar fracturetoughness of the UD fibre tows is related to that of the NCF laminate and the conceptof an homogenised blanket-level translaminar fracture toughness was introduced.

A multiple length/time-scale framework for the virtual testing of large compositestructures is presented. Such framework hinges upon a novel Mesh SuperpositionTechnique (MST) and a novel set of Periodic Boundary Conditions named MultiscalePeriodic Boundary Conditions (MPBCs).

The MST is used for coupling different areas of the composite structure modelledat different length-scales and whose discretizations consist of different element types.Unlike using a sudden discretization-transition approach, the use of the MST elimin-ates the undesirable stress disturbances at the interface between differently-discretizedsubdomains and, as a result, it for instance correctly captures impact-induced damagepattern at a lower computational cost.

The MPBCs apply to reduced Unit Cells (rUCs) and enable the two-scale (solid-to-shell) numerical homogenization of periodic structures, including their bending andtwisting response. The MPBCs allow to correctly simulate the mechanical responseof periodic structures using rUCs (same results as if conventional UCs were used),thus enabling a significant reduction of both modelling/meshing and analysis CPUtimes.

The developments detailed above are finally brought together in a realistic engin-eering application.

3

Keywords: Analytical modelling, Composite structures, Damage, Damage toler-ance, Finite element analysis (FEA), Foams, Large composite structures, Mesh Super-position Technique (MST), Multiple length/time-scale simulation, Multiscale mod-elling, Non-Crimp Fabrics, Numerical homogenization, PBCs, Periodic structures,Sandwich structures, Symmetries, Translaminar fracture toughness, Virtual Testing

4

Declaration of originality

The work presented hereafter is based on research carried out by the author at theDepartment of Aeronautics of Imperial College London and it is all the author’s ownwork except where otherwise acknowledged. No part of the present work has beensubmitted elsewhere for another degree or qualification.

Luigi Gigliotti

5

Copyright declaration

The copyright of this thesis rests with the author and is made available under aCreative Commons Attribution Non-Commercial No Derivatives licence. Researchersare free to copy, distribute or transmit the thesis on the condition that they attributeit, that they do not use it for commercial purposes and that they do not alter,transform or build upon it. For any reuse or redistribution, researchers must makeclear to others the licence terms of this work.

7

Acknowledgements

The work presented in this thesis was carried out under the ”Damage Toleranceof Large Sandwich Structures” project, funded by Airbus UK Ltd, and involvinga collaboration also with SWEREA SICOMP. The author would like to gratefullyacknowledge Airbus UK Ltd for sponsoring this project and making this researchpossible.

I would like to express my deepest gratitude to Prof. Silvestre Pinho for hisguidance and supervision during this project, for his constant support and invaluableexpertise, for having been a role-model in terms of dedication, passion and hard-work,for his bottomless curiosity and, above all, for being a very good friend.

I would like to gratefully acknowledge Mr. Tim Axford (who initiated the ”Dam-age Tolerance of Large Sandwich Structures” project), Mr. Tim Crundwell and Mr.Martin Gaitonde from Airbus UK Ltd for their constant support and help throughoutthis project.

I would like to thank Dr. Renaud Gutkin and Dr. Robin Olsson from SWEREASICOMP for their invaluable help and expertise, for sharing their knowledge and forthe enjoyable time spent together in several occasions throughout this project.

I would like to acknowledge Mr. Benjamin Teich (Airbus), Mr. Christian Mudra(Airbus), Mr. Rainer Tegtmeyer (Airbus) and Dr. Timothy Block (FaserinstituteBremen, now at Nordex SE) for providing materials for analysis, sharing technicaldata and many helpful discussions throughout this project.

I would like to thank the students from the Department of Aeronautics at ImperialCollege London who developed their projects in the scope of this work: Miss ChanaGoldberg, Mr. Mengfei Wang, Mr Brieux LeGall, Mr Pavel Perrotey, Miss Anne-Claire Lemoine, Mr. Vinod Raj, Mr. Ying Zhou and Mr. Cristian Ciuca.

I would like to express my gratitude to all my colleagues and friends from theDepartment of Aeronautics at Imperial College London who made the last threeyear an unforgettable and special period: Soraia, Julian, Stefano, Salvatore, Boyang,

9

Matthew, Gael, Spyros, George, Lorenza, Francesco, Gianmaria, Federico, Tobiasand many others for feedback, discussions and many enjoyable social ventures.

I would like to thank my colleague and best friend during my time at Departmentof Aeronautics at Imperial College London, Dr. Andre Wilmes, for being a true friendduring the last three years, for his passion and brilliance, for the countless late nightsin the office and for the great time we spent together during conferences.

I would like to say a wholehearted thank you to Giulia, for having been with meevery single day, for her patience and for having always believed in me more than Iever did and more than I ever will. Thank you, for everything.

Finally, I would like to thank my family and closest friends for their constant,discrete and warm presence and for having always been on my side.

10

Table of contents

Abstract 3

Declaration of originality 5

Copyright declaration 7

Acknowledgements 9

Table of contents 11

List of figures 17

List of tables 25

List of publications & dissemination 27

Acronyms & abbreviations 31

1 Introduction 35

1.1 Motivation & Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 35

1.2 Outline of this thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

2 Prediction of the post-crushing compressive response of progress-ively crushable sandwich foam cores 41

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

11

2.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.1.2 Models for the post-crushing compressive response of foam ma-terials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

2.1.3 Structure of this chapter . . . . . . . . . . . . . . . . . . . . . . 43

2.2 Model development . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

2.2.1 Monotonic compressive response . . . . . . . . . . . . . . . . . 44

2.2.2 Thickness of the uncrushed and crushed layers in crushablefoam materials . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

2.2.3 Post-crushing compressive response and residual strain . . . . . 51

2.3 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . 53

2.4 Model validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

2.4.1 Characterization of the monotonic and cyclic compressive re-sponse of the Rohacell HERO 71 foam . . . . . . . . . . . . . . 55

2.4.1.1 Material . . . . . . . . . . . . . . . . . . . . . . . . . 55

2.4.1.2 Cyclic compressive tests . . . . . . . . . . . . . . . . . 55

2.4.2 Measured foam properties . . . . . . . . . . . . . . . . . . . . . 56

2.4.3 Model predictions . . . . . . . . . . . . . . . . . . . . . . . . . 56

2.4.3.1 Thickness of the crushed layer . . . . . . . . . . . . . 56

2.4.3.2 Residual strain . . . . . . . . . . . . . . . . . . . . . . 57

2.4.3.3 Post-crushing compressive response . . . . . . . . . . 59

2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

2.5.1 Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . 61

2.5.2 Thickness of the crushed layer . . . . . . . . . . . . . . . . . . 61

2.5.3 Residual strain . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

2.5.4 Post-crushing compressive response . . . . . . . . . . . . . . . . 62

2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

12

3 Translaminar fracture toughness of NCF composites with multiaxialblankets 65

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

3.2 Experimental method . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

3.2.1 Material system . . . . . . . . . . . . . . . . . . . . . . . . . . 66

3.2.2 Specimen and layup configuration . . . . . . . . . . . . . . . . 67

3.2.3 Test method and experimental setup . . . . . . . . . . . . . . . 68

3.3 Data reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

3.3.1 NCF laminate-level translaminar fracture toughness . . . . . . 69

3.3.2 Fibre tow-level translaminar fracture toughness . . . . . . . . . 71

3.3.2.1 Relating the toughness of the individual fibre tows tothe toughness of the laminate . . . . . . . . . . . . . . 71

3.3.2.2 Relating the toughness of the off-axis fibre tows tothat of the 0◦ and 90◦ fibre tows . . . . . . . . . . . . 72

3.3.3 NCF blanket-level translaminar fracture toughness . . . . . . . 73

3.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

3.4.1 Load displacement curves . . . . . . . . . . . . . . . . . . . . . 73

3.4.2 Translaminar fracture toughness . . . . . . . . . . . . . . . . . 74

3.4.2.1 NCF laminate-level translaminar fracture toughness . 74

3.4.2.2 Fibre tow-level translaminar fracture toughness . . . 74

3.4.2.3 NCF blanket-level translaminar fracture toughness . . 76

3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4 Multiple length/time-scale simulation of localized damagein composite structures using a Mesh Superposition Technique 79

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4.1.2 3D solid elements and 2D shell elements for laminated composites 80

13

4.1.3 Multi-dimensional finite elements coupling . . . . . . . . . . . . 80

4.2 Mesh Superposition Technique . . . . . . . . . . . . . . . . . . . . . . 82

4.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

4.2.1.1 Weighting factors computation . . . . . . . . . . . . . 85

4.2.2 Finite Element implementation . . . . . . . . . . . . . . . . . . 86

4.3 Application: Multiple length/time-scale analysis of low-velocity impacton a composite plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

4.3.1 Problem description . . . . . . . . . . . . . . . . . . . . . . . . 88

4.3.2 Finite Element models . . . . . . . . . . . . . . . . . . . . . . . 89

4.3.2.1 Multiple length-scale analysis . . . . . . . . . . . . . . 91

4.3.2.2 Multiple length and time-scale analysis . . . . . . . . 92

4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

4.4.1 Damage prediction MST . . . . . . . . . . . . . . . . . . . . . . 94

4.4.2 Computational efficiency . . . . . . . . . . . . . . . . . . . . . . 94

4.4.2.1 Multiple length-scale analysis . . . . . . . . . . . . . . 94

4.4.2.2 Multiple length/time-scale analysis . . . . . . . . . . 95

4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

4.5.1 Multiple length-scale analysis . . . . . . . . . . . . . . . . . . . 96

4.5.2 Multiple length/time-scale analysis . . . . . . . . . . . . . . . . 98

4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

5 Exploiting symmetries in solid-to-shell homogenization, with applic-ation to periodic pin-reinforced sandwich structures 101

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

5.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5.2.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . 102

5.2.2 Multiscale Periodic Boundary Conditions (MPBCs) . . . . . . 104

5.2.2.1 Physical equivalence and Periodicity . . . . . . . . . . 104

14

5.2.3 Subdomain admissibility . . . . . . . . . . . . . . . . . . . . . . 105

5.2.4 Derivation of the MPBCs . . . . . . . . . . . . . . . . . . . . . 107

5.2.5 Two-scale (solid-to-shell) homogenisation of periodic structures 109

5.3 FE implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.4 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.4.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.4.1.1 Description . . . . . . . . . . . . . . . . . . . . . . . . 112

5.4.1.2 Results and discussion . . . . . . . . . . . . . . . . . . 114

5.4.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

5.4.2.1 Description . . . . . . . . . . . . . . . . . . . . . . . . 116

5.4.2.2 Results and discussion . . . . . . . . . . . . . . . . . . 117

5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6 Virtual Testing of large composite structures: a multiplelength/time-scale framework 121

6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

6.1.1 Virtual testing of large composite structures . . . . . . . . . . . 121

6.1.2 Multiscale coupling . . . . . . . . . . . . . . . . . . . . . . . . . 122

6.1.3 Solid-to-shell homogenization . . . . . . . . . . . . . . . . . . . 123

6.1.4 Structure of this chapter . . . . . . . . . . . . . . . . . . . . . . 124

6.2 Multiple length/time scale framework . . . . . . . . . . . . . . . . . . 124

6.3 Multiple length/time-scale simulation of a large aeronautical component126

6.3.1 Problem description . . . . . . . . . . . . . . . . . . . . . . . . 126

6.3.2 Implicit & Explicit FE submodels . . . . . . . . . . . . . . . . 127

6.3.3 Multiscale explicit FE submodel . . . . . . . . . . . . . . . . . 129

6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

15

7 Conclusions 133

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

7.2 Novelty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

7.3 Impact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

7.4 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

Bibliography 137

A Propagation of stress waves in non-uniform FE meshes using theMST 153

A.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

A.2 FE models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153

A.3 Stress-wave propagation analysis . . . . . . . . . . . . . . . . . . . . . 154

A.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

16

List of figures

2.1 Schematic of a foam specimen under compressive loading. In the un-deformed configuration, the specimen thickness and the cross-sectionalarea are denoted as h0 and A0, respectively; in the deformed configur-ation, under the applied displacement u, the specimen thickness andcross-sectional area are denoted, respectively, as h=h0−u and A≈A0. 44

2.2 Typical nominal stress (σ) v.s. homogenized true strain (〈ε〉) curvefor a crushable foam material under monotonic compressive loading.Three main deformation regimes can be identified: elasticity (E), pro-gressive crushing (C) and densification (D). . . . . . . . . . . . . . . 45

2.3 Piecewise constitutive law for the compressive response of crushablefoam materials: experimental measurements v.s. model predictions(black curves) within the elastic (a), progressive curshing (b) and dens-ification (c) regimes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

a Elastic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

b Progressive crushing. . . . . . . . . . . . . . . . . . . . . . . . . 46

c Densification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

2.4 Determination of the thickness of the uncrushed and crushed foammaterial layers. For any homogenized strain 〈ε〉, the correspondingcompressive stress σ is computed through Equation 2.9. At this stresslevel, the homogenized strains within the uncrushed and the crushedlayers, respectively denoted as 〈εE〉 and 〈εD〉, are computed assumingthat the uncrushed foam material behaves according to Equation 2.2and the crushed material according to Equation 2.7. The thicknesseshu and hc can therefore be computed by exploiting the equilibriumcondition, i.e. σ is constant through the entire foam specimen. . . . . 49

17

2.5 Typical evolution of the normalized (with respect to the initial speci-men thickness h0) thickness of (i) the uncrushed layer hu (red curve),(ii) the crushed layer hc (green curve) and (iii) the entire specimen h

(black dashed curve) as a function of the applied homogenized strain〈ε〉. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

2.6 Typical response of a crushable foam material subjected to a completeunloading-reloading cycle. The nominal stress and homogenized strainat unloading initiation are denoted, respectively, as σUn and 〈εUn〉. Thereloading is assumed to start immediately at σ = 0 (〈ε〉= 〈εRe〉), uponcomplete unloading. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

2.7 Modelling the post-crushing compressive response of crushable foammaterials. The thickness of the uncrushed and crushed layers atunloading initiation are denoted as h

Unu and h

Unc , respectively. For

〈ε〉 ≤ 〈εUn〉, the post-crushing compressive response (orange curve) iscomputed by imposing the equilibrium between the layers of uncrushedand crushed material, while for 〈ε〉 > 〈εUn〉, it is described by Equa-tion 2.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

2.8 Numerical implementation of the proposed model. . . . . . . . . . . . 54

2.9 Thickness of the layer of crushed material: model predictions v.s. ex-perimental measurements. The normalized thicknesses hu and hc arecalculated by analysing the discontinuous strain field within the speci-men during crushing, using the DIC technique. . . . . . . . . . . . . . 57

2.10 Residual strain upon complete unloading: model predictions (solidcurves) against experimental measurements, for the WF51, H100 andHERO 71 foams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

2.11 Percentage reduction of the error in predicting the residual strain uponcomplete unloading 〈εRe〉 using the proposed model (compared to mod-els assuming a linear [governed by the elastic modulus of the undam-aged material] post-crushing compressive response). The average val-ues of the error-reduction obtained for the three foam materials con-sidered in this work are displayed as dashed lines. . . . . . . . . . . . . 58

18

2.12 Comparison of model predictions (black solid curves) against exper-imental measurements for the post-crushing compressive response ofthe PMI Rohacell WF51 (a), PVC Divinycell H100 (b) and PMI Ro-hacell HERO 71 (c) foams (dashed black lines indicate the predictedresponse if a linear behaviour (with the undamaged elastic modulus)was assumed). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

a Rohacell WF51 foam. Experimental data from Flores-Johnsonet al. [45]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

b Divinycell H100 foam. Experimental data from Flores-Johnsonet al. [45]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

c Rohacell HERO 71 foam. . . . . . . . . . . . . . . . . . . . . . 60

3.1 CT specimens nominal dimensions (in mm) and fibre directions. . . . 67

3.2 Test set up with target points and scale. . . . . . . . . . . . . . . . . . 69

3.3 Compliance calibration curves obtained from FE and the MCCmethod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

a Layup A (see Table 3.2). . . . . . . . . . . . . . . . . . . . . . . 70

b Layup B (see Table 3.2). . . . . . . . . . . . . . . . . . . . . . . 70

3.4 Micrograph of a crack propagating across off-axis plies (angle α) in aprepreg-based composite laminate, after [79]. . . . . . . . . . . . . . . 72

3.5 Experimental load (P ) vs. displacement (d) curves for the layups in-vestigated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74



3.6 R-curves for NCF laminate-level translaminar fracture toughness. . . . 75



3.7 NCF laminate-level translaminar fracture toughnesses (initiation val-ues). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

3.8 Fibre tow-level translaminar fracture toughness (initiation value). . . . 76

3.9 NCF blanket-level translaminar fracture toughness (initiation value). . 77

19

4.1 Mesh Superposition Technique. In the reference configuration, domainΩ is decomposed into the subdomains ΩA and ΩB which overlap overthe subdomain Ωs. Dirichlet and von Neumann boundary conditionsare applied over specific portions of the boundary Γ=∂Ω. . . . . . . . 83

4.2 Distances used for the computation of the weighting factors ψA and ψB. 85

4.3 FE meshes of the superposed subdomains ΩA and ΩB. Within theMST region, the two meshes are superposed and the correspondingstiffness and mass matrices scaled, in order to satisfy the conservationof energy principle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

4.4 Impact-induced damage pattern on [0◦3/90◦

3]s laminate. X-rays analysisshows that the main failure modes are represented by a tensile matrixcracking of the distal 0◦

3 sublaminate (green) and by a two-lobe shapeddelamination at the bottom 0◦/90◦ interface (red), after [116]. . . . . . 88

4.5 Schematic of the impacted specimens with definition of the meso-scaleregion (a) and the FL model used to obtain the reference solution (b). 89

a Schematic (not in scale) with layup and dimensions of plateand impactor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

b FL model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

4.6 ST and MST models used for the 3D/3D meso/macro coupling. 3D solidelements and 3D shell elements are displayed, respectively, in blue andgreen, while the cohesive elements are shown in red. . . . . . . . . . . 91

a 3D/3D ST model. . . . . . . . . . . . . . . . . . . . . . . . . . . 91

b 3D/3D MST model. . . . . . . . . . . . . . . . . . . . . . . . . . 91

4.7 ST and MST models used for the 3D/2D meso/macro coupling. 3D solidelements and conventional 2D shell elements are displayed, respectively,in blue and green, while the cohesive elements are shown in red. . . . . 91

a 3D/2D ST model. . . . . . . . . . . . . . . . . . . . . . . . . . . 91

b 3D/2D MST model. . . . . . . . . . . . . . . . . . . . . . . . . . 91

4.8 Multiple length/time-scale models. The 3D/2D ST and MST mod-els are decomposed into an explicit and an implicit sub-model whichinteract through the interface nodes displayed in red. . . . . . . . . . . 93

a ST multi time/length-scale model. . . . . . . . . . . . . . . . . 93

20

b MST multi time/length-scale model. . . . . . . . . . . . . . . . . 93

4.9 Interlaminar damage evolution at the bottom 0◦/90◦ interface. Onthe left, the evolution of the damage variable d as function of thehorizontal distance from the impact location; the x1-coordinate cor-responds to the height where the delamination attains its maximumextension (3.375 mm for (a) and 4.875 mm for (b)). On the right, thedelamination pattern for the five FE models considered. . . . . . . . . 95

a t=2.0 ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

b t=4.0 ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

4.10 Computational efficiency of the MST. The MST models allow the use ofsmaller meso-scale areas and, thus, lower CPU times. For the proposedexample, using the MST the meso-scale area is reduced by nearly 60%,while the CPU time by nearly 23%. . . . . . . . . . . . . . . . . . . . . 96

a Normalized CPU time as function of the normalized meso-scalearea for the 3D/2D ST and MST models. The extension of themeso-scale region is progressively reduced and the data pointindicated with the red cross represents the first configurationfor which artificial damage at the discretization-transition isobserved. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

b On the left, the comparison between the smallest meso-scaleareas required to correctly capture the damage pattern, withthe ST model and the MST model; on the right, the comparisonbetween the corresponding CPU times required by the optimalST and the optimal MST model. . . . . . . . . . . . . . . . . . 96

4.11 Computational efficiency of the MST model coupled with an impli-cit/explicit co-simulation technique. . . . . . . . . . . . . . . . . . . . 97

a Number of DOFs in the explicit sub-model. . . . . . . . . . . . 97

b Number of interface nodes. . . . . . . . . . . . . . . . . . . . . 97

c CPU time reduction when using the MST. . . . . . . . . . . . . 97

5.1 Problem formulation. The macroscopic response of the periodic do-main D3 can be simulated by means of an equivalent shell model,provided the equivalent 2D constituive response of a representativethree-dimensional UC/rUC (subdomain s) is correctly determined. . . 103

21

5.2 Physical equivalence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5.3 Periodicity condition and Unit Cell (UC). . . . . . . . . . . . . . . . . 105

5.4 Geometrical relation between equivalent points. . . . . . . . . . . . . 106

5.5 The six loading cases in shell theory. . . . . . . . . . . . . . . . . . . . 110

a In-plane stretching ε1◦. . . . . . . . . . . . . . . . . . . . . . . . 110

b In-plane stretching ε2◦. . . . . . . . . . . . . . . . . . . . . . . . 110

c Shearing γ12◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

d Bending κ1◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

e Bending κ2◦. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

f Twisting κ12◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

5.6 FE implementation of the MPBCs. The DOFs of equivalent nodes onthe boundary of the UC/rUC are coupled using constraint equations;the external loading are specified through the DOFs of the master nodeM which does not belong to the mesh of the UC/rUC. . . . . . . . . 112

a FE discretization of a generic heterogeneous UC/rUC. . . . . . 112

b Equivalent nodes at the boundary of the UC/rUC. . . . . . . . 112

5.7 Schematic of the sandwich structure with unequal skins considered.The latter can be subdivided into UCs, which in turn can be furthersubdivided into rUCS. For the solid-to-shell homogenization of thisstructure we considered as rUC the subdomain denoted as s (in green).The LCSs of the rUC and of the surrounding subdomains are also shown.113

5.8 Homogenisation of a composite sandwich structure with unequal skins.The computed terms of the K matrix are compared to the analyticalvalues obtained with CLT, for different rUCs sizes (nnodes is the thenumber of nodes in the laminate plane of the rUC). . . . . . . . . . . 115

a Aii terms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

b Bii terms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

c Dii terms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

d Aij terms (i �=j). . . . . . . . . . . . . . . . . . . . . . . . . . . 115

e Bij terms (i �=j). . . . . . . . . . . . . . . . . . . . . . . . . . . 115

22

f Dij terms (i �=j). . . . . . . . . . . . . . . . . . . . . . . . . . . 115

5.9 Schematic of the UC and rUC considered for the analysis of the mech-anical response of a periodic sandwich structure with unequal skins andpin-reinforced core. The UC’s internal symmetries are exploited for thedefinition of the rUC (in green). The red dotted line A-B indicates thepath along which the membrane strains distribution, computed withthe UC and the rUC models, have been compared. . . . . . . . . . . . 117

5.10 Membrane strains distributions in the UC and rUC models along thepath A-B shown in Figure 5.9 under twisting load. Since the discretiza-tions of the UC and rUC models are the same, the strains distributionscomputed with the UC and the rUC models are, as expected, identical. 118

a Membrane direct strain e11. . . . . . . . . . . . . . . . . . . . . 118

b Membrane direct strain e22. . . . . . . . . . . . . . . . . . . . . 118

c Membrane shear strain e12. . . . . . . . . . . . . . . . . . . . . 118

5.11 The MPCs enable the use of rUCs for the analysis of the mechanicalresponse of periodic structures. This translates into a reduction of thetotal CPU time required (a), as a result of the reduction of both themodelling/meshing CPU time tM needed to create the FE model andapply the suitable MPBCs (b) and the analysis CPU time tA necessaryto run the FE analysis (c). . . . . . . . . . . . . . . . . . . . . . . . . 119

a Total CPU time. . . . . . . . . . . . . . . . . . . . . . . . . . 119

b Modelling/meshing CPU time. . . . . . . . . . . . . . . . . . . 119

c Analysis CPU time. . . . . . . . . . . . . . . . . . . . . . . . . 119

6.1 Multiple length/time-scales framework for the virtual testing of largecomposite components. Such framework consists of a Mesh Super-position Technique for coupling differently-discretized subdomains (A)and on the exploitation of symmetries in the solid-to-shell numericalhomogenization of periodic structures. . . . . . . . . . . . . . . . . . . 125

6.2 Schematic (not in scale) of the helicopter rotor blade considered in thisstudy; dimensions are in mm. . . . . . . . . . . . . . . . . . . . . . . . 127

23

6.3 Implicit/explicit submodels of a multiple length/time-scale FE modelof an helicopter rotor blade. The area surrounding the impact locationis analysed using an explicit solver, to exploit its capabilities to betterhandle contact interactions and complex failures modes; the responseof the remaining of the structure (the largest part) is more efficientlysimulated using an implicit solver. . . . . . . . . . . . . . . . . . . . . 128

6.4 Multiscale explicit FE submodel. Different areas of the structure canbe modelled at different length-scales and their coupling be performedusing the MST. The mechanical response of the structure can be cor-rectly captured at all the length-scales of interest (both in terms ofgeometrical details and failure modes), while keeping the computa-tional cost of the analysis to a minimum. . . . . . . . . . . . . . . . . 130

A.1 Three different FE models of the infinite bar: Fully Solid model, Sud-den Transition model and MST model. . . . . . . . . . . . . . . . . . 154

A.2 Applied displacement and boundary conditions for the stress-wavepropagation study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

A.3 Time and space evolution of the applied displacement u(z, t). Theresulting applied displacement is a finite discrete impulse with variableamplitude along the through-the-thickness direction. . . . . . . . . . 155

a Normalized amplitude of the applied displacement as a functionof time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

b Amplitude of the applied displacement as a function of the z-coordinate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

A.4 Time history of the internal energy of Region 1 (a) and Region 2 (b)associated with the propagation of the travelling stress waves along theinfinite bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

a Region 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

b Region 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

24

List of tables

2.1 Measured properties and model calibration parameters for the RohacellWF51 [45], Divinycell H100 [45] and Rohacell HERO 71 foams. . . . . 56

3.1 Nominal membrane properties of the fibre tows in the triaxial NCFblanket [80,81]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

3.2 Layups investigated. The nominal thickness of the laminates is indic-ated as tLam and the 0◦ fibre tows are aligned with the direction of theapplied load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

3.3 Numerical fitting parameters used in the MCC method (units system:kN; mm). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.1 Elastic properties of the plies, cohesive properties of the interfaces, anddensity used by Aymerich et al. [115]. . . . . . . . . . . . . . . . . . . 92

5.1 Elastic properties, density and thickness of the NCF triaxial weave[155,156]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

5.2 Elastic properties and density of the HERO G3 150 foam [154]. . . . . 114

5.3 Sandwich layers’ thicknesses and facesheets layups. . . . . . . . . . . . 114

5.4 Geometrical parameters defining the rUC of the periodic sandwichstructure with unequal skins and pin-reinforced core. . . . . . . . . . . 116

25

List of publications &dissemination

Parts of the work presented in this thesis have been disseminated through writtenpublications, oral presentations as well as which are listed below as of January 2016.

Peer–reviewed journal publications

[p1] L. Gigliotti and S. T. Pinho,“Translaminar fracture toughness of NCF com-posites with multiaxial blankets”, Materials & Design, In Press, 2016.(see Chapter 3)

[p2] L. Gigliotti and S. T. Pinho,“Virtual testing of large composite struc-tures: a multiple length/time-scale framework”, Journal of Multiscale Modelling,vol. 6(3), 2016.(see Chapter 6)

[p3] L. Gigliotti and S. T. Pinho,“Prediction of the post-crushing compressive re-sponse of progressively crushable sandwich foam cores”, Composites: Part A,vol. 80, pp. 148–158, 2016.(see Chapter 2)

[p4] L. Gigliotti and S. T. Pinho,“Exploiting symmetries in solid-to-shell homogen-ization, with application to periodic pin-reinforced sandwich structures”, Com-posite Structures, vol. 132, pp. 995–1005, 2015.(see Chapter 5)

[p5] L. Gigliotti and S. T. Pinho,“Multiple length/time-scale simulation of localizeddamage in composite structures using a Mesh Superposition Technique”, Com-posite Structures, vol. 121, pp. 395–405, 2015.(see Chapter 4)

27

Refereed conference publications

[c1] L. Gigliotti and S. T. Pinho,“Damage tolerance of sandwich foam cores: exper-imental characterization and numerical modelling”, American Society for Com-posites 30th Technical Conference, East Lansing (MI), United States, 28th–30th

September 2015, http://www.egr.msu.edu/asc2015/.

[c2] L. Gigliotti and S. T. Pinho,“Enabling faster structural design: Efficientmultiscale simulation of large composite structures”, ICCM20, 20th InternationalConference on Composite Materials, Copenhagen, Denmark, 19th–24th July 2015,http://www.iccm20.org.

[c3] L. Gigliotti and S. T. Pinho,“A Mesh Superposition Technique for the simulationof the mechanical response of composite materials at multiple length and time-scales”, American Society for Composites 29th Technical Conference - 16th US-Japan Conference on Composite Materials - ASTM-D30 Meeting, San Diego(CA), United States, 8th–10th September 2014, http://www.asc-composites.org/asc2014.

[c4] L. Gigliotti and S. T. Pinho,“On transition regions for the simulation of themechanical response of composite materials at multiple length-scales”, ECCM16,16th European Conference on Composite Materials, Seville, Spain, 22nd–26th June2014, www.eccm16.org.

Refereed conference abstracts

[a1] L. Gigliotti and S. T. Pinho, “Post-crushing response of sandwich foam cores:model identification and virtual testing”, 5th ECCOMAS Thematic Confer-ence on the Mechanical Response of Composites: COMPOSITES 2015, Bris-tol, United Kingdom, 7th–9th September 2015, http://www.bristol.ac.uk/composites/composites2015/.

[a2] S. T. Pinho and L. Gigliotti, “Multiple length/time-scale analysis of largecomposite structures”, 5th ECCOMAS Thematic Conference on the MechanicalResponse of Composites: COMPOSITES 2015, Bristol, United Kingdom, 7th–9th

September 2015, http://www.bristol.ac.uk/composites/composites2015/.

[a3] S. T. Pinho and L. Gigliotti, “Efficient multiple length/time-scale simulation oflarge composite structures”, 18th International Conference on Composite Struc-tures, Lisbon, Portugal, 15th–18th June 2015.

28

[a4] L. Gigliotti and S. T. Pinho, “Energy-absorption of lightweight hybrid-construction sandwich materials: model identification and virtual testing”, 26th

Annual International SICOMP Conference, Goteborg, Sweden, 1st–2nd June2015, http://www.swerea.se/en/sicomp-conference.

[a5] L. Gigliotti and S. T. Pinho, “On transition regions for the simulation ofthe mechanical response of composite materials at multiple length scales”, IW-CMM23, 23rd International Workshop on Computational Mechanics of Materials,Singapore, Singapore, 2nd–4th October 2013.

Invited seminar presentations[s1] S. T. Pinho and L. Gigliotti,“Multiscale analysis of large composite structures”,

N8 HPC Network Event - Multiscale Computational Mechanics at the Universityof Sheffield, Sheffield, United Kingdom, 30th October 2015, http://n8hpc.org.uk/n8-hpc-network-event-multiscale-computational-mechanics/.

[s2] L. Gigliotti and S. T. Pinho,“Multiple length/time scale simulation oflarge composite structures”, NASA Langley Research Center, Hampton (VA),United States, 25th September 2015, http://www.nianet.org/education/continuing_education/2015-seminar/9-25-15-gigliotti/.

29

Acronyms & abbreviations

ASTM . . . . . . . . . . . . . . . American Society for Testing and Materials International

B31 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2-noded, first-order 3D beam element

(Abaqus elements library)

CFL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Courant-Friedrichs-Lewy

CFRP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Carbon Fiber Reinforced Polymer

CIN3D8 . . . . . 8-noded, first-order, one-way infinite 3D continuum infinite element


CLT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Classical Lamination Theory

COH3D8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8-noded, first-order, cohesive 3D element


CPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Central Processing Unit

CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Compact Tension

CSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Co-Simulation Engine

C3D8R . . . . . . . . . 8-noded, first-order, reduced integration 3D continuum element


DOF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Degree Of Freedom

DCB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Double Cantilever Beam

EEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Embedded Element Method

31

ELST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Equivalent Single Layer Theories

ENF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . End Notched Flexure

EU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . European Union

FALCOM . . . . . . . . . . . . Failure, performance and processing prediction for enhanced

design with non-crimp fabric composites

FE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Finite Element

FETI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Finite Element Tearing and Interconnect

FL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fully Local

FS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Fully Solid

GC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Gravouil and Combescure

LaRC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Langley Research Center

LCS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Local Coordinate System

LSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Level Set Method

MCC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Modified Compliance Calibration

MPBCs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multiscale Periodic Boundary Conditions

MPC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Multi-Point Constraint

MST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mesh Superposition Technique

NASA . . . . . . . . . . . . . . . . . . . . . . . . .National Aeronautics and Space Administration

NCF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Non-Crimp Fabric

PBCs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Periodic Boundary Conditions

PMI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Polymethacrylimide

PTFE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Polytetrafluoroethylene

PVC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . PolyVinyl Xhloride

PU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Partition of Unity

rUC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . reduced Unit Cell

RAM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Random Access Memory

R3D4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4-noded, first-order, 3D rigid element

32


SC8R . . . . 8-noded, first-order, reduced integration 3D continuum shell element


ST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sudden Transition

S4R . . . . . . . . . . . . . . . . . . . 4-noded, first-order, reduced integration shell element


S8R5 . . . . . . . . . . . . . . 8-noded, first-order, reduced integration thin shell element


TECABS . . . . . . . . . . . . . . . . . . . . . . . .Technologies for Carbon fibre reinforced modular

Automotive Body Structures

UC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Unit Cell

UD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . UniDirectional

VCCT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Virtual Crack Closure Technique

vs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . versus

1D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . One Dimensional

2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Two Dimensional

3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Three Dimensional

33

Notation & nomenclature

Lower case roman letters

a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .crack length

ai . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . polynomial fitting coefficient of order i

c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . dilatational wave speed

d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . relative displacement, signed distance

di . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . signed distance from surface Si, i = A, B

e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Euler’s number, engineering strain

eij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . membrane strains, i = 1, 2

eji . . . . . . . . . . . . . . . . . .j-th element of the discretization of subdomain Ωi, i = A, B

ei . . . . . . . . . . . . . . . . . . . . . orthogonal basis of the LCS of subdomain s, i, j = 1, 2, 3

ei . . . . . . . . . . . . . . . . . . . . . orthogonal basis of the LCS of subdomain s, i, j = 1, 2, 3

fmax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . minimum frequency

fmin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .maximum frequency

f〈ε〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . applied homogenized true strain history function

gλi

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . fitting function with coefficients λi

h0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . specimen’s thickness (undeformed configuration)

h . . . . specimen’s thickness (deformed configuration), deformable body thickness

h . . . . . . . . . . . . . . . . . . . . normalized specimen’s thickness (deformed configuration)

hc . . . . . . . . . . . . . . . . . . thickness of the fully-crushed layer (deformed configuration)

hc . . . . . . normalized thickness of the fully-crushed layer (deformed configuration)

35

hUnc . . . . . . . . . . . . . . . . . . . . . . . . . thickness of the crushed layer at unloading initiation

(deformed configuration)

hUnc . . . . . . . . . . . . . .normalized thickness of the crushed layer at unloading initiation

hw . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness of the triaxial NCF weave

htf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness of the top facesheet

hbf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness of the bottom facesheet

hc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness of the foam core layer


hcohe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .cohesive elements thickness

hu . . . . . . . . . . . . . . . . . . . . thickness of the uncrushed layer (deformed configuration)

hUnu . . . . . . . . . . . . . . . . . . . . . . . thickness of the uncrushed layer at unloading initiation


hu . . . . . . . . . normalized thickness of the uncrushed layer (deformed configuration)

hUnu . . . . . . . . . . . normalized thickness of the uncrushed layer at unloading initiation


i . . . . . . . . . . . . . . . . . relative signed distance from surface SI , (i, I) = (a,A), (b,B)

k . . . . . . . . . . . . . . . . . . . . . . . . . . . constant of proportionality, numerical loading time

kmax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . maximum value of numerical loading time

kN , kT , kS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . cohesive zone stiffnesses

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . characteristic length of a spatial grid

e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . finite element in-plane size

i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane dimensions of the UC/rUC, i = 1, 2

e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . embedded element size

h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . hosting element size

n . . . . . . . . . . . . . . . . . . . . order of polynomial fitting, number of physical properties

nnodes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . number of nodes in the laminate plane

n〈ε〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .dimensions of vector 〈ε〉s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . subdomain of D3

s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . subdomain of D3, equivalent to s

36

t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness, time

tLam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . laminate’s thickness

tαK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness of the fibre tows oriented at an angle α

within specimens with generic layup K

tK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . thickness of specimens with generic layup K

t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . total CPU time

tA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . analysis CPU time

tM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . modelling/meshing CPU time

tUC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . total CPU time using a UC

tUC

A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . analysis CPU time using a UC

tUC

M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . modelling/meshing CPU time using a UC

t[ST]min . . . . . . . . . . . . . . . . . minimum (normalized) CPU time needed with the ST model

t[MST]min . . . . . . . . . . . . . . . . minimum (normalized) CPU time needed with the MST model

u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . applied displacement

u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . applied displacement

u . . . . . . . . . . . . . . . . . . . . . . . . . . . . displacement field, in-plane projection of vector u

u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . three-dimensional displacement vector

x . . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane 1st coordinate, independent scalar variable

xAi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . coordinates of point A, i = 1, 2, 3

xExp . . . . . . . . . . . . . . . . . . . . . . . . . independent variable vector (experimental measures)

y . . . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane 2nd coordinate, dependent scalar variable

yExp . . . . . . . . . . . . . . . . . . . . . . . . . . . dependent variable vector (experimental measures)

z . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . through-the-thickness coordinate

w . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . out-of-plane displacement

Δa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . crack growth

Δt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . time increment

ΔtImp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . implicit sub-model time increment

Δtstable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . stable time increment

ΔtXpl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . explicit sub-model time increment

Δtmin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . minimum CPU time reduction

37

Upper case roman letters

A . . . . . . . . . . . . . . . . . . . . . specimen’s cross-sectional area (deformed configuration),

generic point in R3

Aij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . extensional stiffnesses

A0 . . . . . . . . . . . . . . . . . . . specimen’s cross-sectional area (undeformed configuration)

Ameso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane area of the meso-scale region

ArUC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane area of the rUC

Atot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . total in-plane area

A[ST]meso . . . . . . . . . . . . minimum (normalized) meso-scale area needed with the ST model

A[MST]meso . . . . . . . . . . . minimum (normalized) meso-scale area needed with the MST model

Bij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . bending-extension coupling stiffnesses

B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . body forces vector

Bi . . . . . . . . . . . . . . . . . . . . . . . . . . . . body forces vector within subdomain Ωi, i = A, B

Bs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . body forces vector within transition subdomain

Bej

i. . . . . . . . . . . . . . . . . . shape functions derivatives matrices of element ej

i , i = A, B

B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . three-dimensional deformable body

B3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . three-dimensional deformable body

C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . compliance

C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . compliance matrix

C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Courant number

Cmax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .maximum Courant number

C0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . set of all continuous function

Ci . . . . . . . . . . . . fourth-order material elasticity tensor associated to subdomain Ωi

Dij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . bending stiffnesses

D2 . . . . . . . . . . . . . . . . . . . . . .domain occupied by a two-dimensional deformable body

D3 . . . . . . . . . . . . . . . . . . . . domain occupied by a three-dimensional deformable body

E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . elastic modulus (undamaged material)

E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Green-Lagrange strain tensor

38

Ef . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . foam material Young’s modulus

Ei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Young’s modulus in direction i = 1, 2, 3

EL . . . . . . . . . . . . . . . . . . . . . . tangent modulus of the fully-crushed material at σD = 0

EL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . tangent modulus at the onset of densification

Epin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . pins Young’s modulus (fibre direction)

F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . deformation gradient

Gij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . shear modulus in direction i, j = 1, 2, 3

GIc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . interlaminar fracture toughness (Mode I)

GIIc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . interlaminar fracture toughness (Mode II)

GIIIc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . interlaminar fracture toughness (Mode III)

GαIc . . . . . . . . . . . . . . . . . . . . . . . . . translaminar fracture toughness of off-axis fibre-tows

(at an angle α with respect to the 0◦ fibre-tows)

GKIc . . . . . . . . . . . . . . . . . . laminate-level translaminar fracture toughness of a laminate

with generic layup K

GLamIc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . laminate-level translaminar fracture toughness

GNCFIc . . . . . . . . . . . . . . . . . . . . . . . . . . . NCF blanket-level translaminar fracture toughness

I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . identity matrix

I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 × 3 identity matrix

K . . . . . . . . . . . . finite element stiffness matrix, shell stiffness matrix, ABD matrix

KTr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . transition finite element stiffness matrix

Ksej

i. . . . . . . . . . . stiffness matrix of the element ej

i within the MST region, i = A, B

L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . panel’s length (y-direction)

Lmeso . . . . . . . . . . . . . . . . . . . . . . . . . in-plane length (y-direction) of the meso-scale region

Mi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . bending moment per unit/length, i = 1, 2

M12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . twisting moment per unit/length

Msej

i. . . . . . . . . . . . . . mass matrix of the element ej

i within the MST region, i = A, B

N, T, S . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .cohesive zone strengths

Ni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .direct force per unit/length, i = 1, 2

N12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . shear force per unit/length

39

Nej

A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .shape functions matrices of element ej

i , i = A, B

O . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . origin of the LCS of subdomain s

O . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . origin of the LCS of subdomain s

P . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . applied load, generic point in R3

Pc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . critical applied load

S . . . . . . . . . . restraint matrix, vector of resultant forces/moments per unit-length

Rn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . real coordinate space

S . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Piola-Kirchhoff stress tensor

Si . . . . . . . . . . . . . . . . . . Piola-Kirchhoff stress tensor within subdomain Ωi, i = A, B

Ss . . . . . . . . . . . . . . . . . . . . .Piola-Kirchhoff stress tensor within transition subdomain

Si . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . surface, i = A, B

T . . . . . . . . . Piola traction vector, in-plane projection of transformation matrix T

T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . three-dimensional transformation matrix

U . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . internal energy

X . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . position vector (reference configuration)

Xc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . compressive (crushing) strength

XL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lock-up/densification strength

W . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . panel’s width (x-direction)

Wext . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . external work

Wint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . internal work

Wkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . inertial work

Wmeso . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane width (x-direction) of the meso-scale region

ΔAmeso . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . minimum meso-scale area reduction

Lower case greek letters

α . . . . . . . . . . . . . . . . . . ply orientation angle, ratio embedded/hosting element sizes

α, β, χ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .MCC fitting parameters

β . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . in-plane stitching angle

40

εi◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .membrane direct strains

γ12◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . membrane shear strain

κ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fourier series terms

κ1◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . bending curvatures

κ12◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . twisting curvature

εc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . homogenized true strain at crushing initiation

〈ε〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .homogenized true strain

εL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . lock-up/densification true strain

〈εE〉 . . . . . . . . . . . . . . . . . . . . . . . . . . homogenized true strain within the uncrushed layer

〈εD〉 . . . . . . . . . . . . . . . . . . . . . . . .homogenized true strain within the fully-crushed layer

εL . . . . . . . . . . . . . . homogenized true strain in the fully-crushed material at σD = 0

〈εRe〉 . . . . . . . . . . . . . . . . . . . homogenized residual true strain upon complete unloading

〈εUn〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .homogenized true strain at unloading initiation

〈ε〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .applied homogenized true strain discrete vector

〈ε〉Exp . . . . . . . . . . . . . .homogenized true strain discrete vector - experimental measures⟨ε

Exp

Re

⟩. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . experimentally measured residual strain⟨

εNum

Re

⟩. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . computed residual strain

γ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . loading reversal factor

λi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . fitting function coefficients

λ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . fitting function coefficients vector

νij . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Poisson’s ratio in direction i, j = 1, 2, 3

νf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . foam material Poisson’s ratio

π . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . pi

ξM . . error in predicting the residual strain 〈εRe〉 using the model proposed in § 2

ξL . . . . . . . . . . error in predicting the residual strain 〈εRe〉 assuming a linear elastic

post-crushing compressive response

ξ◦ . . . . . . . . . . . . . . . . . . . . resultant mid-surface membrane strains/curvatures vector

ρ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . material mass density

ρ∗ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . foam material nominal relative density

41

ρD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . fully-crushed material density

ρi . . . . . . . . . . . . . . . . . . . . . . . . material mass density within subdomain Ωi, i = A, B

ρf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . foam material density

ρL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . density of the fully-crushed material at σD = 0

ρPMI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . polymethacrylimide density

ρs . . . . . . . . . . . . . . . . . . . . . . . . . . . material mass density within transition subdomain

σ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . nominal stress

σC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . nominal stress - progressive crushing regime

σD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .nominal stress - densification regime

σE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .nominal stress - elastic regime

σExpC . . . . . . . . . . . . . . . . . . . . . . . . . . nominal stress discrete vector (progressive crushing)

- experimental measures

σExpD . . . . . . . . . . . . . . . . . . . . . . . . nominal stress discrete vector (densification crushing)

- experimental measures

σUn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . nominal stress at unloading initiation

σ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . nominal stress discrete vector

ψi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . weighting factors, i = A, B

ϕ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . through-the-thickness stitching angle

φ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . diameter of pin-reinforcements

ω . . . . . . . . . . . . . . . . . . . . domain occupied by a three-dimensional deformable body


Δξ . . . . . . . . . . . . . . . . . . reduction of the error in predicting the residual strain 〈εRe〉

Upper case greek letters

Ei . . . . . . . . . . . . . . . . . . . . . set of elements with support in subdomains Ωi, i = A, B

EB . . . . . . . . . . . . . . . . . . . . . set of elements with support in subdomains Ωi, i = A, B

Esi . . . . . . . set of elements with support in the superposed portions of subdomains

Ωi, i = A, B

42

Γ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . boundary of domain Ω

Γt . . . . . . . . . . . . . subset of Γ where von Neumann boundary conditions are applied

Γu . . . . . . . . . . . . . . . . . . subset of Γ where Dirichlet boundary conditions are applied

Λ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Least Square Method operator

Πj . . . . . . . . . . . . . . . . . . . . . . . spatial distribution of physical properties, j∈{1, . . . , n}Ω . . . . . . . . . . . . . . . . . . . . domain occupied by a three-dimensional deformable body

(reference configuration)

Ωi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . subdomains, i = A, B

Ωs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .transition subdomain, MST subdomain

(Ωi) . . . . . . . . . . . . . . . . . . . . . . . . . . .non overlapping portion of subdomain Ωi, i = A, B

Ωej

i. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . integration domain of element ej

i

Mathematical operators, symbols and accents

= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is equal to

≈ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is approximately equal to

≡ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is equivalent/congruent to∧= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is physically equivalent to

> . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is strictly greater than

< . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is strictly less than

≥ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is greater than or equal to

≤ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is less than or equal to

�= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is not equal to

± . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .plus/minus

ln(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .natural logarithm, i.e. ln (e•) = •log(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . common logarithm, i.e. log (10•) = •exp(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .exponential function, i.e. exp(•) = e•

sin(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . sine operator

cos(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .cosine operator

43

√(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . square root operator

(•)T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . tensor transpose on first two indices

(•)t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . tensor transpose on first two indices

s〈•〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . volume average operator over subdomain s

s• . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .first-order fluctuation term over subdomain s

s˜• . . . . . . . . . . . . higher-order (2nd and higher) fluctuation terms over subdomain s

〈 • 〉 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . averaged quantity

˙(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . first time derivative operator

¨(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . second time derivative operator

d(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . infinitesimal operatord(a)d(b)

. . . . . . derivative of vector/tensor quantity a relative to all dependencies in the

vector/tensor quantity b, adhering to the numerator layout convention∫(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . integral operator∑(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .summation operator

Δ(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . finite difference operator

∂(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . boundary operator

δ(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .virtual variation

∇(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Del operator

(•)◦ . . . . . . . . . . . . . . . . . . . . . . . . . . . . quantity evaluated at the midsurface of a domain

| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . such that

⇒ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . implies

⇐ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . if

∀ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . for all

∃ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . exists

∧ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . and

‖ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is parallel to

[ a , b ] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . closed interval / domain with endpoints a and b

{ a1, a2, · · · , an } . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . unordered n-tuple

· . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vector dot product

44

: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . double contraction operator

⊗ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . tensor product operator

‖•‖ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . norm of a tensor

|•| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .absolute value

∈ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . member of/in

/∈ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . not member of/in

⊂ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is a subset of

⊆ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is a subset of/coincident with

∩ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . intersection of two sets

∅ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . empty set

\ . . . . . . . . . . . . . . . . . . . . . . . . . . . . relative complement, i.e. B \ A = {x ∈ B | x /∈ A}⋃(•) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . union operator (set theory)

45

Chapter 1

Introduction

1.1 Motivation & Objectives

Composite sandwich structures consist of two composite facesheets separated by alow-density core designed to sustain the transverse shear and through-the-thicknessloads. The resulting structure combines the favourable in-plane stiffness and strengthof the composite facesheets with superior bending stiffness-to-weight ratios.

Therefore, composite sandwich structures are widely regarded as a cost/weight-effective alternative to conventional composite stiffened panels. Nowadays, compos-ite sandwich structures are extensively used for lightweight applications in varioussectors, including the aeronautical [1,2], naval [3] and transport industries [4,5]. Un-fortunately, damage tolerance of composite sandwich structures remains a criticalissue; numerous studies demonstrate that, owing to the low bending stiffness of thefacesheets, composite sandwich structures are highly susceptible to localized through-the-thickness loads such as those occurring during low-velocity impact [6–9]. This isparamount for foam-cored sandwich structures which, primarily for this reason, havenot yet found common application in the aeronautical industry, where honeycombcores are still preferred.

This work aims to develop reliable analytical and numerical tools for the design ofdamage-tolerant advanced foam-cored composite sandwich structures for aerospaceapplications. To this purpose, the focus of this work has been dedicated to the topicsdetailed below.

Post-crushing behaviour of foam materialsIn foam-cored sandwich structures, the local damage induced by a low-velocityimpact event consists of a region of crushed core accompanied by a residual dent

47

Chapter 1

in the impacted facesheet. Such residual dent results from the combinationof: (i) the residual local stress field underneath the load introduction pointand (ii) the extent of damage in both the impacted facesheet and the crushedcore [10, 11]. In addition, the residual dent constitutes a possible source offurther damage-growth upon subsequent reloading, and it is shown to severelyreduce the residual local stiffness and strength of the sandwich structure.

Thus, for damage-tolerant design of foam-cored composite sandwich structures,it is of paramount importance to accurately predict:

(i) the residual after-crushing strain in the foam core, as it contributes to thefinal depth of the residual after-impact dent;

(ii) the post-crushing compressive behaviour of the foam core, as it contrib-utes to the residual local stiffness and strength of the impacted sandwichstructure.

Translaminar fracture toughness of carbon NCF compositesAs a result of the increasing share of composite materials in sectors wherecycle times and manufacturing costs significantly impact the final product cost,the development of inexpensive and automated production methods is crucial[12–15].

The need for cost-effective alternatives to conventional prepreg-based compos-ites led for instance to the development of Non-Crimp Fabric (NCF) composites[16–18]. When compared to their prepreg-based counterpart, NCF compositesoffer higher deposition rates, reduced labour time, higher degree of tailorabilityand improved impact properties [19,20]. Therefore, NCF composites are widelyregarded as one of the most promising technologies for both aerospace [21, 22]and automotive [23, 24] structural composites, including composite sandwichstructures. The growing industrial interest towards NCF composites led to twoEU-funded research projects: FALCOM [25,26] and TECABS [27,28], respect-ively for aerospace and automotive applications.

Failure in NCF composites can be predicted using physically-based failure cri-teria [25, 29, 30]; physically-based failure criteria require, as input data, homo-genised (fibre tow- and blanket-level) properties which can be measured mostlyfrom standard tests. Specifically, translaminar fracture toughnesses are para-mount for the damage-tolerant design of composite structures.

The main challenges related to this subject addressed in this thesis are:

(i) to measure experimentally the translaminar fracture toughness of NCFcomposites;

(ii) to analytically relate the translaminar fracture toughness of off-axis fibretows/NCF blankets to that of axially-loaded fibre tows/NCF blankets.

48

Introduction

Multiscale virtual testing of large composite structuresTo address complex industrial structural challenges, such as the design of thecomposite central wing box of the Airbus A380 [31], virtual testing methodshave been exploited [32]. The extensive use of virtual testing based on nonlin-ear FE analyses is envisioned to be a key-aspect towards an increased confid-ence in the real-scale and expensive structural tests required for certification;furthermore, virtual testing provides useful insight into the likelihood, causesand consequences of structural failure [33–35]. However, to be fully establishedin structural design and certification, virtual testing methods need to be valid-ated against all level of structural testing, from the coupon-level (e.g. materialspecimens) to the system-level (e.g. wing or fuselage) [36].

Within this framework, the virtual testing of large-scale composite structures forindustrial applications entails significant challenges; these are primarily ascrib-able to the inherently multiscale nature of composite materials and, as a result,to their highly complex failure modes [37]. For the efficient structural design oflarge composite components, their virtual testing often requires that differentparts of the structure are modelled at multiple length- and time-scales, eventu-ally even using different physics.

Hence, it is crucial to develop:

(i) suitable techniques for coupling areas of the structure modelled at differentlength-scales, i.e. discretized using different finite element types;

(ii) numerical methods to efficiently compute equivalent homogenized proper-ties to be used in both 2D FE models and in the lower-scale subdomainsof multiscale FE models of large composite components.

1.2 Outline of this thesis

This thesis is divided into several chapters which address the challenges expressed in§ 1.1 one by one, and then finally brings together the developments made. Because thelatter involve several fields (e.g. numerical and experimental), the literature relevantfor each contribution is reviewed inside the respective chapter, thereby, putting thedevelopments made directly and effectively in their context.

Chapter 2 Prediction of the post-crushing compressive response of progressivelycrushable sandwich foam coresIn this chapter, a novel analytical model for predicting the post-crushingcompressive response of progressively crushable sandwich foam cores is

49

Chapter 1

presented. The calibration of the model is performed using experimentalmeasurements obtained exclusively from standard monotonic compress-ive tests. Therefore, the need for performing time-consuming compressivetests including multiple unloading-reloading cycles is avoided. Model pre-dictions have been validated against experimental measurements availablefor three different foam materials. The model is shown to accurately pre-dict the thickness of the crushed material layer during progressive crushingand the residual after-crushing strain (the latter with a maximum error of12.1%). The proposed model is capable of predicting the residual after-crushing strain with a significantly smaller error (error-reduction over 56%)than existing models, whose calibrations require the same experimentalmeasurements as the present model. The results presented in this chapterdemonstrate the relevance of the proposed model for a damage-tolerantdesign of foam-cored composite sandwich structures.

Chapter 3 Translaminar fracture toughness of NCF composites with multiaxialblanketsIn this chapter, the translaminar initiation fracture toughness of a carbon-epoxy Non-Crimp Fabric (NCF) composite laminate was measured usinga Compact Tension (CT) test. The translaminar fracture toughness of theindividual UD fibre tows was related to that of the NCF laminate and theconcept of an homogenised blanket-level translaminar fracture toughnesswas introduced. Using an approach developed for UD-ply prepreg compos-ites, it is demonstrated that the translaminar fracture toughness of off-axisfibre tows/NCF blankets can be analytically related to that of axially-loaded fibre tows/NCF blankets with a difference between experimentally-measured and predicted values lower than 5%.

Chapter 4 Multiple length/time-scale simulation of localized damagein composite structures using a Mesh Superposition TechniqueIn this chapter, a Mesh Superposition Technique (MST) for the progressivetransition between differently-discretized subdomains is proposed and thekey-aspects of its implementation in an FE code are presented. The inter-faces between these subdomains are replaced by transition regions wherethe corresponding meshes are superposed. The MST is applied to the mul-tiple length/time-scale analysis of a low-velocity impact of a projectile ona composite plate. Unlike when using a sudden discretization-transitionapproach, the use of the MST eliminates the undesirable stress disturb-ances at the interface between differently-discretized subdomains and, as a

50

Introduction

result, it correctly captures the impact-induced damage pattern at a lowercomputational cost. Finally, the MST is coupled with an implicit/explicitco-simulation technique for a multiple time/length-scale analysis. The res-ults indicate that, if the length-scale transition is performed using the pro-posed MST instead of a sudden discretization-transition, the CPU timecan be nearly halved.

Chapter 5 Exploiting symmetries in solid-to-shell homogenization, with applicationto periodic pin-reinforced sandwich structuresIn this chapter, a novel set of Periodic Boundary Conditions namedMultiscale Periodic Boundary Conditions (MPBCs) that apply to reduced

Unit Cells (rUCs) and enable the two-scale (solid-to-shell) numerical ho-mogenization of periodic structures, including their bending and twistingresponse, is presented and implemented in an FE code. Reduced UnitCells are domains smaller than the Unit Cells (UCs), obtained by exploit-ing the internal symmetries of the UCs. When applied to the solid-to-shellhomogenization of a sandwich structure with unequal skins, the MPBCsenable the computation of all terms of the fully-populated ABD mat-rix with negligible error, of the order of machine precision. Furthermore,using the MPBCs it is possible to correctly simulate the mechanical re-sponse of periodic structures using rUCs (retrieving the same results as ifconventional UCs were used), thus enabling a significant reduction of bothmodelling/meshing and analysis CPU times. The results of these ana-lyses demonstrate the relevance of the proposed approach for an efficientmultiscale modelling of periodic materials and structures.

Chapter 6 Virtual Testing of large composite structures: a multiple length/time-scaleframeworkThis chapter illustrates a multiple length/time-scale framework for the vir-tual testing of large composite structures. Such framework hinges uponthe Mesh Superposition Technique (MST) presented in Chapter 4 (usedfor coupling areas of the structure modelled at different length-scales) andupon the solid-to-shell numerical homogenization which exploits the in-ternal symmetries of Unit Cells (UCs), presented in Chapter 5. The rel-evance and key-aspects of the multiple length/time-scale framework aredemonstrated through the analysis of a real-sized aeronautical compositecomponent.

51

Chapter 1

The overall conclusions of this work and suggestions for further work are sum-marised in chapter 7.

52

Chapter 2

Prediction of the post-crushingcompressive response ofprogressively crushable sandwichfoam cores

2.1 Introduction

2.1.1 Motivation

Composite sandwich structures consist of two composite facesheets separated by alow-density core designed to sustain the transverse shear and through-the-thicknessloads; nowadays, composite sandwich structures are extensively used for lightweightapplications in various sectors, including the aeronautical [1,2], naval [3] and transportindustries [4, 5].

Nevertheless, numerous studies demonstrate that, owing to the low bending stiff-ness of the facesheets, sandwich structures are highly susceptible to localized through-the-thickness loads such as those occurring during low-velocity impact [6–9]. In foam-cored sandwich structures, the local damage caused by a low-velocity impact eventconsists of a region of crushed core accompanied by a residual dent in the impactedfacesheet [38, 39]. This residual dent results from the combination of the residuallocal stress field underneath the load introduction point and the extent of damage inboth the impacted facesheet and the crushed core. Moreover, the residual dent rep-resents a possible source of further damage-growth upon subsequent reloading, and it

53

Chapter 2

is shown to severely reduce the residual local stiffness and strength of the sandwichstructure [10,11,40,41].

Therefore, for damage-tolerant design of foam-cored composite sandwich struc-tures, it is of paramount importance to accurately predict (i) the residual after-crushing strain in the foam core, as it contributes to the final depth of the residualafter-impact dent and (ii) the post-crushing compressive behaviour of the foam core,as it contributes to the residual local stiffness and strength of the impacted sandwichstructure.

2.1.2 Models for the post-crushing compressive response of foammaterials

The typical approaches for modelling the post-crushing compressive response of foammaterials can be gathered into three categories, in decreasing order of accuracy: (i)phenomenological models, (ii) models assuming a linear behaviour governed by adegraded elastic modulus for the elastic regime of the post-crushing response, and(iii) models assuming a linear behaviour governed by the elastic modulus of theundamaged foam material for the elastic regime of the post-crushing response:

(i) Numerous phenomenological models describing the compressive behaviour offoam materials under monotonic compression can be found in literature [42–44].In these works, analytical stress-strain relationships are proposed, whose cal-ibration coefficients are determined to best fit the experimental measurements.The same stress-strain relationships (with modified calibration coefficients) canthen be used for modelling the post-crushing compressive response.

However, this implies that the appropriate calibration coefficients need to bedetermined (usually by least-square fitting to experimental measurements) atany residual strain level. Therefore, the calibration of such phenomenologicalmodels requires that the cyclic compressive response of the investigated foammaterial is experimentally characterized. Unfortunately, although a high num-ber of unloading-reloading cycles would allow for a more accurate descriptionof the post-crushing compressive response, this would lead to practically unaf-fordable testing times, particularly if the investigated foam material exhibits ahighly non-linear behaviour when subjected to cyclic compressive loading.

(ii) Flores-Johnson et al. [45] suggested to model the elastic regime of the post-crushing compressive response of foam materials assuming that the latter exhibita linear behaviour; the authors also presented an analytical model to predictthe degradation of elastic modulus as a function of the residual strain. However,

54

Prediction of the post-crushing compressive response of progressively crushablesandwich foam cores

this model requires, as input, the evolution of the residual strain upon completeunloading as a function of the applied strain at unloading; thus, compressivetests with multiple unloading-reloading cycles are needed for its calibration.

(iii) The most simplistic approach consists in modelling the elastic regime of thepost-crushing compressive response of foam materials assuming that the latterexhibit a linear behaviour governed by the elastic modulus of the undamagedfoam material. Although in many cases inaccurate, this approach is frequentlyused for practical applications [46–49], since its calibration requires to experi-mentally characterize the response of the foam material only under monotoniccompressive loading.

2.1.3 Structure of this chapter

In this chapter, we present a novel analytical model for predicting the post-crushingcompressive response of crushable foams; the model is developed such that its cal-ibration can be performed using exclusively data obtained from standard monotoniccompressive tests. This chapter is organized as follows: the proposed analytical modelis described in § 2.2, while modelling predictions are compared against available ex-perimental data in § 2.4. The results of this comparison are presented and discussed,respectively, in § 2.4.3 and § 2.5. Finally, conclusions are drawn in § 2.6.

2.2 Model development

The following assumptions of the analytical model presented in this thesis need tohighlighted:

(i) the model assumes strain localisation; the underlying hypothesis is that thesandwich foam core is sufficiently thick to be larger than a characteristic ma-terial length-scale, which for ductile foams is typically of the order of few cells;strain localisation may not occur in too thin foam cores;

(ii) the model assumes that the foam core deformation consists of a crushing frontnormal to the loading direction; this is not necessarily true for very thick sand-wich foam cores or for neat foam specimens;

(iii) the model neglects material visco-elasticity and visco-plasticity (strain-ratesensitivity), as well as the effect of air pressure;

(iv) the objective of the model is to predict only the reloading compressive behaviour.

55

Chapter 2

2.2.1 Monotonic compressive response

Let us consider a prismatic crushable foam specimen of thickness h0 and cross-sectional area A0 in the undeformed configuration, as shown in Figure 2.1. Fur-thermore, let the specimen be loaded in compression (under displacement-control),with u being the applied displacement. The resulting true strain in the specimen(homogenized across its entire thickness) is calculated as

〈ε〉=ln(

h0

h0−u

)= ln(

h0

h

), (2.1)

where h = h0 − u is the specimen thickness in the deformed configuration and

〈•〉 =1h

h∫0

(•)dz is the average operator along h, with z indicating the through-the-

thickness direction. Following the definition given in Equation 2.1, throughout thischapter, compressive strains are assumed to be positive.

In this work, the response of crushable foams under monotonic compression ismodelled through a piece-wise continuous constitutive law, relating the nominal com-pressive stress σ to the homogenized strain 〈ε〉. The use of nominal stresses, ratherthen true stresses, is supported by the negligible lateral expansion exhibited by typ-

Undeformedcon¯guration

Deformedcon¯guration

h0h0h

A0 A

z z

u

Figure 2.1: Schematic of a foam specimen under compressive loading. In the un-deformed configuration, the specimen thickness and the cross-sectionalarea are denoted as h0 and A0, respectively; in the deformed config-uration, under the applied displacement u, the specimen thickness andcross-sectional area are denoted, respectively, as h=h0−u and A≈A0.

56


h"i

¾

D C

E

Strain localization

Lock-up

Figure 2.2: Typical nominal stress (σ) v.s. homogenized true strain (〈ε〉) curve fora crushable foam material under monotonic compressive loading. Threemain deformation regimes can be identified: elasticity (E), progressivecrushing (C) and densification (D).

ical foam materials when loaded in compression [50–52], i.e. A0 ≈A (see Figure 2.1);moreover, in the remaining of this chapter the compressive stresses are assumed tobe positive.

The proposed constitutive law individually describes the three main deformationregimes (see Figure 2.2), i.e. initial elasticity (E), progressive crushing (C) and finaldensification (D), as follows.

E: Before crushing initiates, although a narrow region of nonlinear behaviour iscommonly observed immediately prior to strain localization (see Figures 2.1and 2.3a), a linear elastic behaviour is assumed, i.e.

σE = E〈ε〉 , (2.2)

where E is the elastic modulus of the undamaged foam material, and σE is thenominal compressive stress in the elastic regime. In Figure 2.3a, Xc indicatesthe compressive strength (stress at strain localization) of the foam material.Due to the inhomogeneity of foam materials at the micro-scale (cell-scale), Xc

is, effectively, the compressive strength of the weakest layer of material.

57

Chapter 2

h"i

Xc

E

Xc

E

¾E

"c

Crushing initiation

(a) Elastic.

Xc

XL

Xc

E

"L

¾C

h"i

(b) Progressive crushing.

"L~"

L

EL

~EL h"i

¾D

XL

(c) Densification.

Figure 2.3: Piecewise constitutive law for the compressive response of crushablefoam materials: experimental measurements v.s. model predictions(black curves) within the elastic (a), progressive curshing (b) and dens-ification (c) regimes.

The assumption of linearity leads to an underestimation of the strain at crush-ing initiation, as shown in Figure 2.3a, where εc is the homogenized strain atcrushing initiation; however, since the elastic response of the foam material isconfined to very small levels of strain, the effect of the linear idealization isnegligible.

Furthermore, the response of most ductile foams at the cell-level is characterisedby pronounced strain-hardening prior to softening; nonetheless, as the strainhardening is large, plastic strains tend to be very small. Therefore, prior tostrain localisation, it is reasonable to assume that the modulus is equal to theoriginal elastic modulus.

C: The nominal stress during progressive crushing, denoted as σC , is describedthrough a n-th order polynomial, i.e.

σC =n∑

i=0ai

[〈ε〉− Xc

E

]i, (2.3)

where ai are the polynomial coefficients. The coefficients a0 and a1 in Equa-tion 2.3, can be computed by imposing the continuity condition at 〈ε〉=

Xc

Eand

58


〈ε〉=εL , i.e

⎧⎪⎪⎪⎨⎪⎪⎪⎩σC

(〈ε〉 =

Xc

E

)=Xc

σC(〈ε〉 = εL)=XL

=⇒

⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩a0 =Xc

a1 =E(XL − Xc)

EεL −Xc−

n∑i=2

ai

[εL − Xc

E

](i−1) ,

(2.4)where εL is the lock-up or densification strain (at which the crushing processis concluded and the foam material is fully-crushed), while the correspondingnominal stress level (for 〈ε〉 = εL) is denoted as XL and commonly referredto as the lock-up or densification strength. The remaining coefficients ai, withi∈{2,n}, can be obtained through a least-square fitting [53] against the availableexperimental data.

The order n of the polynomial depends on the strain hardening during progress-ive crushing; a 3rd-order polynomial is generally sufficient to accurately fit theexperimental data.

D: It is assumed that the tangent modulus and the density of the fully-crushed foammaterial are related through a power law, i.e.

dσD

d〈ε〉 = k (ρD)μ =⇒ dσD

d〈ε〉 = EL

[ρD(〈ε〉)

ρL

]μ, (2.5)

where σD is the nominal stress in the fully-crushed material, ρD is its density andk is a proportionality constant. Furthermore, EL and ρL denote, respectively,the tangent modulus and the density of the fully-crushed material at σD =0.

By imposing ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩σD(〈ε〉 = εL)=XL

dσD

d〈ε〉∣∣∣∣〈ε〉=εL

= EL

, (2.6)

then Equation 2.5 can be reformulated as

σD =[

EL

μ−XL

][EL

EL −μXLexp[μ(〈ε〉−εL)] − 1

], (2.7)

where EL is the tangent modulus of the foam material at the onset of densi-fication (〈ε〉 = εL). The value of the exponent μ can be determined through aleast-square fitting against the available experimental data. The strain valueindicated as εL in Figure 2.3c, corresponding to the homogenized strain in the

59

Chapter 2

fully-crushed material at σD =0, can be computed from Equation 2.7 as

εL =εL − 1μ

ln[

EL

EL −μXL

]. (2.8)

Therefore, the model assumes a non-linear elastic behaviour of the fully-crushedmaterial; furthermore, the model approximates the response of the fully crushedmaterial as fully-reversible elastic.

To summarize, the compressive response of crushable foams can be describedthrough the following piece-wise constitutive law:

σ =

⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

σE (as defined in Equation 2.2) forXc

E≤ 〈ε〉≤ Xc

E

σC (as defined in Equation 2.3) forXc

E<〈ε〉<εL

σD (as defined in Equation 2.7) forXc

E≤ 〈ε〉 ≥ εL

. (2.9)

2.2.2 Thickness of the uncrushed and crushed layers in crushablefoam materials

Let us consider a specimen of crushable foam material subjected to a homogenizedcompressive strain 〈ε〉, as defined in Equation 2.1; the corresponding compressivestress σ (constant across the entire specimen thickness) is computed according toEquation 2.9.

In the most general case, it is possible to identify, across the specimen thickness,two distinct regions which co-exist at the equilibrium stress σ: (i) a layer of uncrushedfoam material and (ii) a layer of fully-crushed foam material. The thickness of theselayers, respectively indicated as hu and hc in Figure 2.4, are related by

h(〈ε〉)=hu(〈ε〉) + hc(〈ε〉) . (2.10)

The accurate prediction of hu and hc, at any level of applied homogenized strain〈ε〉, is of paramount importance to determine the residual mechanical properties ofthe foam material.

The homogenized strains within the layers of uncrushed and crushed foam mater-ial, respectively denoted as 〈εE〉 and 〈εD〉, are different (see Figure 2.4). In this work,in order to compute 〈εE〉 and 〈εD〉, it is assumed that the uncrushed material behaves

60


according to Equation 2.2, while the crushed material according to Equation 2.7. Atany level of homogenized strain 〈ε〉, the continuity condition reads

h exp [〈ε〉] = hu exp [〈εE〉] + hc exp [〈εD〉] . (2.11)

Figure 2.5 shows the typical evolution of the normalized thicknesses hu =hu

h0(red

curves), hc =hc

h0(green curves) and h =

h

h0(dashed black curves). Similarly to the

previous section, the evolution of hu and hc with 〈ε〉 is individually derived for theelastic (E), progressive crushing (C) and densification (D) deformation regimes, asfollows.

E: For 〈ε〉 ≤ Xc

E, the entire specimen is uncrushed and behaves according to Equa-

tion 2.2. According to Figure 2.5, the normalized thicknesses hu and hc aretherefore

hu(〈ε〉)=h(〈ε〉)=1

exp [〈ε〉] and hc(〈ε〉)=0 . (2.12)

¾

h"i

hu

hc

h"Ei h"

Di

E D C

Figure 2.4: Determination of the thickness of the uncrushed and crushed foam ma-terial layers. For any homogenized strain 〈ε〉, the corresponding com-pressive stress σ is computed through Equation 2.9. At this stress level,the homogenized strains within the uncrushed and the crushed layers,respectively denoted as 〈εE〉 and 〈εD〉, are computed assuming that theuncrushed foam material behaves according to Equation 2.2 and thecrushed material according to Equation 2.7. The thicknesses hu and hccan therefore be computed by exploiting the equilibrium condition, i.e.σ is constant through the entire foam specimen.

61

Chapter 2

0.0

1.0

0.0

1.0

h"iXc

E

"L

hu hc

h

E

C D

Figure 2.5: Typical evolution of the normalized (with respect to the initial specimenthickness h0) thickness of (i) the uncrushed layer hu (red curve), (ii)the crushed layer hc (green curve) and (iii) the entire specimen h (blackdashed curve) as a function of the applied homogenized strain 〈ε〉.

C: ForXc

E<〈ε〉<εL , the homogenized strains 〈εE〉 and 〈εD〉 are computed, respect-

ively from Equation 2.2 and Equation 2.7, as

〈εE〉=σ

Eand 〈εD〉=

1μ

ln[

EL −μXL

EL

(1+

μσ

EL −μXL

)]. (2.13)

By replacing the latter into Equation 2.11 and exploiting Equation 2.10, aftersome mathematical manipulations, the expressions for the normalized thick-nesses hu and hc (see Figure 2.5) come, respectively, as

hu(〈ε〉)=1

exp[〈ε〉] ·[ exp[〈εD〉]−exp[〈ε〉]

exp[〈εD〉]−exp[〈εE〉]]

and

hc(〈ε〉)=1

exp[〈ε〉] ·[ exp[〈ε〉]−exp[〈εE〉]

exp[〈εD〉]−exp[〈εE〉]] . (2.14)

D: For 〈ε〉≥εL , the entire specimen is crushed and behaves according to Equation 2.7.The normalized thicknesses hu and hc are therefore (see Figure 2.5)

hu(〈ε〉)=0 and hu(〈ε〉)=h(〈ε〉)=1

exp [〈ε〉] . (2.15)

62


2.2.3 Post-crushing compressive response and residual strain

The typical σ −〈ε〉 curve for a crushable foam subjected to a complete unloading-reloading cycle is shown in Figure 2.6, where the homogenized strain and nominalstress at unloading initiation are respectively denoted as 〈εUn〉 and σUn . Crushablefoams exhibit an hysteretic behaviour when subjected to cyclic loading, i.e. the un-loading and reloading paths do not coincide. Such hysteretic behaviour is associatedto the inherent viscous response of the parent material [54–58].

The correct determination of the homogenized residual strain upon complete un-loading (indicated as 〈εRe〉 in Figure 2.6) is of paramount importance, along withthe detailed prediction of the post-crushing compressive response of the foam ma-terial (’Reloading’ path in Figure 2.6), for the damage-tolerant design of foam-coredsandwich structures. Within this context, the accurate modelling of the unloadingresponse is, comparatively, of minor relevance; therefore, in this work, it is assumedfor simplicity that the unloading curve coincides with the reloading curve.

Let a crushable foam specimen be loaded in compression to the homogenizedstrain 〈ε〉=〈εUn〉, as shown in Figure 2.7, and subsequently unloaded to the generichomogenized strain 〈ε〉<〈εUn〉. Depending on whether the homogenized strain 〈εUn〉

Reloa

din

g

h"Unih"

Rei

¾Un

Unlo

adin

g

¾

h"i

Unloading initiation

Figure 2.6: Typical response of a crushable foam material subjected to a completeunloading-reloading cycle. The nominal stress and homogenized strainat unloading initiation are denoted, respectively, as σUn and 〈εUn〉. Thereloading is assumed to start immediately at σ = 0 (〈ε〉 = 〈εRe〉), uponcomplete unloading.

63

Chapter 2

¾

h"Ei

hUnu

h"i

h"Uni

¾Un

h"Di

hUnc

E D C

Figure 2.7: Modelling the post-crushing compressive response of crushable foammaterials. The thickness of the uncrushed and crushed layers at unload-ing initiation are denoted as h

Unu and h

Unc , respectively. For 〈ε〉 ≤ 〈εUn〉,

the post-crushing compressive response (orange curve) is computed byimposing the equilibrium between the layers of uncrushed and crushedmaterial, while for 〈ε〉 > 〈εUn〉, it is described by Equation 2.9.

is in the elastic (E), in the progressive crushing (C) or in the densification (D)deformation regime, the expression of the reloading curve and the residual strain atσ =0 can be derived as follows.

E: If 〈εUn〉 ≤ Xc

E, the material is entirely uncrushed at unloading initiation. There-

fore, its post-crushing compressive response can be described by Equation 2.9and 〈εRe〉=0.

C: IfXc

E< 〈εUn〉 < εL , the thicknesses of the uncrushed and crushed layers at un-

loading initiation, respectively denoted as hUnu and h

Unc in Figure 2.7, can be

computed using Equation 2.14. For 〈ε〉 < 〈εUn〉, following the approach pro-posed in the previous section, the post-crushing compressive response of thefoam material (orange curve in Figure 2.7) can be computed by imposing theequilibrium between the uncrushed and crushed layers of foam material, i.e.

1exp [〈ε〉] =hUn

u · exp[

σUn −σ

E

]+ hUn

c ·[

EL +μ (σUn −XL)EL +μ(σ−XL)

]μ. (2.16)

where hUnu =

hUnu

h0and hUn

c =h

Unc

h0are, respectively, the normalized thickness of the

uncrushed and crushed layers at unloading initiation.

64


For 〈εUn〉< 〈ε〉<εL , the post-crushing compressive response of the foam mater-ial can be described by Equation 2.9. The residual strain 〈εRe〉 is computed bysolving Equation 2.16 in 〈ε〉 for σ =0, i.e.

〈εRe〉 = − ln[hUn

u · exp(

σUn

E

)+ hUn

c ·(

EL + μ(σUn −XL)EL −μXL

)μ]. (2.17)

D: If 〈εUn〉 ≥ εL , the material is entirely crushed at unloading initiation; its post-crushing compressive response is described by Equation 2.9 and 〈εRe〉= εL (seeFigure 2.3c).

2.3 Numerical implementation

An overview of the numerical implementation of the proposed model is provided inFigure 2.8; the use of array programming (e.g MATLAB) significantly reduces therunning time and simplifies the numerical implementation. The following remarksshould be highlighted:

(i) the input to the model is the applied homogenized strain, expressed as a C0

continuous function f〈ε〉(k), where the real variable k ∈[0, kmax

]acts as a nu-

merical loading time. The code then samples f〈ε〉 to define a discrete vector ofapplied homogenized strain 〈ε〉, with n〈ε〉 strain values. As output, the modelcalculates the corresponding nominal stress vector σ;

(ii) the model requires, as input, selected properties of the foam material{E, EL , Xc, XL , εL}, as well as the experimental σ−〈ε〉 curve within the crushingregime

[〈ε〉Exp , σExp

C

]and the densification regime

[〈ε〉Exp , σExp

D

];

(iii) the coefficients{

ai

}n

i=0and μ are calculated by means of a nonlinear least-

square fitting of experimental data; such fitting is denoted in Module II.2 andII.3 with the function Λ, defined as

λ = Λ([

xExp, yExp], y = g

λi(x))

. (2.18)

The arguments of Λ are the experimental data to be fitted[xExp, yExp

](where

xExp and yExp are the experimental values of, respectively, the independentvariable x and the dependent variable y) and the analytical expression of thefitting function y = g

λi(x), with coefficients λi; the latter (gathered in the vector

λ in Equation 2.18) represent the outcome of Λ;

65

Chapter 2

I. De¯nition of input variables

I.2: Material properties I.3: Experimental data

II. Preliminary calculations

II.1: De¯ne applied homogenized strain vector

III. Compressive response

III.0: Set step counter and compute maximum applied homogenized strain

III.1a: Main loading curve

if h"ij ¸ h"ij¡1 and h"ij ¸ h"Uni if h"ij < h"Unij = j + 1 and h"Uni = max

³nh"ioj0

´

) ~"L ="L¡1

¹ln

·EL

EL¡¹XL

¸

¾j=Eh"ij

¾j=

nXi=0

ai

·h"ij¡

XcE

¸i

)³hu

j=

1

exp[h"ij]¢"

exp[h"Dij]¡exp[h"ij ]exp[h"Dij ]¡exp[h"Eij ]

#,

³hc

j=

1

exp[h"ij ]¢"

exp[h"ij ]¡exp h"Eij ]exp[h"Dij]¡exp[h"Eij ]

#

III.1b: Cyclic loading curve

III.1a(C): if Xc=E<h"ij<"L (Crushing regime)

III.1a(D): if h"ij¸"L (Densi¯cation regime)

¾j=

·EL

¹¡XL

·EL

EL¡¹XL

exp[¹(h"ij¡"L)]¡ 1

¸

) h"Eij=

¾jE

, h"Dij=

1

¹ln

·EL¡¹XL

EL

μ1+

¹¾jEL¡¹XL

¶¸

with a0 = Xc ; a1 =E(XL ¡Xc)E"L¡Xc

¡nXi=2

ai

·"L¡XcE

¸(i¡1)

(Unloading within elastic regime)

¾j=Eh"ij )³hu

j=

1

exp[h"ij];³hc

j=0

)³hu

j=0 ;

³hc

j=

1

exp[h"ij]

) h"Eij=h"ij ; h"

Dij=0

) h"Eij=

¾jE; h"

Dij=

1

¹ln

·EL¡¹XL

EL

μ1+

¹¾jEL¡¹XL

¶¸

I.1: Numerical variables

III.1b(C): if Xc=E<h"Uni<"L(Unloading within crushing regime)

¾Un

=

nXi=0

ai

·h"

Uni¡Xc

E

¸i

) hUnu =1

exp[h"Uni] ¢

26664·1+

¹

EL

(¾Un¡X

L)

¸1=¹¡exp[h"

Uni]·

1+¹

EL

(¾Un¡XL)

¸1=¹¡exp

h¾UnE

i37775 ;

if j < nh"iif j < nh"i

¹ = ¤

Ãhh"iExp ;¾Exp

D

i; ¾D =

·EL

¹¡XL

·EL

EL¡¹XL

exp[¹(h"i¡"L)]¡ 1

¸!nai

oni=0

= ¤

Ãhh"iExp ;¾Exp

C

i; ¾C =

nXi=0

ai

·h"i¡Xc

E

¸i!

hh"iExp ;¾Exp

C

i;hh"iExp ;¾Exp

D

i

j = 0 )h"i0 = 0 ¾0 = 0

E;EL;Xc;XL

; "L

¢k =

¯kmaxnh"i ¡ 1

¯n;nh"i; fh"i(k); kmax

II.2: De¯ne compressive response coe±cients - Progressive crushing regime

II.3: De¯ne compressive response coe±cients - Densi¯cation regime

III.1a(E): if h"ij · Xc=E (Elastic regime) III.1b(E): if h"Uni · Xc=E

III.1b(D): if h"Uni ¸ "

L

(Unloading within densi¯cation regime)

¾j=

·EL

¹¡XL

·EL

EL¡¹XL

exp[¹(h"ij¡"L)]¡ 1

¸!

)³hu

j=0 ;

³hc

j=

1

exp[h"ij]

hUnc =1

exp[h"Uni]¢

26664 exp[h"Uni]¡exp

h¾UnE

i·1+

¹

EL

(¾Un¡X

L)

¸1=¹¡exp

h¾Un

E

i37775

h"i =nfh"i(j ¢¢k)

oj=0

nh"i¡1

¾j=¾

¯1

exphh"ij

i=hUnu ¢exp

·¾Un¡¾E

¸+ hUnc ¢

·EL+¹ (¾

Un¡X

L)

EL+¹(¾¡XL)

¸¹

)³hu

j=hUnu ¢exp

·¾Un¡¾jE

¸,

³hc

j=hUnc ¢

·EL+¹ (¾

Un¡X

L)

EL+¹(¾j¡XL)

¸¹

)³hu

j=

1

exp[h"ij];³hc

j=0

Figure 2.8: Numerical implementation of the proposed model.

66


(iv) the homogenized strains 〈εE〉 and 〈εD〉, as well as the normalized thicknesseshu and hc, become discrete vectors (indicated in upright bold in Figure 2.8).

2.4 Model validation

2.4.1 Characterization of the monotonic and cyclic compressive re-sponse of the Rohacell HERO 71 foam

2.4.1.1 Material

The material analysed experimentally in this work is the polymeric foam ROHA-CELL HERO 71 by Evonik Industries, a fully-isotropic and closed-cell rigid PMI(polymethacrylimide) foam. The foam has a nominal relative density ρ∗ = 0.0625,computed as the ratio of the foam density ρf = 75 kg/mm3 and the PMI densityρPMI =1200 kg/mm3 [51].

2.4.1.2 Cyclic compressive tests

Six specimens with nominal in-plane dimensions equal to 20.0×20.0 mm2 were cutfrom 16.3 mm-thick panels using a wire saw. For each panel, a reference surface wasidentified so that every specimen cut from the same panel could be identically orientedwhen tested. The specimens were conditioned in accordance with Procedure C of theASTM D5229 standard [59].

The crushing response of the HERO G3 foam was analysed using flatwise com-pressive tests as prescribed by the ASTM C365-57 standard [60]. The specimenswere positioned between two flat steel plates and were tested in compression usingan INSTRON 5969 servo-hydraulic machine with a 50 kN load cell at a displacementrate u=0.50 mm/min, corresponding to an engineering strain rate e≈511 μs−1. Theloading platens were coated with silicon spray (PTFE) to minimize frictional effects.

To characterize the response of the HERO foam under compressive cyclic load-ing, four quasi-static complete unloading-reloading cycles were performed during thecrushing at intervals of 0.15 engineering strain (e ≈ 511 μs−1 during both unloadingand subsequent reloading).

67

Chapter 2

2.4.2 Measured foam properties

The model predictions are compared against experimental results for three differentfoam materials, two of these from literature (PVC closed-cell Divinycell H100 [45]and PMI Rohacell WF51 [45]) and a third one characterised specifically for thisvalidation (PMI Rohacell HERO 71, see § 2.4.1). Table 2.1 summarizes the measuredproperties of the foam materials investigated, along with the parameters required forthe calibration of the proposed model.

Table 2.1: Measured properties and model calibration parameters for the RohacellWF51 [45], Divinycell H100 [45] and Rohacell HERO 71 foams.

E EL Xc XL εL ai, i∈ [1, 3] μ

WF51 44.58 21.7120 1.0943 0.9753 0.5235

⎡⎣−5.712112.9557−6.1876

⎤⎦ 5.1454

(MPa) (MPa) (MPa) (MPa) (-) (MPa) (-)

H100 35.70 27.8949 1.5000 3.0068 0.5407

⎡⎣ 22.4890−89.0717130.6767

⎤⎦ 4.9839


HERO 71 39.0701 18.0100 0.9716 2.1067 0.5268

⎡⎣ 8.8664−20.670426.5687

⎤⎦ 4.6577


When a clear onset of densification cannot be easily identified on the σ−〈ε〉 curve,the determination of XL will be semi-arbitrary. This choice will, admittedly, have aninfluence on the predicted reloading compressive modulus. Unlike the PMI foams,for the PVC Divinycell H100 foam, the onset of densification is not characterized bya clear knee in the σ−〈ε〉 curve; therefore, following Arezoo et al. [52], for this foammaterial, the lock-up strain εL is defined as the homogenized compressive strain atσ =2Xc.

2.4.3 Model predictions

2.4.3.1 Thickness of the crushed layer

Experimental measurements for the thickness of the crushed layer as function ofthe applied homogenized strain 〈ε〉 are available only for the HERO 71 foam. Thethickness of the crushed layer was determined through the analysis of the strainfields across the foam specimens during progressive crushing; such strain fields areobtained using the Digital Image Correlation (DIC) technique. Figure 2.9 compares

68


these against model predictions. Since the thickness of the uncrushed and crushedlayers are, for any applied stain 〈ε〉, related by Equation 2.10, only experimentalmeasurements of hc(〈ε〉) are shown in Figure 2.9.

2.4.3.2 Residual strain

Figure 2.10 compares model predictions against experimental measurements for theresidual strain 〈εRe〉 as a function of the strain at unloading initiation 〈εUn〉, for thethree foam materials considered in this work. Furthermore, the reduction of the error

0.00.0

0.00.0

h"i

hu

hc

1.0

h

E

C D

Figure 2.9: Thickness of the layer of crushed material: model predictions v.s. ex-perimental measurements. The normalized thicknesses hu and hc arecalculated by analysing the discontinuous strain field within the speci-men during crushing, using the DIC technique.

69

Chapter 2

(compared to models assuming a linear elastic post-crushing response) in predictingthe residual strain 〈εRe〉 and denoted as Δξ, is shown in Figure 2.11.

h"Uni

h"Rei

h"Ui

Rei Exp. data - WF51

Exp. data - HERO 71

Mod. pred. - WF51

Mod. pred. - HERO 71Mod. pred. - H100

Exp. data - H100

Figure 2.10: Residual strain upon complete unloading: model predictions (solidcurves) against experimental measurements, for the WF51, H100 andHERO 71 foams.

100

60

80

40

h"Uni

20

0

¢»(%)

WF51

HERO 71H100

WF51

HERO71

H100

Figure 2.11: Percentage reduction of the error in predicting the residual strainupon complete unloading 〈εRe〉 using the proposed model (comparedto models assuming a linear [governed by the elastic modulus of theundamaged material] post-crushing compressive response). The aver-age values of the error-reduction obtained for the three foam materialsconsidered in this work are displayed as dashed lines.

70


The error-reduction Δξ is defined as

Δξ =ξM − ξL

ξL, (2.19)

where ξL and ξM represent, respectively, the errors in predicting the residual strain,assuming a linear elastic post-crushing compressive response (elastic modulus of theundamaged material) and using the model proposed in this work. These errors canbe defined as

ξL =

∣∣∣∣∣∣1 − E ·〈εUn〉 − σUn

E ·⟨εExp

Re

⟩∣∣∣∣∣∣ and ξM =

∣∣∣∣∣∣1 −⟨

εNum

Re

⟩⟨

εExpRe

⟩∣∣∣∣∣∣ , (2.20)

where the residual strain⟨

εNum

Re

⟩is computed according to Equation 2.17 and

⟨ε

Exp

Re

⟩is the experimentally measured residual strain. In Figure 2.11, the average values ofΔξ for the WF51, H100 and HERO 71 are indicated.

2.4.3.3 Post-crushing compressive response

In Figure 2.12, model predictions for the post-crushing compressive response arecompared against experimental results (Figure 2.12a for the WF51 foam, Figure 2.12bfor the H110 foam and Figure 2.12c for the HERO 71). Here, the unloading branchesare not displayed for clarity. Moreover, the predictions obtained with the analyticalmodel presented in this chapter are shown as solid black curves, while dashed blackcurves indicate the predicted post-crushing compressive behaviour if a linear response(with the undamaged elastic modulus) is assumed.

¾(MPa)

Linear unloading-reloading(Undamaged elastic modulus)

h"i

(a) Rohacell WF51 foam. Experimental data from Flores-Johnson et al. [45].

71

Chapter 2

h"i

¾(MPa)

Linear unloading-reloading(Undamaged elastic modulus)

(b) Divinycell H100 foam. Experimental data from Flores-Johnson et al. [45].

¾(MPa)¾

(MPa)Linear unloading-reloading(Undamaged elastic modulus)

h"i

(c) Rohacell HERO 71 foam.

Figure 2.12: Comparison of model predictions (black solid curves) against exper-imental measurements for the post-crushing compressive response ofthe PMI Rohacell WF51 (a), PVC Divinycell H100 (b) and PMI Ro-hacell HERO 71 (c) foams (dashed black lines indicate the predictedresponse if a linear behaviour (with the undamaged elastic modulus)was assumed).

72


2.5 Discussion

2.5.1 Model calibration

From § 2.2, it follows that the calibration of the present model is performed using ex-perimental measurements obtained exclusively from standard monotonic compressivetests. Consequently, since the post-crushing response is not replicated (as for exampleusing phenomenological models, see § 2.1) but rather predicted, the need for carryingout time-consuming compressive tests including multiple unloading-reloading cyclesis avoided. Therefore, the time required for the experimental characterization of foammaterials can be significantly reduced.

Indicatively, the experimental characterization of the cyclic compressive responseof the HERO G3 71 foam (four unloading-reloading cycles) described in § 2.4.1 re-quires approximately 44 additional minutes per specimen (≈ +157 %) compared tothe case when only the monotonic compressive response is characterized.

2.5.2 Thickness of the crushed layer

Figure 2.9 exhibits a good agreement between the predicted and measured thicknessof the crushed layer, with the maximum error being approximately equal to 13.6 %).Furthermore, unlike any other model available in literature that we are aware of, theformulation proposed in § 2.2 captures the effect of strain hardening during crushingon the variation of the crushed layer thickness.

2.5.3 Residual strain

Figure 2.10 shows that the present model accurately predicts the value of residualstrain 〈εRe〉 upon complete unloading (σ =0) as a function of the strain at unloadinginitiation 〈εUn〉. The residual strain is predicted with a maximum error of about1.5 %, 8.2 % and 12.1 %, respectively, for the Rohacell WF51, Divinycell H100 andRohacell HERO 71 foams.

According to the formulation presented in § 2.2.3, the predicted residual straindoes not account for visco-elastic/plastic strain relaxation effects. For a meaning-ful comparison, in the experiments, both those available in literature for the WF51and H100 foams [45], and those performed as part of the work presented in thischapter (§ 2.4.1), the compressive reloading starts immediately upon complete un-

73

Chapter 2

loading (σ =0), thus not allowing for any visco-elastic/plastic strain relaxation effectsto take place.

When compared to models assuming a linear elastic behaviour (governed by theelastic modulus of the undamaged material) for both unloading and subsequent re-loading, as shown in Figure 2.11, the proposed model is shown to strongly reduce theerror in predicting the residual strain 〈εRe〉. The latter is reduced, on average, byabout 56.6 %, 70.8 % and 68.3 % for the WF51, H100 and HERO 71 foams, respect-ively.

Although the specific values summarized in this section are, admittedly, dependenton the chosen foam materials, they confirm the capabilities of the approach proposedin this chapter.

2.5.4 Post-crushing compressive response

The comparison of model predictions against experimental measurements (Fig-ure 2.12) shows the capabilities of the present model in predicting the post-crushingcompressive response of crushable foam materials. Noticeably, the formulation pro-posed in § 2.2.3 accounts separately for both the linear-elastic contribution of theuncrushed material layer and the nonlinear contribution of the crushed material layerto the overall post-crushing compressive response. Thus, unlike for other modelswhose calibration requires the same experimental measurements, the model proposedin this chapter captures the increasingly nonlinear response for large values of strainat unloading initiation.

2.6 Conclusions

In this chapter, a novel analytical model for predicting the post-crushing compressiveresponse of crushable foams is presented. The calibration of the model is performedusing experimental measurements obtained exclusively from standard monotonic com-pressive tests; therefore, the need for performing time-consuming compressive testsincluding multiple unloading-reloading cycles is avoided and the effective testing timesignificantly reduced.

Model predictions were validated against experimental measurements for threedifferent foam materials (two from the literature, and one originally presented in thischapter). The model is shown to accurately predict the thickness of the crushedmaterial layer during progressive crushing (maximum error about 13.6%) and, in ad-dition, the residual after-crushing strain (maximum error ranging from approximately

74


1.5% for the WF51 foam to 12.1% for the HERO 71 foam). If compared to otheranalytical models whose calibration requires the same experimental measurements,the present model can predict the residual after-crushing strain with a significantlysmaller error, i.e. with an error-reduction ranging from approximately 56.6% for theWF51 foam to 70.8% for the H100 foam. Furthermore, it is shown that the proposedmodel is able to capture the characteristic features of the post-crushing compress-ive response of crushable foam materials, such as the increasingly nonlinear responseexhibited by the latter for large values of strain at unloading initiation.

The results presented in this chapter demonstrate the relevance of the proposedmodel for damage-tolerant design of foam-cored composite sandwich structures.

75

Chapter 3

Translaminar fracture toughnessof NCF composites withmultiaxial blankets

3.1 Introduction

As a result of the increasing share of composite materials in sectors where cycle timesand manufacturing costs significantly impact the final product cost, the developmentof inexpensive and automated production methods is crucial [12–15].

The need for cost-effective alternatives to conventional prepreg-based compositesled to the development of Non-Crimp Fabric (NCF) composites [16–18]. When com-pared to their prepreg-based counterpart, NCF composites offer higher depositionrates, reduced labour time, higher degree of tailorability and improved impact prop-erties [19,20,25]. Therefore, NCF composites are widely regarded as one of the mostpromising technologies for both aerospace [21, 22] and automotive [23, 24] structuralcomposites. The growing industrial interest towards NCF composites led to two EU-funded research projects: FALCOM [25, 26] and TECABS [27, 28], respectively foraerospace and automotive applications.

Owing to the complex micro-structure of NCF composites, Finite Element (FE)models have been extensively used to investigate their mechanical response at differ-ent length-scales [61–64]. Physically-based failure criteria for NCF composites havebeen proposed by several authors [25, 65]. Nonetheless, these criteria have not yetgained general acceptance due to the extremely high number of input parametersthey require, as well as the difficulties in measuring the latter.

77

Chapter 3

Alternatively, state-of-the-art physically-based criteria developed for conventionalunidirectional (UD) composites [66–69] can, in principle, be applied to the analysis ofNCF composites. However, these criteria do not account for the transverse orthotropyof NCF composites [70]; hence, they require further developments to be capable ofaccurately predicting relevant failure mechanisms in NCF composites subjected tocomplex 3D stress states. Molker et al. [29,30] have proposed a novel set of physically-based failure criteria for NCF composites, based on LaRC05 [71], which account fortheir transverse orthotropy with an additional failure mode.

Physically-based failure criteria can predict failure at the ply-level and require, asinput data, homogenised ply properties which can be measured mostly from standardtests. Particularly, translaminar fracture toughnesses are paramount for the damage-tolerant design of composite structures. Numerous studies have been carried outto characterise the translaminar fracture toughness of UD-ply prepreg composites,e.g. glass/epoxy laminates [72, 73], E-glass fibre-reinforced epoxy laminates [74] andcarbon/epoxy laminates [75,76], as well as of woven composites [77,78]. However, tothe knowledge of the author, no work has been published on the measurement of thetranslaminar fracture toughness of NCF composites.

To address this, the translaminar fracture toughness of a carbon-epoxy NCF com-posite laminate with triaxial ([45◦/0◦/−45◦]) blankets is experimentally measured inthis work. The translaminar fracture toughness of both the individual UD fibre towsand of the triaxial NCF blanket are determined and the concept of a homogenisedNCF blanket-level translaminar fracture toughness was introduced. Furthermore, us-ing an approach developed for UD-ply prepreg composites [79], it is demonstratedthat the translaminar fracture toughness of off-axis fibre tows/NCF blankets can beanalytically related to that of axially-loaded fibre tows/NCF blankets.

The present chapter is organized as follows: the experimental method and the datareduction scheme for the analysis of experimental results are described, respectively,in § 3.2 and § 3.3; experimental results are presented and discussed in § 3.4. Finally,the main conclusions are drawn in § 3.5.

3.2 Experimental method

3.2.1 Material system

The material used in this work is a triaxial NCF composite produced by SaertexGmbH consisting of Toho Tenax HTS fibres and a polyester knitting yarn, infusedwith Hexcel RTM6 epoxy resin. The layup of the triaxial NCF blankets, expressed

78

Translaminar fracture toughness of NCF composites with multiaxial blankets

in terms of the UD fibre tows, is [45◦/0◦/ − 45◦], and their nominal thickness isequal to 0.375 mm (the thickness of the individual fibre tows is 0.125 mm). Thenominal membrane properties of the individual fibre tows are listed in Table 3.1;here, subscripts 1 and 2 denote longitudinal and transverse direction of the fibretows.

Table 3.1: Nominal membrane properties of the fibre tows in the triaxial NCFblanket [80,81].

E1 [GPa] E2 [GPa] G12 [GPa] ν12 [−]130.00 9.00 4.50 0.26

3.2.2 Specimen and layup configuration

Compact Tension (CT) specimens [75, 82], with dimensions shown in Figure 3.1 andlayups provided in Table 3.2 were cut using a CNC water-jet cutter. The notchesof the specimens were machined using a diamond coated disk-saw to guarantee anaccurate and sharp crack tip [83].

0±+45±¡45±

P

P

Á 8

14 a0 = 26 25t

2860

Figure 3.1: CT specimens nominal dimensions (in mm) and fibre directions.

In Table 3.2, each layup is expressed as:

• a tow-level layup, defined considering the orientations of the individual UD fibretows within the triaxial NCF blanket;

79

Chapter 3

• a blanket-level layup, defined homogenising the triaxial NCF blankets as UDlayers oriented as their 0◦ fibre tows.

The translaminar fracture toughness of the NCF laminates is denoted as GAIc

(laminates with layup A) and GBIc (laminates with layup B). Furthermore, we define:

• a tow-level translaminar fracture toughness (i.e. the translaminar fracturetoughness of the individual UD fibre tows within the NCF blankets), denotedas G0

Ic and G45Ic respectively for the 0◦ and 45◦ fibre tows;

• a blanket-level translaminar fracture toughness (i.e the translaminar fracturetoughness of the entire NCF blankets homogenised as UD layers), denoted asGNCF

Ic .

Table 3.2: Layups investigated. The nominal thickness of the laminates is indicatedas tLam and the 0◦ fibre tows are aligned with the direction of the appliedload.

LayupLayup ID

tLam Purpose[mm]

Tow-level Blanket-level of layup

A 6.0 [(45◦/0◦/−45◦)s]8 [0◦]16 GNCFIc , G0

Ic, G45Ic

B 6.0 [(90◦/45◦/0◦2/−45◦/90◦)s]4 [(45◦/−45◦)s]4 G0

Ic, G45Ic

3.2.3 Test method and experimental setup

At least five CT specimens were tested for each layup indicated in Table 3.2 usingan Instron machine with a 20 kN load cell; the applied displacement rate was equalto 0.5 mm/min. A video strain gage system was used to measure and record therelative displacement d of two target points drawn on the surface of the specimens(see Figure 3.2). Load measurements were recorded via the Instron load frame andsynchronized with the relative displacement of the two target points measured by thevideo strain gage system.

80


Top target point

Bottom target point

Figure 3.2: Test set up with target points and scale.

3.3 Data reduction

3.3.1 NCF laminate-level translaminar fracture toughness

The modified compliance calibration (MCC) method [84] was used to calculate theNCF laminate translaminar fracture toughness. Unlike other data reduction schemes,the MCC method does not require optical measurement of the crack length, thereforereducing the operator-dependence of the results. In addition, for the analysis of lam-inates with different ply orientations, not using the optically measured crack positionon the surface is important as the external plies of the specimen often do not reflectthe actual crack front within the specimen during crack propagation, i.e. the crackfront is not necessarily uniform across the specimen thickness.

The MCC method requires the elastic compliance C of the CT specimen to bedetermined at several values of the crack length a. For each of the layups in Table 3.2,an FE model of a half CT specimen (exploiting symmetry) was created in Abaqus [85].Square 8-noded (S8R5) shell elements with side e = 0.5 mm were used. The shapeof the initial notch is not explicitly modelled, as the stress intensity factor is notsignificantly affected by the morphology of the initial opening [86].

A 1 N load was applied at the position of the loading pin. The compliance calib-ration curve C vs. a was obtained in 0.5 mm increments of the initial crack length(across the whole potential crack growth length). The C vs. a data were fitted witha function of the form [87]

C(a) = (αa + β)χ , (3.1)

81

Chapter 3

where α, β and χ were calculated to best fit the experimental data for each layup;the values of α, β and χ for the two layups investigated in this work are provided inTable 3.3. The compliance calibration curves obtained with FE and the correspondingMCC method fitting curves are shown in Figure 3.3a (Layup A) and Figure 3.3b(Layup B).

Table 3.3: Numerical fitting parameters used in the MCC method (units system:kN; mm).

Layup ID α β χ

A −8.944 × 10−2 4.639 -2.192

B −9.338 × 10−2 4.796 -2.205

Therefore, an effective crack length aeff can be determined using the elastic com-pliance computed from the load vs. displacement curve as

aeff =C

1χ − β

α. (3.2)

0

0.2

0.4

0.6

0.8

1

20 30 40

C [mm/kN]

[mm]

FE

Series1

FE

MCC curve (Eq. 1)

a

(a) Layup A (see Table 3.2).

0

0.2

0.4

0.6

0.8

1

20 30 40

C [mm/kN]

[mm]

Series2

Series1

FE

MCC curve (Eq. 1)

a

(b) Layup B (see Table 3.2).

Figure 3.3: Compliance calibration curves obtained from FE and the MCC method.

82


Finally, the translaminar fracture toughness of each layup can be calculated as

GLamIc =

P 2c

2t

dC

da, (3.3)

where Pc is the measured load that propagates the crack.

3.3.2 Fibre tow-level translaminar fracture toughness

3.3.2.1 Relating the toughness of the individual fibre tows to the tough-ness of the laminate

The fracture toughness of the individual UD fibre tows can be obtained from that ofthe NCF laminate using a rule of mixture [75, 88]. Thus, the translaminar fracturetoughness of the layups investigated in this work, and respectively denoted as GA

Ic

and GBIc, can be expressed as

GAIc =

t0A

tAG0

Ic +t45AtA

G45Ic and (3.4)

GBIc =

t0B

tBG0

Ic +t45BtB

G45Ic +

t90BtB

G90Ic . (3.5)

where:

• G0Ic, G45

Ic and G90Ic are, respectively, the translaminar fracture toughness of the

0◦, 45◦ and 90◦ fibre tows;

• t0K, t45

K and t90K are, respectively, the total thicknesses of the 0◦, 45◦ and 90◦ fibre

tows within specimens with generic layup K.

Therefore, assuming that the intralaminar fracture toughness G90Ic is negligible

when compared to the translaminar fracture toughness of the 0◦ and 45◦ fibre tows(G90

Ic << G0Ic, G45

Ic ) [75], G0Ic and G45

Ic can be obtained, respectively, as

G0Ic =[

tAt45B

t0At45

B − t45A t0

B

]·GA

Ic −[

t45A tB

t0At45

B − t45A t0

B

]·GB

Ic , (3.6)

G45Ic =

[t0AtB

t0At45

B − t45A t0

B

]·GB

Ic −[

tAt0B

t0At45

B − t45A t0

B

]·GA

Ic . (3.7)

83

Chapter 3

3.3.2.2 Relating the toughness of the off-axis fibre tows to that of the 0◦

and 90◦ fibre tows

Teixeira [79] investigated the crack propagation across off-axis plies in prepreg-basedCFRPs and showed that a crack propagates across 45◦ plies through a combinationof tensile fibre failure (as in 0◦ plies under translaminar tension (dashed blue curvesFigure 3.4)) and splits in between fibres (as in 90◦ plies under intralaminar tension(dashed red curves in Figure 3.4)). Thus, the translaminar fracture toughness ofoff-axis plies can, in principle, be expressed as a function of the 0◦ plies translaminarfracture toughness and of the 90◦ plies intralaminar toughness.

Therefore, from geometrical considerations, the translaminar fracture toughnessof off-axis plies (at an angle α with respect to the 0◦ fibre tows), denoted as Gα

Ic, canbe estimated as [79]

GαIc = cos(α) · G0

Ic + sin(α) · G90Ic . (3.8)

In the present work, we use the same relation for off-axis fibre tows in an NCF ar-chitecture. Therefore, rearranging Equation 3.8 and neglecting G90

Ic , the translaminarfracture toughness of the 0◦ plies can be independently calculated from the translam-

Tensile fibre failure

Splits between fibre

0 ±

®

(G 0Ic

) (G 90Ic

)

Macroscopic crack-propagation

direction

Figure 3.4: Micrograph of a crack propagating across off-axis plies (angle α) in aprepreg-based composite laminate, after [79].

84


inar fracture toughnesses of layups A and B, i.e.

G0Ic =

⎡⎢⎢⎣ tK

t0K +

√2

2·t45

K

⎤⎥⎥⎦·GKIc , K ∈ {A,B} . (3.9)

3.3.3 NCF blanket-level translaminar fracture toughness

Finite Element models of NCF composite laminates are often created by modellingthe multi-axial NCF blankets as a single layer of material with homogenised proper-ties. Therefore, physically-based failure criteria as those reviewed in § 3.1 require, inaddition to the translaminar fracture toughness of the individual UD fibre tows, alsothe NCF blanket-level translaminar fracture toughness.

The translaminar fracture toughness of the triaxial NCF blanket investigated inthis work can be directly evaluated from the measured translaminar fracture tough-ness of layup A, i.e.

GNCFIc = GA

Ic (3.10)

or, following the approach detailed in § 3.3.2.2, from the measured translaminarfracture toughness of layup B, i.e.

GNCFIc =

√2 ·GB

Ic . (3.11)

3.4 Results and discussion

3.4.1 Load displacement curves

Figure 3.5 shows the experimental load vs. displacement curves for layup A (Fig-ure 3.5a) and layup B (Figure 3.5b). All the CT specimens tested exhibited a stick-slip crack-growth during testing. Initial failure of the CT specimens was taken asthe first significant load-drop in the load-displacement curves. Final failure of theCT specimens corresponded to compressive failure near the edge opposite to the loadapplication (last significant load-drop in the load-displacement curves).

85

Chapter 3

0

2

4

6

8

10

12

0 1 2 3 4

d [mm]

P [kN]

Pmax = 10 :2 (2 :6% )


0

2

4

6

8

10

0 1 2 3 4

d [mm]

P [kN]

Pmax = 9 :3 (3 :3% )


Figure 3.5: Experimental load (P ) vs. displacement (d) curves for the layups in-vestigated.

3.4.2 Translaminar fracture toughness

3.4.2.1 NCF laminate-level translaminar fracture toughness

Figure 3.6 shows the R-curves for layup A (Figure 3.6a) and layup B (Figure 3.6b).Fracture toughness initiation values are defined as the intersection between the dashedlines at an angle and the vertical axes; although an R-curve effect could be inferred, nomeaningful propagation values were obtained as a result of the premature compressivefailure of the CT specimens. The average values of the translaminar initiation fracturetoughness of layup A and B are provided in Figure 3.7; the corresponding coefficientsof variation are provided in brackets.

3.4.2.2 Fibre tow-level translaminar fracture toughness

The average values of the translaminar fracture toughness for the UD fibre towsare shown in Figure 3.8; the corresponding coefficients of variation are provided inbrackets. The leftmost and rightmost columns indicate the fracture toughness ofthe 0◦ and 45◦ fibre tows obtained using, respectively, Equations 3.6 and 3.7, i.e.from both GA

Ic and GBIc. The second and third columns (from the left) indicate the

fracture toughness of the 0◦ fibre tows predicted independently from GAIc and GB

Ic,using Equation 3.9.

86


0

100

200

300

400

0 2 4 6 8 10

Initiation

¢a [mm]

GLamIc [kJ=m2]


0

50

100

150

200

0 2 4 6 8 10

Initiation

¢a [mm]

GLamIc [kJ=m2]


Figure 3.6: R-curves for NCF laminate-level translaminar fracture toughness.

0

50

100

150

GLamIc [kJ=m2]

G AIc

G BIc

128 (4 %)

91 (7 %)

Figure 3.7: NCF laminate-level translaminar fracture toughnesses (initiation val-ues).

Quantitatively, the difference between the average value of G0Ic computed using

both GAIc and GB

Ic (see Equation 3.6) and that computed exclusively from GAIc is approx-

imately equal to 5 %: furthermore, the difference between the value of G0Ic computed

using both GAIc and GB

Ic and that computed exclusively from GBIc is approximately

equal to 3 %. With regards to the scatter, because Equation 3.9 uses informationfrom both the 0◦ fibre tows and the 45◦ fibre tows, while Equation 3.6 uses only

87

Chapter 3

information from the 0◦ fibre tows, the scatter in the measurements is significantlyreduced by using the analytical model expressed in Equation 3.8.

Therefore, knowing the translaminar fracture toughness of the 0◦ fibre tows theapproach outlined in § 3.3.2.2 allows accurate prediction of the translaminar fracturetoughness of off-axis fibre tows. Hence, only the translaminar fracture toughnessof the 0◦ fibre tows may be needed to estimate the translaminar fracture tough-ness of laminates with complex layups (with several differently-oriented off-axis fibretows). Although further verification is needed for other values of the angle α (seeEquation 3.8), this result is particularly relevant for the design of NCF compositelaminates to be used in large-scale structural applications.

3.4.2.3 NCF blanket-level translaminar fracture toughness

Figure 3.9 shows the average values of the translaminar fracture toughness of theNCF blanket computed according to Equation 3.10 and Equation 3.11 (approximatedifference of 2 %); the corresponding coefficients of variation are provided in brackets.This result confirms the validity of the approach detailed in § 3.3.2.2 also for theprediction of the translaminar fracture toughness of off-axis NCF blankets.

0

50

100

150

200

GIc [kJ=m2]

G 0Ic

G 0Ic

G 0Ic

G 45Ic

Eq. 6 Eq. 7 Eq. 9 (K=A)

Eq. 9 (K=B)

166 (25 %)

158 (3 %)

161 (7 %)

108 (22 %)

Figure 3.8: Fibre tow-level translaminar fracture toughness (initiation value).

88


0

50

100

150

G NCFIc

[kJ=m 2 ]

128 (4 %)

129 (9 %)

(Eq. 10) (Eq. 11)

Figure 3.9: NCF blanket-level translaminar fracture toughness (initiation value).

3.5 Conclusions

In this work, the translaminar initiation fracture toughness of a carbon-epoxy NCFcomposite laminate with triaxial ([45◦/0◦/−45◦]) blankets was measured using a Com-pact Tension (CT) test; no meaningful propagation values could be determined, as aresult of premature compressive failure of the CT specimens. The translaminar frac-ture toughness of the individual UD fibre tows was related to that of the NCF laminateand the concept of an homogenised blanket-level translaminar fracture toughness wasintroduced.

In this work, translaminar fracture toughness values were computed neglectingthe effect of possible delaminations (inter- and intra-blanket), see Section 3.2.1. SinceNCF blankets are stacked such that adjacent fibre tows belonging to different blanketshave the same orientation (for both Layup A and Layup B), inter-laminar delamina-tions are prevented. Moreover, the transverse stitching yarns inhibit the propagationof intra-laminar delaminations.

Using an approach developed for UD-ply prepreg composites [79], it is demon-strated that the translaminar fracture toughness of off-axis fibre tows/NCF blanketscan be analytically related to that of axially-loaded fibre tows/NCF blankets. Thepercentage difference between the values obtained experimentally and those predictedusing such analytical approach is, for the material system investigated in this work,lower than 5%. Furthermore, using the analytical model allows to reduce significantlythe scatter in the measurements.

89

Chapter 3

Therefore, the translaminar fracture toughness of laminates with complex layups(with several differently-oriented off-axis fibre tows and off-axis NCF blankets) can beaccurately estimated from the translaminar fracture toughness of axially-loaded fibretows/NCF blankets. This result is highly relevant for the design of NCF compositelaminates to be used in large-scale structural applications.

90

Chapter 4

Multiple length/time-scalesimulation of localized damagein composite structures using aMesh Superposition Technique

4.1 Introduction

4.1.1 Motivation

Numerical simulation of the mechanical response of large composite components oftenrequires that different parts of the structure are modelled at different scales, eventuallyeven using different physics. Multiscale modelling techniques can achieve the requiredlevel of accuracy in each part of the model while maintaining the computational timeto a minimum. Within this framework, the development of suitable techniques forcoupling areas of the structure discretized using different element types is crucial.

However, the coupling of subdomains discretized with finite elements of differentphysical dimension/formulation can introduce artificial stresses at the shared bound-aries [89, 90]. Therefore, the stress field within the structure and its mechanicalresponse may not be correctly simulated; as a result, in problems involving failure,the damage pattern might not be faithfully replicated. Additionally, in dynamic prob-lems, the interfaces between differently-discretized subdomains may artificially reflectstress waves [91,92].

91

Chapter 4

4.1.2 3D solid elements and 2D shell elements for laminated com-posites

The typical finite elements used for modelling composite structures can be gatheredinto three categories: (i) three-dimensional solid elements, (ii) two-dimensional shellelements and (iii) three-dimensional shell elements.

To correctly capture the bending response, when using 3D solid elements, severalelements are required through-the-thickness of the laminate, typically one per ply.Furthermore, due to their propensity to shear-locking, their aspect-ratio must bekept close to unity. From these two considerations, it follows that models relyingonly on 3D solid elements are often computationally unaffordable.

Two-dimensional shell finite elements, based on Equivalent Single Layer Theories(ELST) [93, 94], represent a computationally more convenient alternative for model-ling thin-walled structures. However, the simplifying assumptions of ELST are, often,too restrictive for the analysis of laminated composites (for instance, the continuityassumption made for the displacement field and its derivatives leads to continuousout-of-plane shear strains). To overcome this limitation, 2D shell elements basedon Layer-wise theories [95] might be considered. For the latter, a piecewise con-tinuous through-the-thickness displacement field is assumed and, therefore, strain-discontinuity at the plies interfaces can be modelled.

An intermediate approach between 2D shell elements and 3D solid elements isrepresented by the degenerated three-dimensional shell elements whose formulationis based on two-dimensional kinematic constraints [96, 97]. This type of elementsrepresents a compromise between the computational efficiency of 2D shell elementsand the modelling flexibility of 3D solid elements. Moreover, being shear-locking free,the aspect-ratio limitations of the 3D solid elements are overcome. Nonetheless, theunderlying two-dimensional formulation does not provide an accurate description ofthe out-of-plane displacement and interlaminar stress fields, critical for delaminationevaluations.

4.1.3 Multi-dimensional finite elements coupling

From the the previous section, it follows that 3D solid elements should be used fordiscretizing the area of the structure, denoted as local subdomain, where significantthree-dimensional stress fields are likely to occur, e.g. indentation/impact locationsand where interlaminar damage ought to be modelled. For the remaining portion of

92

Multiple length/time-scale simulation of localized damage in composite structuresusing a Mesh Superposition Technique

the structure, referred to as global subdomain (generally the largest portion), either3D or 2D shell elements are preferable.

Several local/global approaches for coupling differently-discretized subdomainshave been proposed, such as mixed-dimensional coupling approaches based on multi-point constraint (MPC) [89]. The MPC equations can be obtained by equating thework done by the local and global subdomains at the shared interface; a perturbationsolution can then be exploited to determine the stress distribution across the thicknessof the lower-dimension elements [89]. Compared to the Shell-to-Solid coupling optionavailable in the finite element software Abaqus [98], the method is shown to elimin-ate the undesirable stress disturbances at the interface between differently-discretizedsubdomains in the case of isotropic elastic material. However, the local/global coup-ling must be carried at a sufficient distance from boundaries and/or discontinuitiesin the model [99]. In addition, the accurate derivation of the MPC equations forcomposite materials may become impracticable [100].

An alternative approach hinges on the use of ad-hoc transition elements as sug-gested by Davila [90]. In this work, two transition elements have been formulated,i.e. one based on the Mindlin-Reissner [101,102] kinematic assumptions and a secondone based on a higher-order theory similar to that developed by Tessler [103]. Thestiffness matrices KTr of the transition elements are constructed from the stiffnessmatrices K of traditional three-dimensional elements and a restraint matrix R. Res-ults show that, using the lower-order transition elements, it is possible to accuratelymodel the stress fields within both the local and global models. However, withinthe transition elements themselves, a stress boundary layer where the interlaminarnormal and shear stresses are largely overestimated, is obtained. Transition elementsbased on an higher-order shell theory strongly mitigate such spurious stresses [90].

Alternatively, the global and local subdomains can be coupled using an uncoupledglobal/local approach where the displacement fields computed using the global modelproduce the boundary conditions for the local one. The mechanical responses of thelocal and global models are simulated separately and the sequence of global/localanalyses can be either run once (see, as an example, the Submodelling procedureimplemented in Abaqus [85]) or iteratively, until the force/momentum convergence isreached at the global/local interface.

Reinoso et al. compared the Submodelling technique to Shell-to-Solid coupling[104,105]. Here, the two approaches have been applied to the global/local FE analysisof debonding failure at the skin-stringer joint within an aeronautical component. Theanalysis was limited to the first stages of damage propagation, and the size of thelocal models was chosen to be sufficiently larger than the final expected damaged

93

Chapter 4

zone. As the applied load and therefore the damage area grow, stress disturbanceswere observed at the local/global boundaries for both approaches.

The computational advantage provided by generic global/local approaches is en-hanced if the size of the three-dimensional local model is kept to a minimum. However,the transition between subdomains discretized using finite elements with different di-mension/formulation should be located such that the stress field in the local subdo-main has decayed to the form assumed in the formulation of the finite elements usedin the global subdomain. Further, the global/local transition should be sufficientlydistant from any perturbation such as boundaries or damaged zones.

The influence of the distance of the global/local transition from the delamin-ation front was investigated by Krueger et al. for the cases of delaminated testspecimens [106, 107] as well as skin/stringer debonding in a composite aircraft com-ponent [108,109]. Results, provided in terms of the mixed-mode strain energy releaserates computed using the VCCT, demonstrate that, for an acceptable agreement withthe reference solution, a minimum length of which the local subdomain needs to beextended, both ahead and in the wake of delamination front, can be identified.

Bridging methods for coupling continuum models with atomistic models [110,111],as well as discrete element models and finite element models [112–114], have beenpresented in the literature. In the approach described by Belytschko and Xiao [110,111], coupled subdomains overlap at the shared interface as in the method presentedhere.

In this chapter, a new local/global coupling approach based on a Mesh Super-position Technique (MST) is proposed. The interfaces between subdomains whoseFE discretizations consist of different element types are replaced by transition regionswhere the corresponding discretizations are superposed. The theoretical details of theMST, as well as the key aspects of its FE implementation, are presented in § 4.2. TheMST is applied to the multiple length/time-scales analysis of a low-velocity impacton a composite plate [115, 116] in § 4.3. Results of this analysis are presented anddiscussed, respectively, in § 4.4 and § 4.5. Finally, conclusions are drawn in § 4.6.

4.2 Mesh Superposition Technique

4.2.1 Theory

In the literature, techniques involving superposed elements occupying the same phys-ical space have been used to address problems in which meshing is challenging, such

94


as in meso-scale modelling of woven composites [117–120]. In the MST we proposehere, we also make use of superposed elements, albeit in a more generic way.

Let us consider a deformable body B occupying the domain ω(t) ∈ R3, where

with R3 we denote the three-dimensional euclidean space, for any instant of time

t > 0 (Figure 4.1). The domain occupied by the deformable body at instant t = 0is indicated with Ω and, the position vector of each point P ∈ B in the referenceconfiguration is denoted as X.

bA

A

s s

B

bB

u

¡u

X

P

¡tT

N

¡

e2

e3

e1

bA

bB

s=

S

Figure 4.1: Mesh Superposition Technique. In the reference configuration, domainΩ is decomposed into the subdomains ΩA and ΩB which overlap overthe subdomain Ωs. Dirichlet and von Neumann boundary conditionsare applied over specific portions of the boundary Γ=∂Ω.

Let the boundary of Ω be Γ=∂Ω, and assume that Dirichlet boundary conditionsare applied over the subset Γu ⊂ Γ, while von Neumann boundary conditions areapplied over the subset Γt ⊂Γ.

We consider domain Ω to be subdivided into subdomains ΩA and ΩB, whoseintersection is Ωs ≡ ΩA ∩ ΩB �= ∅. Additionally, the subdomains ΩA and ΩB can beconveniently decomposed, as shown in Figure 4.1, into an overlapping portion Ωs anda non-overlapping one defined as

Ωi =Ωi\Ωs for i=A, B. (4.1)

Let ρ(X) and ρ(X)u(X,t) be, respectively, the material mass density and theinertia forces per unit volume. The applied body force and the Piola traction vectorare denoted by B(X,t) and T(X,t) respectively. Given the applied Dirichlet and von

95

Chapter 4

Neumann boundary conditions, the principle of virtual work can be expressed as∫Ω

S :δE dΩ +∫Ω

ρu·δu dΩ=∫Ω

B·δudΩ +∫Γt

T·δu dΓ , (4.2)

where δE(X,t) is the virtual variation of the Green-Lagrange strain tensor, and S(X,t)is the second Piola-Kirchhoff stress tensor.

Within subdomain Ωs, where univocal definitions of the second Piola-Kirchhoffstress tensor S(X,t) and of the material density ρ(X) do not exist, the terms of theprinciple of virtual work can be expressed as the linear combination of the contribu-tions of the overlapping subdomains. These contributions are scaled by non-negativescalar-valued weighting factors ψi(X), continuous over the subdomain Ωs.

The internal virtual work δWint, the external virtual work δWext and the inertialvirtual work δWkin can, therefore, be expressed as:

δWint(u,δu) =∫Ωs

∑i∈{A,B}

[ψi(X)Si

]:δE dΩ +

∑i∈{A,B}

[ ∫Ωi

Si :δE dΩ]

,

δWext(u,δu) =∫Ωs

∑i∈{A,B}

[ψi(X)Bi

]·δu dΩ +

∑i∈{A,B}

[ ∫Ωi

Bi · δu dΩ]

+∫Γt

T·δu dΓ ,

δWkin(u,δu) =∫Ωs

∑i∈{A,B}

[ψi(X)ρi

]u·δu dΩ +

∑i∈{A,B}

[ ∫Ωi

ρiu·δu dΩ]

,

(4.3)with Si and ρi being the second Piola-Kirchhoff stress tensors and the material massdensities of subdomains Ωi; the applied body forces to subdomains Ωi are denotedby Bi.

From Equation 4.3 it follows that, within the subdomain Ωs, an equivalent secondPiola-Kirchhoff stress tensor Ss and and equivalent material mass density ρs can bedefined as

Ss(X,t)=∑

i∈{A,B}ψi(X)Si(X,t) and ρs(X) =

∑i∈{A,B}

ψi(X)ρi(X). (4.4)

Similar considerations hold for the applied body forces Bi. To satisfy the con-servation of energy principle, the collection of weighting factors ψi(X) must be apartition of unity of subdomain Ωs, i.e.

∑i∈{A,B}

ψi(X)=1 ∀ X∈Ωs. (4.5)

96


4.2.1.1 Weighting factors computation

The weighting factors ψA(X) and ψB(X) are computed for every point P (X) :X∈Ωs

as function of the relative distance from two generic surfaces SA and SB, see Figure 4.2,that do not necessarily coincide with ΓA and ΓB, i.e. the interfaces shared by Ωs withΩA and ΩB, respectively.

e2

e3

e1

b¡A

b¡B

SA

SB

¡

P

A

B

bA

s

bB

d a

b

dA

dB

A

Figure 4.2: Distances used for the computation of the weighting factors ψA and ψB.

The procedure to compute ψA(X) and ψB(X) consists of the following steps:

(i) Using the Level Set Method [121], evaluate the signed distances dA and dB ofpoint P (X) from surfaces SA and SB. Letting the closest points to P (X) onSA and SB be denoted, respectively, as A(XA) and B(XB) (see Figure 4.2), dA

and dB can be expressed as

dA =‖XA−X‖ and dB =‖XB −X‖ . (4.6)

(ii) Compute distances d, a and b, shown in Figure 4.2, as

d=‖XB −XA‖ , a=|(XB −X) · (XB −XA)|

d, b=

|(XA−X) · (XB −XA)|d

(4.7)

(iii) Assign weighting factors ψA(X) and ψB(X) as:

ψA(X)=

⎧⎪⎪⎨⎪⎪⎩0 ⇐ dA >da

d⇐ dA, dB <d

1 ⇐ dB ≥d

and ψB(X)=

⎧⎪⎪⎪⎨⎪⎪⎪⎩0 ⇐ dB >db

d⇐ dA, dB <d

1 ⇐ dA ≥d

(4.8)

97

Chapter 4

From the definitions provided in Equations 4.6-4.8, we infer that ψA(X) andψB(X) are continuous over Ωs, and that the Partition of Unity condition (see Equa-tion 4.5) is automatically satisfied.

4.2.2 Finite Element implementation

To derive the finite element equations associated to the problem formulated in Equa-tion 4.2, let subdomains ΩA and ΩB be discretized with elements ej

A and ejB , respect-

ively, and let the corresponding integration domains be Ωej

Aand Ω

ejB

, such that

⋃ej

A ∈EA

Ωej

A≡ ΩA and

⋃ej

B ∈EB

Ωej

B≡ ΩB , (4.9)

where EA ={

ejA

}and EB =

{ej

B

}indicate the sets of elements with support in

subdomains ΩA and ΩB as shown in Figure 4.3.

bEA

bEB

EsA

=

MST regionEsB

S SS

S=

Figure 4.3: FE meshes of the superposed subdomains ΩA and ΩB. Within the MSTregion, the two meshes are superposed and the corresponding stiffnessand mass matrices scaled, in order to satisfy the conservation of energyprinciple.

Proceeding as in § 4.2.1, EA and EB can be further decomposed as EA = EA∪EsA

and EB = EB ∪ EsB, where EA and EB are the sets of elements whose support is,

respectively, in subdomains ΩA and ΩB. In addition, EsA ⊂ EA and Es

B ⊂ EB denote

98


the sets of elements which of overlap over the transition region (indicated in Figure 4.3as MST region), i.e.

⋃ej

A ∈EA

Ωej

A≡ ΩA ,

⋃ej

B ∈EB

Ωej

B≡ ΩB and

⋃ej

A ∈EsA

Ωej

A≡⋃

ejB ∈Es

B

Ωej

B≡ Ωs. (4.10)

From Equation 4.3, the stiffness matrices of elements within the MST region, i.e.ej

A ∈EsA and ej

B ∈EsB, are given by

Ksej

A=∫Ω

ejA

[BT

ejA

: ψA(X)(I ⊗ SA + FCAFT

): B

ejA

]dΩ

ejA

and

Ksej

B=∫Ω

ejB

[BT

ejB

: ψB(X)(I ⊗ SB + FCBFT

): B

ejB

]dΩ

ejB

.(4.11)

where Bej

Aand B

ejB

are the shape functions derivatives matrices of finite ele-ments ej

A and ejB , respectively. In addition, Ci is the fourth-order material elasticity

tensor associated to subdomains Ωi, F is the deformation gradient and I identifiesthe identity matrix.

The values of ψA(X) and ψB(X) can be assigned to each element through suitablenodal values, and the element’s shape functions may be used for their evaluation atevery internal point. Alternatively, if the variation of ψA(X) and ψB(X) inside eachelement is neglected, then the centroid value can be used to scale the constitutiveproperties uniformly for the entire element. Since ψA(X) and ψB(X) are computedas function of the elements’ coordinates in the undeformed configuration, their valuesdo not change as the body deforms.

Following a similar approach to the one used to define the stiffness matrices inEquation 4.11, the mass matrices of elements ej

A ∈ EA and ejB ∈ EB are computed as

Msej

B=∫Ω

ejA

ψA(X)ρANTej

AN

ejA

dΩej

Aand

Msej

B=∫Ω

ejB

ψB(X)ρBNTej

BN

ejB

dΩej

B.

(4.12)

99

Chapter 4

4.3 Application: Multiple length/time-scale analysis oflow-velocity impact on a composite plate

4.3.1 Problem description

The MST is applied for simulating the low-velocity impact response of [0◦3/90◦

3]scross-ply laminates made of HS160/REM graphite/epoxy [115, 116]. The specimenswere rectangular with a 65.0 mm × 87.5 mm in-plane area and a nominal thicknessh of 2.0 mm. The specimens were impacted with a 2.3 kg hemispherical impactor12.5 mm in diameter and were simply supported on a steel plate with a rectangular(45.0 mm × 67.5 mm) opening.

An impact energy of 3.1 J was used in the experiments [115, 116]; the impact-induced damage consists of an initial tensile matrix crack in the distal 0◦ layers,followed by a two-lobe shaped delamination at the bottom 0◦/90◦ interface as shownin Figure 4.4.

Tensile matrix crack

Two-lobe shaped

delamination

11.75 mm

11.75 mm

0°

Figure 4.4: Impact-induced damage pattern on [0◦3/90◦

3]s laminate. X-rays analysisshows that the main failure modes are represented by a tensile matrixcracking of the distal 0◦

3 sublaminate (green) and by a two-lobe shapeddelamination at the bottom 0◦/90◦ interface (red), after [116].

100


4.3.2 Finite Element models

Composite structures subjected to low-velocity impact of small objects experiencelocalized damage confined within the area surrounding the impact location [122].Therefore, to accurately simulate such localized damage, while maintaining the re-quired computational cost to a minimum, different parts of the structure need to bemodelled at different length-scales. More precisely, the area surrounding the impactlocation can be conveniently modelled at the meso-scale while the remaining of thestructure being modelled at the macro-scale [123]. In the following, we will refer tothe area modelled at the meso-scale, either as meso-scale region or local region, asopposed to the remaining area, referred to as the macro-scale region or global region.The in-plane area of the meso-scale region will be denoted as Ameso = Lmeso×Wmeso

while, Atot =L×W is the panel’s in-plane area, see Figure 4.5a. In the work presentedhere, the panel’s in-plane area Atot was kept constant, while the width Wmeso andthe length Lmeso of the meso-scale region were varied.

R

IMPACTOR

h

W=2

L=2

x

y

z

x

0±

N

Bottomface

yTopface

J²z

x-symmetryplane

y-symmetryplane

0±3

90±6

90±

Bottom 0±=90±

interface

Distal 0±3sublaminate

0±3

Lmeso=2

Wmeso=2

Meso-scaleregion

(a) Schematic (not in scale) with layup anddimensions of plate and impactor.

xz

y

0±00±90±

`ehcohe

3D solidelements

Cohesiveelements

Topface

(b) FL model.

Figure 4.5: Schematic of the impacted specimens with definition of the meso-scaleregion (a) and the FL model used to obtain the reference solution (b).

Within the meso-scale region, the differently-oriented sublaminates are individu-ally modelled, as well as the interfaces between them. To correctly capture the three-dimensional stress field at the impact location, 3D solid elements are used in this

101

Chapter 4

region. The remaining of the structure, modelled at the macro-scale, is convenientlydiscretized using either 3D or conventional 2D shell elements.

FE models of the impacted composite plates were built using the finite elementssoftware Abaqus (6.12-1) [85], through the in-built Python scripting facilities [124].Because of the symmetry, only one quarter of the specimen was modelled and onlythe portion of the specimen within the rectangular opening was considered, as schem-atically shown in Figure 4.5a. The impactor’s deformations were assumed negli-gibly small and, therefore, it was modelled using rigid surface elements (R3D4); theimpactor-plate contact interaction was simulated by surface-to-surface contact ele-ments between the impactor surface and the top surface of the plate.

The reference solution of the problem in § 4.3.1 was obtained with a fully localmodel, denoted as FL model and shown in Figure 4.5b. Here, the plies are uniformlydiscretized using eight-noded reduced-integration solid elements (C3D8R) and co-hesive elements were inserted at the bottom 0◦/90◦ interface over the entire plate’sin-plane area Atot.

The MST was applied for both a 3D/3D coupling and a 3D/2D coupling. Intotal, four multi-scale models were built: (i) two models with a sudden discretization-transition between the meso-scale and the macro-scale region, denoted as 3D/3D ST(Figure 4.6a) and 3D/2D ST (Figure 4.7a) and (ii) two models where the suddendiscretization-transition of the ST models is replaced by a transition region overwhich the discretizations are superposed, and referred to as 3D/3D MST (Figure 4.6b)and 3D/2D MST (Figure 4.7b). A constant width of four elements in the laminate’splane was considered for the MST region.

To model interlaminar failure, a layer of cohesive elements (COH3D8) was insertedat the bottom 0◦/90◦ interface while two rows of cohesive elements, placed at thesymmetry plane parallel to the 0◦-plies direction, were used to simulate tensile matrixcracking (see Figures 4.5b-4.7b). A mixed-mode bilinear traction-separation law wasused to simulate the softening and fracture response [125, 126]. A quadratic stress-based criterion and a linear energy-based criterion were assumed, respectively, fordamage initiation and propagation.

The material properties used by Aymerich et al. [115] in their simulations areshown in Table 4.1. Although they did not use the typical assumption of transverseisotropy, we used the same properties in our simulations so that meaningful compar-ison can be drawn. The cohesive properties of the interfaces were calibrated againstexperimental data from static Mode I (DCB) and Mode II (ENF) fracture tests onunidirectional laminates [115].

102


3D solidelements

xyz

0±00±90± 3D shell

elements

Topface

(a) 3D/3D ST model.

MST

0±00±90±

xyz

3D shellelements

3D solidelements

Topface

(b) 3D/3D MST model.

Figure 4.6: ST and MST models used for the 3D/3D meso/macro coupling. 3D solidelements and 3D shell elements are displayed, respectively, in blue andgreen, while the cohesive elements are shown in red.

Topface

xyz

0±00±90±

3D solidelements

2D shellelements

(a) 3D/2D ST model.

MST

xyz

0±00±90±

3D solidelements

2D shellelements

Topface

(b) 3D/2D MST model.

Figure 4.7: ST and MST models used for the 3D/2D meso/macro coupling. 3D solidelements and conventional 2D shell elements are displayed, respectively,in blue and green, while the cohesive elements are shown in red.

4.3.2.1 Multiple length-scale analysis

Within the meso-scale region (see Figure 4.5a), composite plies were discretized us-ing eight-noded reduced-integration solid elements (C3D8R) and delamination initi-

103

Chapter 4

Table 4.1: Elastic properties of the plies, cohesive properties of the interfaces, anddensity used by Aymerich et al. [115].

E1 E2 =E33 G12 = G13 = G23 ν12 = ν13 = ν23 ρ

93.7 GPa 7.45 GPa 3.97 GPa 0.261 1600 kg/m3

kN N GIc120 GPa/mm 30 MPa 520 J/m2

kS = kT S = T GIIc = GIIIc43 GPa/mm 80 MPa 970 J/m2

ation/propagation at the bottom 0◦/90◦ interface was taken into account using co-hesive elements. Within the macro-scale region, for the 3D/3D models, eight-nodedreduced integration 3D shell elements (SC8R) were used; the latter were replaced byconventional four-noded reduced-integration 2D shell elements (S4R) for the 3D/2Dmodels. Both 3D and 2D shell elements were not integrated though-the-thicknessduring the analysis; therefore, no additional integration points were needed.

For all models, a constant in-plane element length e = 0.25 mm was adopted forthe entire structure and, a thickness hcoh

e =20 μm was used for the cohesive elements.Each of the 0◦

3 and 90◦3 sublaminates were discretized with two elements through their

thickness when three-dimensional (solid or shell) elements were used.

Simply supported boundary conditions were enforced either by constraining thez-displacement of the nodes initially lying on the edges of rectangular supportingopening (3D/3D models), or through specific MPC equations prescribed at the shellelements’ nodes on the outer edges of the panel (3D/2D models). Further MPCequations were prescribed at the shell-to-solid interfaces including, for the 3D/2DMST model, the nodes within the MST regions.

Simulations were performed using Abaqus/Explicit with a constant time stepΔt = 10 ns (Δtstable ≈ 16.22 ns) and the total analysis time was equal to 4 ms. Nomass scaling was considered during the simulations and enhanced hourglass controlwas used for reduced-integration elements.

4.3.2.2 Multiple length and time-scale analysis

In addition to the multiple length-scale approach introduced in § 4.3.2.1, a multi-solver technique can be exploited. In the latter, different solvers can be used tosimulate the mechanical response of different portions of the structure, dependingon where they are expected to provide the most computationally efficient solution.Considering the problem in § 4.3.1, material failure and complex contact interactionsat the impact location are best analysed using FE solvers based on explicit time-

104


Interface nodes Implicit sub-model

Implicit sub-model

Explicit sub-model

Interface nodes Explicit sub-model

Meso/Macro

Meso-scale region

Macro-scale region

(ST)

(a) ST multi time/length-scale model.

Implicit sub-model

Interface nodes Explicit sub-model

Meso/Macro

Interface nodes Implicit sub-model

Explicit sub-model

Meso-scale region

Macro-scale region

(MST)

(b) MST multi time/length-scale model.

Figure 4.8: Multiple length/time-scale models. The 3D/2D ST and MST models aredecomposed into an explicit and an implicit sub-model which interactthrough the interface nodes displayed in red.

integration schemes, e.g. Abaqus/Explicit, while the elastic behaviour of light andstiff components can, more efficiently, be simulated with FE solvers using implicittime-integration schemes, e.g. Abaqus/Standard.

The use of different solvers leads to the definition of multiple time-scales at whichthe structural response is simulated, i.e. an explicit/micro time-scale and an impli-cit/macro time-scale. The former is characterized by a high number of short andrelatively inexpensive time-steps, while a reduced number of larger time-steps arerequired when implicit integrators are used, due to their unconditional stability.

Within this framework, the 3D/2D ST and MST models were decomposed into anexplicit and an implicit sub-model as shown in Figure 4.8. The 3D/2D coupling,either using a sudden discretization-transition (Figure 4.8a) or the proposed MST(Figure 4.8b), is carried out within the explicit sub-model.

Simulations were performed using Abaqus/Explicit for the explicit sub-model witha constant time-step ΔtXpl = 10 ns (Δtstable ≈ 16.22 ns) while Abaqus/Standard wasused for the implicit sub-model, with a constant time-step of ΔtImp =1 ms; the totalanalysis time was also equal to 4 ms. Enhanced hourglass control was considered forreduced-integration elements and no mass scaling was used within Abaqus/Explicit.

105

Chapter 4

The explicit and implicit analyses are coupled using the GC method [127, 128],a staggered method built within the framework of the FETI method [129, 130] andimplemented in Abaqus [85] in the form of a Co-Simulation Engine (CSE). The velo-city continuity is prescribed at the interface nodes (displayed in red in Figure 4.8) bymeans of Lagrange multipliers. Although the detailed description of the GC methodis beyond the purpose of this chapter, it is worth noting that its computational effi-ciency is primarily dependent on the number of degrees of freedom (DOFs) belongingto the explicit sub-model and on the number of interface nodes.

4.4 Results

4.4.1 Damage prediction MST

The interlaminar damage pattern at the bottom 0◦/90◦ interface, simulated with theFE models detailed in § 4.3.2.1, and using Ameso = 28.25 × 6.00 mm2, is shown fort=2.0 ms and t=4.0 ms, in Figure 4.9a and Figure 4.9b, respectively.

4.4.2 Computational efficiency

4.4.2.1 Multiple length-scale analysis

The CPU time t associated to the 3D/2D ST and to the 3D/2D MST models as function

of the normalized meso-scale areaAmeso

Atot, is provided in Figure 4.10a. Results are

normalized with respect to the CPU time of the FL model t[FL].

Starting from the reference configuration (FL model), the extension of the meso-scale region was progressively reduced until interlaminar damage was observed at themeso/macro-scale interface. The minimum (normalized) meso-scale areas required tocorrectly replicate the damage pattern when using either the ST or the MST model, de-noted in Figure 4.10a respectively as A

[ST]meso and A

[MST]meso, are compared in Figure 4.10b,

together with the corresponding (normalized) CPU times t[ST]min and t

[MST]min . Addition-

ally, in Figure 4.10b we evaluated the benefits of the MST through the reductionof meso-scale area ΔAmeso and the resulting reduction of CPU time Δtmin achievedwhen using the MST instead of a sudden discretization-transition as:

ΔAmeso =A

[MST]meso − A

[ST]meso

A[ST]meso

and Δtmin =t[MST]min − t

[ST]min

t[ST]min

. (4.13)

106


0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

d

FL

3D/3D ST

3D/2D ST

x1=Wmeso

3D/3D MST

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

d

x

Wmeso

x

W3D/2D MST

ST STMST MST

3D/3D coupling 3D/2D coupling

FL

0± 0± 0± 0± 0±

x1x1

(a) t=2.0 ms.

0

0.2

0.4

0.6

0.8

1

0 0.2 0.4 0.6 0.8 1

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

d

d

FL

3D/3D ST

3D/2D ST

3D/3D MST

3D/2D MSTWmeso

ST STMST MST

3D/3D coupling 3D/2D coupling

FL

0± 0± 0± 0± 0±

xx

x1x1

x1=Wmeso

(b) t=4.0 ms.

Figure 4.9: Interlaminar damage evolution at the bottom 0◦/90◦ interface. On theleft, the evolution of the damage variable d as function of the hori-zontal distance from the impact location; the x1-coordinate correspondsto the height where the delamination attains its maximum extension(3.375 mm for (a) and 4.875 mm for (b)). On the right, the delamina-tion pattern for the five FE models considered.

4.4.2.2 Multiple length/time-scale analysis

The ST and MST models used for the multiple length/time-scale analysis have, re-spectively, a meso-scale area equal to A

[ST]meso and to A

[MST]meso, i.e. the smallest meso-scale

areas for which damage is correctly replicated with the sudden-discretization trans-

107

Chapter 4

0%

20%

40%

60%

80%

100%

0% 20% 40% 60% 80% 100%%

Ameso=Atot

bA[MST]meso =22% bA[ST]

meso=56%

bt[MST]min =49%

bt[ST]min =64%bt

60

ST

MST

FLt=t [FL]

(a) Normalized CPU time as function of thenormalized meso-scale area for the 3D/2DST and MST models. The extension ofthe meso-scale region is progressively re-duced and the data point indicated withthe red cross represents the first config-uration for which artificial damage at thediscretization-transition is observed.

1 2

Ameso=Atot

Ameso=Atot

Minimum Minimum t=t [FL]

bA[MST]meso

bA[ST]meso bt [MST]min

bt [ST]min

0%

20%

40%

60%

80%

100%

\ 0%

20%

40%

60%

80%

100%

t=t [FL]

¢Ameso¼¡60% ¢tmin¼¡23%

(b) On the left, the comparison between thesmallest meso-scale areas required to cor-rectly capture the damage pattern, withthe ST model and the MST model; on theright, the comparison between the corres-ponding CPU times required by the op-timal ST and the optimal MST model.

Figure 4.10: Computational efficiency of the MST. The MST models allow the use ofsmaller meso-scale areas and, thus, lower CPU times. For the proposedexample, using the MST the meso-scale area is reduced by nearly 60%,while the CPU time by nearly 23%.

ition approach and using the proposed MST, respectively. The number of DOFs inthe explicit sub-model and the number of interface nodes for the ST model and forthe MST model are compared, respectively, in Figures 4.11a and 4.11b.

In Figure 4.11c, the absolute CPU time reductions |Δtmin| achievable with theMST are compared for the cases of multiple length-scale analysis and of multiplelength/time-scale analysis.

4.5 Discussion

4.5.1 Multiple length-scale analysis

The stress field disturbances resulting from the local/global coupling can lead tounrealistic interlaminar failure at the meso/macro-scale interface, particularly whena sudden discretization-transition approach is adopted (ST models), see Figure 4.9.

108


0.E+00

1.E+05

2.E+05

3.E+05

1 2

# DOFs in explicit subdomain

¼¡43%

ST MST

(a) Number of DOFs in theexplicit sub-model.

0

50

100

150

200

1 2

¼¡20%

ST MST

# Interface nodes

(b) Number of interfacenodes.

0%

20%

40%

60%

80%

100%

1 2No Co-sim.

With Co-sim.

¼23%

¼49%

j¢tminjCPU time reduction

(c) CPU time reductionwhen using the MST.

Figure 4.11: Computational efficiency of the MST model coupled with an impli-cit/explicit co-simulation technique.

From the evolution of the damage variable d at t=2.0 ms (Figure 4.9a), it followsthat, the stress disturbances are more pronounced in the case of 3D/2D couplingrather than for the 3D/3D case. Moreover, when compared to the 3D/3D models, theextension of the lobe-shaped delamination is slightly underestimated when using the3D/2D models (Figure 4.9a). This characteristic is probably related to the differentbending response of 3D solid elements and 2D shell elements.

If the discretization-transition is sufficiently close to the impact location, as thetwo lobe shaped delamination grows, the latter interacts with the artificial failureat the discretization-transition, resulting in excessively large delamination areas (seeFigure 4.9b). The delaminated area predicted with the 3D/3D ST model is muchlarger than that obtained using the 3D/2D ST model.

The CPU time associated to the ST and MST models decreases linearly with theextension of the meso-scale region (see Figure 4.10a) and, as a result of the highernumber of elements and of MPC equations required at the global/local coupling, forequal values of Ameso the CPU time associated to the ST models is always lower thanthat of the MST models. However, as the distance of the discretization-transitionfrom the impact location decreases, the damage pattern is not correctly replicatedand failure at the global/local interface is simulated when using the ST models. InFigure 4.10a, the first configuration for which the delamination pattern is not correctlypredicted with the ST models is marked with a red cross. With the proposed MST,

109

Chapter 4

it is possible to use smaller meso-scale regions (ΔAmeso ≈ −60%) while accuratelymodelling the damage pattern, as well as lower CPU times (Δtmin ≈−23%) as shownin Figure 4.10b.

Although these specific values are, admittedly, dependent on the model’s sizeand on the specific analysis performed, they confirm the capabilities of the approachproposed in this chapter. The apparently large discrepancy between the reductionof the meso-scale area and the corresponding computational gain (60% vs. 23%) isjustifiable by the higher number of finite elements and of MPC equations requiredwhen using the MST.

4.5.2 Multiple length/time-scale analysis

The GC co-simulation technique has been successfully applied to study the responseof reinforced concrete structures subjected blast [131] and earthquake loading [132].Additional applications include the simulation of tire/road interaction during fullvehicle durability tests [133,134]. Recently, Brun et al. [135] analysed a flat compositestiffened panel subjected to localised loads. In Brun et al.’s work [135], damage was notmodelled and the time-scale transition was not combined with an effective couplingbetween different length scales.

To the best of the authors’ knowledge, the example proposed in this chapter repres-ents the first attempt towards the multiple length/time-scale modelling of compositestructures experiencing localized damage.

When the MST is combined with the implicit/explicit co-simulation for a multiplelength/time-scale analysis, the possibility to minimize the meso-scale area becomeseven more attractive. Interestingly, within the framework of a multiple length/time-scale analysis, the meso-scale area reduction provided by the MST, as opposed to asudden discretization-transition approach, allows to decrease of the number of DOFsin the explicit sub-model and of interfaces nodes, by about 43% and 20%, respectively(see Figures 4.11a and 4.11b). As a result, the CPU time reduction achieved whenusing the MST increases, from the 23% obtained for the multiple length-scale analysis,to about 49% in the case of multiple length/time-scale analysis.

4.6 Conclusions

In this chapter, a Mesh Superposition Technique (MST) for a progressive element-type transition between differently-discretized subdomains is proposed.

110


The MST is applied to the multiple length/time-scale analysis of a compositeplate subjected to low-velocity impact. The area of the structure that, in order tocorrectly capture the damage pattern, needs to be modelled at the smallest length-scale can be significantly reduced when compared to a local/global model with asudden discretization-transition (approximately 60% for the proposed example) and,therefore, a computational cost saving might be achieved (approximately 23% CPUtime reduction).

Finally, the MST was coupled with an implicit/explicit co-simulation techniquefor a multiple time/length-scale analysis. The results indicate that, if the length-scaletransition is performed using the proposed MST instead of a sudden discretization-transition, the CPU time can be approximately halved (49%).

111

Chapter 5

Exploiting symmetries insolid-to-shell homogenization,with application to periodicpin-reinforced sandwichstructures

5.1 Introduction

For a faster structural design of large composite components, the numerical analysisof their mechanical response requires different parts of the structure to be modelledat multiple length/time-scales [136]. In addition, due to the heterogeneous micro-structure of composite materials, the applicability of conventional FE analyses basedon 3D (solid) models is limited to small components; 2D (shell) FE models withequivalent homogenised properties are preferable for larger components.

Within this framework, numerical homogenization of periodic structures and ma-terials represents a powerful numerical tool. The latter commonly involves the ana-lysis of Unit Cell (UC) models in which the structure is explicitly resolved at the lowestlength-scale; several studies on the existence and size of suitable UCs are available inliterature [137,138]. Furthermore, numerous works concerning the correct applicationof Periodic Boundary Conditions (PBCs) to UCs can be found, e.g. [139–142] andreferences therein.

113

Chapter 5

However, for several practical cases, e.g. textile composites and pin-reinforcedcomposite sandwich structures, the topological complexity of the representative UCsmay lead to unaffordable modelling/meshing and analysis CPU times. Therefore, theinternal symmetries of the UCs should, whenever they exist, be exploited to reducethe analysis domain, thus enabling a reduction of both the modelling/meshing andanalysis CPU times.

Domains smaller than the UCs, obtained by exploiting the internal symmetries ofthe latter, are referred to as reduced Unit Cells (rUCs). Several works can be founddiscussing the determination of suitable rUCs (and corresponding appropriate PBCs)for UD composites [143], particle-reinforced composites [144] and textile composites[145–147].

The PBCs proposed in these studies allow the determination of the homogenised3D elasticity tensor of the investigated structure. Nevertheless, for the numericalanalysis of large components using equivalent shell models, it is of interest to determ-ine the homogenised shell constitutive response of the structure using high-fidelitythree-dimensional UCs or, preferably, rUCs.

To address this, we present a novel set of PBCs named Multiscale Periodic Bound-ary Conditions (MPBCs), that represents the first set of PBCs that apply to rUCsand enable the direct two-scale (solid-to-shell) numerical homogenization of periodicstructures, including their bending and twisting response.

The proposed MPBCs are formulated in § 5.2.2. Details on their use withinthe context of a solid-to-shell homogenization and on their FE implementation areprovided in § 5.2.5 and § 5.3, respectively. The MPCBs are applied to the solid-to-shell homogenization of a sandwich structure with unequal skins in § 5.4.1 and to theanalysis of the mechanical response of a periodic sandwich structure with unequalskins and pin-reinforced core in § 5.4.2. Conclusions are drawn in § 5.5.

5.2 Theory

5.2.1 Problem formulation

Consider a three-dimensional deformable body B3 occupying the domain D3∈R3, as

shown in Figure 5.1, and assume that one of the dimensions of D3, designated as thethickness (denoted as h), is much smaller than the other two. Furthermore, let D3

be periodic in the plane normal to the thickness direction.

114

Exploiting symmetries in solid-to-shell homogenization, with application to periodicpin-reinforced sandwich structures

h

Equivalent shell model

s

s

@sh

Three-dimensional deformable body B3

D3

Equivalent 2Dconstitutive response

D2

Three-dimensional UC/rUC

2666666664

A11 A12 A16 j B11 B12 B16A22 A26 j B22 B16

sym A66 j sym B66¡ ¡ ¡ ¡ ¡ ¡ ¡

j D11 D12 D16sym j D22 D16

j sym D66

3777777775

Two-scale (solid-to-shell)homogenization

Figure 5.1: Problem formulation. The macroscopic response of the periodic domainD3 can be simulated by means of an equivalent shell model, provided theequivalent 2D constituive response of a representative three-dimensionalUC/rUC (subdomain s) is correctly determined.

For reasons of numerical efficiency, the macroscopic mechanical response of B3 ismore conveniently simulated by means of an equivalent shell model with homogenizedproperties contained in D2∈R

2. The corresponding shell constitutive response can beexpressed by the relation R=Kξ◦, where R is the vector of resultant forces/momentsper unit-length acting on the structure, ξ◦ is the resultant mid-surface membranestrains/curvatures vector and the K matrix (referred to, in lamination theory, as theABD matrix) contains the equivalent shell stiffness terms.

Therefore, the accurate determination of the entire K matrix, including all bend-ing and twisting terms, as well as the shear-extension, bending-extension and bending-twisting coupling terms, is of paramount importance for the use of equivalent shellmodels of large composite components.

115

Chapter 5

The homogenized 2D response of a three-dimensional periodic structure can beevaluated through the analysis of three-dimensional UCs/rUCs modelled at the meso-scale, provided the appropriate PBCs are applied to its boundaries. These PBCs,referred to in the following as Multiscale PBCs (MPBCs), are prescribed such thatthe three-dimensional UC/rUC behaves as it were contained in an infinite body B3

(in domain D3) whose macroscopic response can be analysed using shell theory (indomain D2).

5.2.2 Multiscale Periodic Boundary Conditions (MPBCs)

The formulation of the MPBCs is derived following the equivalence framework pro-posed by de Carvalho et al. [146,147].

5.2.2.1 Physical equivalence and Periodicity

Consider a deformable body B3 occupying the domain D3 ∈R3 and let s, s⊆D3 be

two generic subdomains within D3 as shown in Figure 5.2. Let subdomains s and s

have, respectively, Local Coordinate Systems (LCSs) [O,ei]s and [O,ei]s, where O andO are the origins of the LCSs, while ei and ei indicate the bases defining the LCSs,with i = {1,2,3}. Furthermore, assume a generic spatial distribution of n physicalproperties Πj with j∈{1, . . . , n} and let each of these to be expressed as tensors Πj

s

and Πjs, respectively in [O,ei]s and [O,ei]s.

Definition 1 Two subdomains s and s are physically equivalent, i.e. s∧= s, if

∀ A∈s ⇒ ∃ A∈ s

∣∣∣∣ (OA)s≡(OA)s ∧ Πjs (A)≡Πj

s(A). (5.1)

In Equation 5.1, as well as in the remaining of this work, the subscripts refer tothe LCS in which the vectors and tensors are expressed, when it is required to do so.

Definition 2 A domain D3∈R3 is periodic if it can be reconstructed by tessellation

of non-overlapping physically equivalent (see Equation 5.1) subdomains with parallelLCSs as shown in Figure 5.3, i.e. if

sk∧= sl ∧ [O,ei]sk

‖ [O,ei]sl, ∀ k �= l . (5.2)

The smallest subdomain verifying Equation 5.2 is referred to as Unit Cell (UC).

116


s

s

s

@s

@sD3

(OA)s´(OA)s ^¦js (A)´¦j

s(A)

s

¦js(A)

¦js(A)

^

e2

O

e1

e3

A

e3

e2e1

O

A

Figure 5.2: Physical equivalence.

¯e3

e1¯ee3

ee111111¯Os e2 ¯ssss ¯

e3

e2

e1

reduced Unit Cell (rUC)

Unit Cell (UC)

e2

¯e3e1

Periodic body B3

s

¯

O

s

¯O

D3

@s

A

A

A

A

@s

Figure 5.3: Periodicity condition and Unit Cell (UC).

5.2.3 Subdomain admissibility

Let [O, ei]s and [O, ei]s be, respectively, the LCS of subdomains s and s, as shown inFigure 5.4. Let T={Tij} be the transformation matrix between [ei]s and [ei]s, where

Tii =∂ei

∂eiwith i ∈ {1,2,3} and Tij = 0 if i �= j; furthermore, denote the corresponding

projection in the 1,2 plane as T, i.e. T={Tij} with i, j ∈{1, 2}.

From Figure 5.4, the position of two equivalent points A ∈ s and A ∈ s can berelated by

(OA) = (OO) + T(OA) . (5.3)

117

Chapter 5

e3e3

¯

s Xs

XsXs

s

e3

e2

e1

e2

e1A

A

X[O]ss

O

O

Figure 5.4: Geometrical relation between equivalent points.

In the remaining of this work, we will indicate the three-dimensional displacementvector as u={u, v, w}t and the corresponding projection in the 1,2 plane as u={u, v}t.The displacement gradient at the point A∈s can be decomposed as

∇u(A) = s〈∇u〉 + s∇u(A) , (5.4)

where s〈•〉 =1Vs

∫Vs

(•)dV denotes the volume average operator while, s•=• − s〈•〉

represents the first-order fluctuation term over subdomain s. Equation 5.4 is validonly if the the displacement vector and the gradient operator are expressed in the LCSof the domain where the volume average is taken and the displacement fluctuationterm is evaluated.

Definition 3 A generic subdomain s is defined as admissible for the analysis of aperiodic structure under a given loading s〈∇u〉 if for any other subdomain sk, withtransformation matrix Tk, there is a γk =±1 such that

s〈∇u〉 = γkTk

[sk〈∇u〉

]Tk ∀ k , (5.5)

where γk = ±1 is the loading reversal factor (correspondent to subdomain sk) used toenforce the equivalence between the stress/strain fields of physically equivalent subdo-mains [146].The subdomain, smaller than the UC, verifying the condition expressedin Equation 5.5 is referred to as reduced Unit Cell (rUC), see Figure 5.3.

118


5.2.4 Derivation of the MPBCs

Within subdomain s, the displacement field u(A) can be expressed in the form of atruncated Taylor series as

u(A) = u(O) + s〈∇u〉(OA) + s ˜u(A) , (5.6)

where u(O) is the displacement at the origin of the LCS of s while s˜• indicates thehigher-order (2nd and higher) fluctuation components of the field • over subdomains.

Let s and s be two adjacent subdomains satisfying both the physical and loadingequivalence (Equations 5.1 and 5.2) and consider two equivalent points A ∈ s andA ∈ s. If A is at the boundary of s, i.e. A ∈ ∂s, then also A ∈ ∂s (see Figure 5.3);therefore, according to Equation 5.6, the displacement at these two points is

u(A)=u(O) + s〈∇u〉(OA) + s˜u(A) , (5.7)

u(A)=u(O) + s〈∇u〉(OA) + s˜u(A) . (5.8)

The fluctuation components of the displacement field at two equivalent points arerelated, as derived by de Carvalho et al. [146,147], by the expression

s˜u(A) = γT[s˜u(A)

]. (5.9)

Thus, if Equation 5.8 is pre-multiplied by γT and subtracted to Equation 5.7,we obtain the equation to apply PBCs for the analysis of a generic subdomain (thesuperscript s for the volume average operator is omitted for readability):

u(A)−γTu(A)=[I−γT]u(O)−〈∇u〉T(OO)+[〈∇u〉T−γT〈∇u〉](OA) , (5.10)

where I is the 3×3 identity matrix. Therefore, Equation 5.10 is applied to all pointsbelonging to the boundary of subdomain s and the external loading can be appliedby prescribing the terms of 〈∇u〉. At this point, let one of the dimensions of D3,designated as the thickness and denoted as h, be much smaller than the other two, andassume that D3 is periodic in the plane normal to the thickness direction. Under theseassumptions, the thickness of subdomain s∈D3 is also equal to h. The displacementfield at the point A ∈ s, i.e. u(A) = {u(A), v(A), w(A)}t, can be decomposed intoan in-plane component u(A)={u(A), v(A)}t and an out-of-plane displacement term

119

Chapter 5

w(A), i.e.

u(A) − s ˜u(A) = u◦(A)−xA3 ∇w◦(A) ∀ xA

3 ∈[−h

2,h

2

], (5.11)

w(A) − s ˜w(A) = w◦(A) , (5.12)

where xA3 is the out-of-plane coordinate of point A and the subscript ◦ denotes the

quantities evaluated at the mid-surface of subdomain s. Within such framework,Equation 5.10 can be reformulated for the in-plane and out-of-plane displacement,respectively as

u(A)−γTu(A) =[〈∇u◦〉−

(T33xA

3

)〈∇∇w◦〉

]T[(OA)◦−(OO)◦

]+[I − γT] u(O) +

(5.13)

− γT[s〈∇u◦〉 − xA

3 〈∇∇w◦〉](OA)◦ ∀ xA

3 ∈[−h

2,h

2

],

w(A)−γT33w(A)=[(OA)◦−(OO)◦

]tT[

〈∇∇w◦〉2

]T[(OA)◦−(OO)◦

]+[1−γT33] w(O) +

(5.14)

− γT33(OA)t◦

[〈∇∇w◦〉

2

](OA)◦ ∀ xA

3 ∈[−h

2,h

2

],

where I is the 2×2 identity matrix and xA3 is the out-of-plane coordinate of point A.

The membrane deformation tensor ∇u◦ and the curvature tensor ∇∇w◦ are defined,respectively, as

∇u◦ =

⎡⎢⎢⎢⎢⎣∂u◦∂x1

∂u◦∂x2

∂v◦∂x1

∂v◦∂x2

⎤⎥⎥⎥⎥⎦ and ∇∇w◦ =

⎡⎢⎢⎢⎢⎢⎣∂2w◦∂x2

1

∂2w◦∂x1∂x2

∂2w◦∂x1∂x2

∂2w◦∂x2

2

⎤⎥⎥⎥⎥⎥⎦ . (5.15)

At this point, the following aspects of the proposed MPBCs can be highlighted:

(i) the proposed formulation is valid for both conventional UCs and, more rel-evantly, for rUCs. Secondly, no limitations on the undeformed shape of theUC/rUC are assumed, therefore, complex UC/rUC geometries can be con-sidered;

(ii) no limitations on the deformed shape of the UC/rUC were hypothesized. Hence,the plane-remains-plane boundary conditions [148,149], which have been provento be over-constraining [139], as well as to violate the stress-strain periodicitycondition [140], are avoided;

120


(iii) since no periodicity exists in the thickness direction, the MPBCs are appliedonly on the lateral faces of the UC/rUC; in addition, as a result of the plane-stress assumption of shell theory, the upper and lower faces of the UC/rUC aretraction-free;

(iv) the MPBCs are expressed in terms of the relative displacement between equi-valent points belonging to the boundary of the UC/rUC. This approach is sup-ported by the results presented in [142] which demonstrate that imposing equi-librium and/or compatibility equations as part of the boundary conditions, seee.g. [150], is not only unnecessary but it also violates the minimum total poten-tial energy principle;

(v) no a priori restrictions to the nonlinear behaviour at the lowest length-scalewere made (provided there is no localization [151]);

(vi) The presented MPBCs can be used for the two-scale (3D-2D) homogenizationof any microstructure, regardless of its complexity, and all loading componentscan be prescribed to the UC/rUC, including bending and twisting.

5.2.5 Two-scale (solid-to-shell) homogenisation of periodic struc-tures

Provided the MPBCs are applied to all points on the lateral boundary of the three-dimensional UC/rUC, the 2D constitutive response, i.e. all terms of the equivalentshell stiffness matrix K, can be determined by sequentially subjecting the UC/rUCto the six conventional loadings of shell theory, i.e. membrane stretching, membraneshearing, bending and twisting, see Figure 5.5. These are applied by specifying thecomponents of 〈∇u◦〉 for the macroscopic membrane deformations, and of 〈∇∇w◦〉 forthe bending and twisting curvatures. In the remaining, all quantities are evaluatedin the LCS of the UC/rUC; the corresponding subscripts and superscripts will beomitted hereafter for convenience.

The components of 〈∇u◦〉 and of 〈∇∇w◦〉 can be related to the membrane strainsand curvatures adopted in shell theory as follows:

ε1◦ =:⟨

∂u◦∂x1

⟩, ε2◦ =:

⟨∂v◦∂x2

⟩, γ12◦ =:

⟨∂u◦∂x2

⟩+⟨

∂v◦∂x1

⟩

κ1◦ =:⟨

∂2w◦∂x2

1

⟩, κ2◦ =:

⟨∂2w◦∂x2

2

⟩, κ12◦ =: 2

⟨∂2w◦

∂x1∂x2

⟩ . (5.16)

where ε1◦ and ε2◦ are the membrane direct strains, γ12◦ is the membrane shear strain,κ◦

1 and κ◦2 are the bending curvatures and κ◦

12 is the twisting curvature.

121

Chapter 5

e3

e2

e1

(a) In-plane stretching ε1◦.

e3

e2

e1

(b) In-plane stretching ε2◦.

e3

e2

e1

e3

e2

e1

(c) Shearing γ12◦ .

e3

e2

e1

e3

e2

e1

(d) Bending κ1◦.

e3

e2

e1

e3

e2

e1

(e) Bending κ2◦.

e3

e2

e1

e3

e2

e1

(f) Twisting κ12◦ .

Figure 5.5: The six loading cases in shell theory.

The 2D equivalent strains/curvatures vector ξ◦ =[ε1◦, ε2◦, γ12◦ , κ1◦, κ2◦, κ12◦

]tcanbe related to the vector of the equivalent forces/moments per unit lengthR = [N1, N2, N12, M1, M2, M12]t, through the tensor equation R = Kξ◦ or, in itsexpanded form: ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

N1

N2

N12

−−M1

M2

M12

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦=

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

A11 A12 A16 | B11 B12 B16

A22 A26 | B22 B16

sym A66 | sym B66

− − − − − − −| D11 D12 D16

sym | D22 D16

| sym D66

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣

ε◦1

ε◦2

γ◦12

−−κ◦

1

κ◦2

κ◦12

⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦(5.17)

The K matrix can be determined according to either of the two procedures:

Procedure A- by prescribing a non-zero value to the i-th term of ξ◦ (while fixing theremaining terms to zero), all terms of R are, in the most general case, non-zero.Therefore, from Equation 5.17 the i-th column of the K matrix is proportional

122


to the vector of generalized forces that is obtained from the solution of the FEproblem. Repeating the same procedure for all six components of ε◦, the entireK matrix can be computed;

Procedure B- by prescribing a non-zero value to the i-th term of ξ◦ (while lettingthe remaining terms unconstrained), all terms but the i-th of R are equal tozero. Hence, the i-th column of the compliance matrix C = K−1 (ξ◦= CR) isproportional to the vector of generalized displacements that is obtained fromthe solution of the FE problem. Repeating the same procedure for all sixcomponents of ξ◦, the entire C matrix can be computed, and by inverting thelatter, the K matrix can be determined.

5.3 FE implementation

Concerning the FE implementation of the proposed MPBCs, let us consider a dis-cretized UC/rUC as the one shown in Figure 5.6a. The MPBCs, expressed for thecontinuum case in the form of Equations 5.13 -5.14, are enforced through constraintequations (CEs) acting on equivalent nodes on the boundary of the UC/rUC (see e.g.nodes A and A in Figure 5.6b). Since they are based on a node-to-node coupling,the application of the MPBCs is ideally suited for a periodic discretization along theboundary of the UC/rUC.

The terms of ξ◦ are prescribed by specifying the degrees of freedom (DOFs) ofmaster nodes that do not belong to the mesh of the UC/rUC. In principle, a singlemaster node M with at least six DOFs is sufficient for this purpose. Regarding theterms of the generalized forces vector R, these simply correspond to the reactionforces/moments at node M. Therefore, assuming that the i-th component of ξ◦ isassociated to the j-th DOF of node M, the i-th term of R corresponds to the j-threaction force/moment at node M.

The MPBCs have been implemented in the Finite Element (FE) package Abaqus(6.12) [85] through the in-built Python scripting facilities [98]. The simulations havebeen performed using Abaqus/Standard, and linear geometric behaviour was used inthe examples that follow.

123

Chapter 5

e3

e2

e1

e3

e2

(a) FE discretization of a generic heterogen-eous UC/rUC.

e1

e2

e3

e2

A

A

M

(b) Equivalent nodes at the boundary of theUC/rUC.

Figure 5.6: FE implementation of the MPBCs. The DOFs of equivalent nodes onthe boundary of the UC/rUC are coupled using constraint equations;the external loading are specified through the DOFs of the master nodeM which does not belong to the mesh of the UC/rUC.

5.4 Numerical Examples

5.4.1 Validation

5.4.1.1 Description

The aim of this validation example is to demonstrate that the current formulationcan predict the fully-populated shell stiffness matrix K, in a situation for which theexact analytical solution is known, even if there are no real in-plane symmetries. Todo so, we consider a composite sandwich structure with unequal skins, for which thevalues of the corresponding fully-populated K matrix can be analytically computedby means of Classical Lamination Theory (CLT) [152].

The considered sandwich structure can be ideally subdivided in UCs, as schemat-ically shown in Figure 5.7; however, given the structure’s homogeneity at the ply-level,such subdivision is a purely mathematical abstraction rather than a result of peri-odicity. Furthermore, let the selected UC (see Figure 5.7) be subdivided into rUCswhich, in the most general case, have different LCSs. As for the definition of theUCs, for this specific example, the existence of rUCs results from the definition of

124


their LCSs rather than from the exploitation of UC’s internal symmetries. Nonethe-less, albeit the solid-to-shell homogenization of the considered sandwich structurerequires neither the analysis of UCs nor of rUCs, this example is useful to validatethe proposed MPBCs for a problem whose exact solution is known analytically.

The facesheets material is composed of multiple layers of Saertex triaxial weavewith Toho Tenax HTS carbon fibres and a polyester knitting yarn, infused withRTM6 epoxy resin [153]. The layup of the individual Non-Crimp Fabric (NCF)weaves is [45◦/0◦/−45◦] and their material properties are listed in Table 5.1 alongwith their thickness hw. The core material is the closed-cell polymetacrylimide foamROHACELL HERO 150 [154], whose elastic properties and density are provided inTable 5.2.

`2

`1

0±

s

¯e3 e1

e2¯e3 e1

e2UC

UC

rUC

0

¯e3e1

e2

s UCC

90±

¯e1 e2

e3

htf

hbfhc

e22

sse1e2

e3

s

¯e3 e1

e2

rUC¯e3e1

e2

UUC

¯e3 e1

e2

e1e2

e3

e1e2

e3

Figure 5.7: Schematic of the sandwich structure with unequal skins considered. Thelatter can be subdivided into UCs, which in turn can be further sub-divided into rUCS. For the solid-to-shell homogenization of this struc-ture we considered as rUC the subdomain denoted as s (in green). TheLCSs of the rUC and of the surrounding subdomains are also shown.

Table 5.1: Elastic properties, density and thickness of the NCF triaxial weave [155,156].

E1 E2 E33 G12 G13 G2354.7 GPa 23.3 GPa 9.0 GPa 24.0 GPa 3.0 GPa 2.9 GPa

ν12 ν13 ν23 hw ρ

0.73 0.26 0.45 0.375 mm 1600 kg/m3

125

Chapter 5

Table 5.2: Elastic properties and density of the HERO G3 150 foam [154].

Ef νf ρf127.5 MPa 0.33 150 kg/m3

The thicknesses of the individual layers of the sandwich structure, along with thelayups for the top and bottom facesheets are provided in Table 5.3. Here, ht

f andhb

f are, respectively, the thickness of the top and bottom facesheet, and hc is thethickness of the foam core (see Figure 5.7); furthermore, the individual NCF weaveorientations are defined relatively to their 0◦-plies.

The solid-to-shell homogenization was performed using the rUC schematicallyshown in Figure 5.7 (denoted as subdomain s), along with its surrounding subdo-mains. The latter are assumed to have LCSs differing from the one of the modelledrUC, such that the corresponding transformation matrices are not equal to the unitymatrix, i.e. Ti �= I.

The facesheets and the core were discretized using eight-noded reduced integrationsolid elements (C3D8R) with enhanced hourglass control, and a constant element sizee = 0.375 mm was adopted; the interfaces between different plies and between thefacesheets and the core were not explicitly modelled, i.e. they share nodes. A size-sensitivity analysis was performed by considering rUCs with different in-plane areaArUC =12, while maintaining a constant aspect ratio

1

2=1. The number of elements

in the laminate plane was varied from 1 to 1600.

5.4.1.2 Results and discussion

For the sandwich structure considered, all components of its equivalent shell stiffnessmatrix K are non-zero. The terms of the latter, computed with the two-scale homo-genization approach described in § 5.2.5, are compared against the reference valuesobtained using CLT.

The relative percentage error for all terms of the K matrix are plotted in Fig-ure 5.8 as function of the number of nodes in the laminate plane indicated withnnodes. All terms of the K matrix, including the bending and twisting terms, as wellas the shear-extension, bending-extension and bending-twisting coupling terms, are

Table 5.3: Sandwich layers’ thicknesses and facesheets layups.

htf hb

f hc Top skin Layup Bottom skin Layup3.0 mm 3.0 mm 10.5 mm [0◦/30◦/60◦, 90◦]s [15◦/45◦/75◦,−15◦]s

126


computed with negligible error, of the order of machine precision. Furthermore, thiserror exhibits, for all terms, a satisfactory convergence trend as the size of the rUCincreases.

Error[10¡6%]

log(nnodes)

0

2

4

6

8

0 1 2 3 4

A_11

A_22

A_66

A11

A22

A66

(a) Aii terms.

-10

0

10

20

0 1 2 3 4

B11

B22

B66

Error[10¡6%]

log(nnodes)

B11

B22

B66

(b) Bii terms.

Error[10¡6%]

-2

0

2

4

0 1 2 3 4

D11

D22

D66

D11

D22

D66

log(nnodes)

(c) Dii terms.

Error[10¡6%]

0

4

8

12

0 1 2 3 4

A_12

A_16

A_26

A12

A26

A16

log(nnodes)

(d) Aij terms (i �=j).

-10

0

10

20

0 1 2 3 4

B12

B16

B26

Error[10¡6%]

B12

B16

B26

log(nnodes)

(e) Bij terms (i �=j).

Error[10¡6%]

0

2

4

6

0 1 2 3 4

D12

D16

D26

D12

D16

D26

log(nnodes)

(f) Dij terms (i �=j).

Figure 5.8: Homogenisation of a composite sandwich structure with unequal skins.The computed terms of the K matrix are compared to the analyticalvalues obtained with CLT, for different rUCs sizes (nnodes is the thenumber of nodes in the laminate plane of the rUC).

127

Chapter 5

5.4.2 Application

5.4.2.1 Description

The aim of this application example is to demonstrate that the current formulationenables to correctly simulate the mechanical response of periodic structures usingrUCs (retrieving the same results as if conventional UCs were used but at a fractionof the CPU time). To this purpose, the MPBCs have been applied for the analysisof the mechanical response of a periodic composite sandwich structure with unequalskins and pin-reinforced core.

The unreinforced sandwich configuration corresponds to that described in§ 5.4.1.1. The reinforcements consist of CFRP pins made of dry intermediate mod-ulus (IM7 12K) [157] carbon fibres infiltrated with RTM6 epoxy resin during themanufacturing of the sandwich panel [158]; the resulting elastic modulus in the fibredirection is Epin =126 GPa.

The pins are inserted through the foam core (without going through or penetratingthe facesheets) in a periodic pattern, see the UC in Figure 5.9. As shown in thelatter, the UC can be reconstructed by tessellation of smaller subdomains (rUCs), byexploiting its internal symmetries; in this example, the analysis is performed usingthe rUC highlighted in green in Figure 5.9. The geometrical parameters of the rUCare provided in Table 5.4. Here, ϕ and β are, respectively, the through-the-thicknessand in-plane stitching angles (see Figure 5.9); 1 and 2 are the in-plane dimensionsof the rUC and φ is the diameter of the pin-reinforcements.

The pin-reinforcements were discretized with two-noded linear beam elements(B31); to simplify the meshing of the pin-reinforced core, the pins’ finite elements wereembedded within the foam core exploiting the Embedded Element Method (EEM)available in Abaqus. This method has been already used for modelling pin-reinforcedcores [159] and textile composites [119]. The optimal ratio between the embedded

beam elements’ size e and the hosting elements’ size h, i.e α =e

hwas determ-

ined so to avoid the characteristic artificial stiffness [160] resulting from the EEMimplementation. In this example, the value of α=1 was adopted.

Table 5.4: Geometrical parameters defining the rUC of the periodic sandwich struc-ture with unequal skins and pin-reinforced core.

1 = 2 φ ϕ β

10.125 mm 1.0 mm 60◦ 45◦

128


`2

`1

s

¯

Á

'

¯e1 e2

e3

0±

¯e3e1

e2

¯e3 e1

e2

¯e3e1

e2

¯e3 e1

e2

s ss

s

¯eee333 e1

e2

s

rUC

UC

B

A

B

A

e2 e1

e3

e2e1

e3

e2

e3

e1

Figure 5.9: Schematic of the UC and rUC considered for the analysis of the mech-anical response of a periodic sandwich structure with unequal skins andpin-reinforced core. The UC’s internal symmetries are exploited for thedefinition of the rUC (in green). The red dotted line A-B indicates thepath along which the membrane strains distribution, computed withthe UC and the rUC models, have been compared.

The sandwich structure was subjected to a twisting loading, i.e. κ12◦ =10−3 rad/mmand the reference solution for this problem was obtained using the entire UC (see Fig-ure 5.9). It is worth-noting that, when using the latter, the MPBCs reduce to thestandard PBCs already used by other authors [161,162].

5.4.2.2 Results and discussion

The in-plane strains (e11, e22, e12) evaluated along the path shown in Figure 5.9 (reddotted line A-B passing through the center of the UC), computed with the rUC andthe UC models are compared, respectively, in Figures 5.10a, 5.10b and 5.10c. Thesestrain values are extracted at each element’s integration point.

129

Chapter 5

-5

-2.5

0

2.5

5

-9 -7 -5 -3 -1 1 3 5 7 9

Series2

Series1

x3 [mm]

e11 ¢104 [¡]UC

rUC

(a) Membrane direct strain e11.

-6

-3

0

3

6

-9 -7 -5 -3 -1 1 3 5 7 9

Series2

Series1

x3 [mm]

e22 ¢104 [¡]UC

rUC

(b) Membrane direct strain e22.

-9

-4.5

0

4.5

9

-9 -7 -5 -3 -1 1 3 5 7 9

Series2

Series1

x3 [mm]

e12 ¢104 [¡]UC

rUC

(c) Membrane shear strain e12.

Figure 5.10: Membrane strains distributions in the UC and rUC models along thepath A-B shown in Figure 5.9 under twisting load. Since the discretiza-tions of the UC and rUC models are the same, the strains distributionscomputed with the UC and the rUC models are, as expected, identical.

130


The total CPU time t associated to the analyses using the rUC and the UCmodels are compared in Figure 5.11a. In addition, results are provided in terms ofboth the modelling/meshing CPU time tM needed for the Python script to createthe FE model and apply the suitable MPBCs (Figure 5.11b), and the analysis CPUtime tA needed to perform the FE analyses (Figure 5.11c). The CPU times providedin Figure 5.11 have been normalized with respect to the values associated with theanalysis performed using the UC model, respectively, tUC , tUC

M and tUC

A .

Since the discretizations of the UC and rUC models are the same, as expected,the membrane strains distributions computed with these models are identical; thisdemonstrate that, through the correct application of the proposed MPBCs, the mech-anical response of periodic structures can be analysed using domains smaller than theUnit Cell, while maintaining the same level of accuracy.

Compared to the UC model, the rUC model has, approximately, 25% of the de-grees of freedom (DOFs) and 50% of the constraint equations (CEs). These attributeslead to a reduction of the total CPU time of about 87% (see Figure 5.11a); the latterresults from the reduction of both the modelling/meshing CPU time (≈ −85%) andthe analysis CPU time (≈−89%), as shown in Figures 5.11b and 5.11c, respectively.

0%

20%

40%

60%

80%

100%

UC rUC

t=tUC

¼13%

(a) Total CPU time.

0%

20%

40%

60%

80%

100%

UC rUC

tM=

tUC

M

¼15%

(b) Modelling/meshing CPUtime.

0%

20%

40%

60%

80%

100%

UC rUC

tA=tUC

A

¼11%

(c) Analysis CPU time.

Figure 5.11: The MPCs enable the use of rUCs for the analysis of the mechanicalresponse of periodic structures. This translates into a reduction of thetotal CPU time required (a), as a result of the reduction of both themodelling/meshing CPU time tM needed to create the FE model andapply the suitable MPBCs (b) and the analysis CPU time tA necessaryto run the FE analysis (c).

131

Chapter 5

Although these specific values are, admittedly, dependent on the model’s sizeand on the specific analysis performed, they confirm the capabilities of the approachproposed in this chapter.

5.5 Conclusions

In this chapter, a novel set of Multiscale Periodic Boundary Conditions (MPBCs)enabling the direct two-scale (solid-to-shell) numerical homogenization of periodicstructures, including their bending and twisting response, using domains smaller thanthe Unit Cells (UCs), named reduced Unit Cells (rUCs), is presented.

Firstly, the validity of proposed MPBCs is demonstrated through the solid-to-shellhomogenization of a composite sandwich structure with unequal skins: all terms ofthe equivalent shell stiffness matrix K are computed with a negligible error, of theorder of machine precision. Secondly, the MPCs are applied to the analysis of themechanical response of a periodic composite sandwich structure with unequal skinsand pin-reinforced core, using a rUC. The results show that the use of the proposedMPBCs allows to correctly simulate the mechanical response of periodic structuresusing rUCs, retrieving the same results as if conventional UCs were adopted.

When compared to MPBCs which do not exploit symmetries, the MPBCs presen-ted in this chapter are shown to achieve a significant reduction in analysis CPU time(approximately 89%), as well as modelling/meshing CPU time (over 85%). Theseresults demonstrate the relevance of the proposed approach for an efficient multiscalemodelling of periodic materials and structures.

132

Chapter 6

Virtual Testing of largecomposite structures: a multiplelength/time-scale framework

6.1 Introduction

6.1.1 Virtual testing of large composite structures

To address complex industrial structural challenges, such as the design of the com-posite central wing box of the Airbus A380 [31], virtual testing methods have beenexploited [32]. The extensive use of virtual testing based on nonlinear FE analysesis envisioned to be a key-aspect towards an increased confidence in the real-scaleand expensive structural tests required for certification; furthermore, virtual testingprovides useful insight into the likelihood, causes and consequences of structural fail-ure [33–35]. However, to be fully established in structural design and certification,virtual testing methods need to be validated against all level of structural testing,from the coupon-level (e.g. material specimens) to the system-level (e.g. wing orfuselage) [36].

Within this framework, the virtual testing of large-scale composite structures forindustrial applications entails significant challenges; the latter are primarily ascrib-able to the inherently multiscale nature of composite materials and, as a result, totheir highly complex failure modes [37]. For the efficient structural design of largecomposite components, their virtual testing often requires that different parts of thestructure are modelled at multiple length- and time-scales, eventually even usingdifferent physics. Therefore, it is crucial to develop:

133

Chapter 6

(i) suitable techniques for coupling areas of the structure modelled at differentlength-scales, i.e. discretized using different finite element types;

(ii) numerical methods to efficiently compute equivalent homogenized properties tobe used in both 2D FE models and in the lower-scale subdomains of multiscaleFE models of large composite components.

6.1.2 Multiscale coupling

The coupling of subdomains discretized with finite elements of different physical di-mension/formulation can introduce artificial stresses at the shared boundaries [89,90].Hence, the stress and strain fields within the structure and its mechanical responsemay not be correctly simulated. This may lead, in problems involving failure, toa low-fidelity damage pattern prediction: furthermore, in dynamic problems, theinterfaces between differently-discretized subdomains might artificially reflect stresswaves [91,92].

Global/local approaches can be categorized as coupled and uncoupled approaches[163]. In the former, the coupling between the local and global models (discretizedusing different finite elements) is enforced through Multi-Point Constraints (MPC).Numerous mixed-dimensional coupling techniques based on MPCs are available inliterature [89, 99, 100]. Compared to standard solid-to-shell coupling, these method-ologies are shown to attenuate the undesirable stress disturbances at the interfacebetween differently-discretized subdomains; nevertheless, the derivation of the suit-able MPC equations for generally orthotropic materials may become computationallyimpracticable [100]. Alternatively, Davila [90] demonstrated that the use of ad-hoctransition elements allows to faithfully resolve the stress and strain fields within boththe local and global models in the proximity of the interfaces.

In uncoupled global/local approaches, the displacements/tractions obtained from ahigher-scale global model are used to prescribe the boundary conditions for the lower-scale local model. The series of higher-scale analysis followed by a lower scale analysiscan be run once (a typical example is the Submodelling technique available in Abaqus[85]), or iteratively until convergence of the forces and moments on the local/globalinterface is obtained. As a result, these approaches are difficult to automate and posesignificant challenges in terms of computational resources; furthermore, they might beinaccurate if damage and failure propagation within the local model significantly affectthe response of the global model. Reinoso et al. [104,105] compared the Submodellingtechnique (uncoupled approach) to Shell-to-Solid coupling (coupled approach). In thelatter works, stress disturbances were observed at the local/global boundaries for both

134

Virtual Testing of large composite structures: a multiple length/time-scaleframework

approaches; as a result, the correct identification of the damage extension might bejeopardized.

To enhance the computational efficiency of multi-scale approaches, it is desirablethat the size of the lower-scale model is kept to a minimum. However, the global/localtransition should be sufficiently distant from any perturbation such as boundaries ordamaged zones which could potentially interact with the stress disturbance at the dis-cretization transition. The influence of the distance of the global/local transition fromthe delamination front was investigated by Krueger et al. for the cases of delaminatedtest specimens [106, 107] as well as skin/stringer debonding in an composite aircraftcomponent [108,109].

The need of having the global/local transition at a sufficient distance from anyboundaries/damaged areas often leads to overly large models at the lower-scales; asa result, the computational efficiency of the analysis might not be optimal. Thus, tofully exploit the computational advantages provided by multiscale approaches, it isof paramount importance to use a local/global coupling technique which:

(i) avoids the presence of artificial stress disturbances, as well as the unrealisticstress wave reflections observed at the boundaries between differently discretizedsubdomains;

(ii) allows to minimize the size of the areas that are required to be modelled at thelower length-scales.

6.1.3 Solid-to-shell homogenization

For the virtual testing of large composite structures, FE (preferably 2D) modelswith homogenised material properties are usually preferred to high-fidelity 3D FEmodels. Furthermore, shell-based FE models with equivalent homogenized propertiesare used in the higher-scale subdomains of multiscale FE models of large compositecomponents.

The numerical homogenization of periodic structures is often carried out throughthe analysis of the mechanical response of Unit Cell (UC) models, where the com-posite microstructure is modelled at the lowest length/scale of interest. Numerousstudies focused on the determination and use of UCs [137, 138], as well as on theappropriate periodic boundary conditions (PBCs) that need to be applied for theiranalysis [139–142]. However, as a result of the increasing complexity of composite mi-crostructures, e.g. textile and NCF composites, pin-reinforced sandwich structures,etc., the computational cost associated to modelling/meshing and analysis might be-come burdensome.

135

Chapter 6

To circumvent this difficulty, several authors have proposed to exploit (wheneverthey exist) internal symmetries of the UCs, thus enabling a significant reduction ofthe analysis domain. Analysis domains smaller than UCs are denoted as reduced UnitCells (rUCs) and have successfully been used for the analysis UD composites [143],particle-reinforced composites [144] and textile composites [145–147]. However, inthese works address the issue of using rUCs for obtaining equivalent 3D homogen-ised properties, without considering the computation of equivalent 2D homogenisedproperties.

The efficient computation of the equivalent homogenized properties to be usedin both 2D FE models and in the higher-scale subdomains of multiscale FE mod-els of large composite components hinges upon the development of a mathematicalframework which:

(i) enables the direct two-scale (solid-to-shell) homogenization of periodic struc-tures, including their bending and twisting response;

(ii) enables the use of rUCs.

6.1.4 Structure of this chapter

In this chapter, we illustrate the capabilities of a multiple length/time-scale frameworkfor the virtual testing of large composite structures. The multiple length/time-scaleframework is described in § 6.2 and its primitives are illustrated and discussed, re-spectively in § 6.3 and § 6.4, through the simulation of a real-sized helicopter rotorblade subjected to a low-velocity impact; conclusions are drawn in § 6.5.

6.2 Multiple length/time scale framework

The multiple length/time-scale framework used in this chapter, and graphically illus-trated in Figure 6.1, consists of (i) the Mesh Superposition Technique (MST) [136](see Chapter 4) for the progressive transition between subdomains discretized withfinite elements of different physical dimension/formulation and (ii) a mathemat-ical framework which exploits UCs internal symmetries in the context of the directsolid-to-shell homogenization of periodic structures [164] (see Chapter 5).

The key-aspects of the MST are:

(i) unlike conventional solid-to-shell coupling techniques based on a suddendiscretization-transition, the MST eliminates the artificial stress disturbancesat the shared boundaries between differently discretized subdomains;

136


Multiple length/time scale framework for the virtual testing of large composite structures

Mesh Superposition Technique (MST)

Progressive element-type transition between differently discretized subdomains

bEA

bEB

EsA

=

MST regionEsB

S SS

Accurate damage pattern prediction Model size & CPU time reduction

0%

50%

100%

Normalized CPU time

Sudden Tr.

Sudden Tr.

MST

MST

-60% Model size

-23% CPU time

Solid-to-shell homogenization exploiting symmetries

Periodic composite structure

Unit Cells (UCs) vs Reduced Unit Cells (rUCs)

Modelling/Analysis CPU time reduction

A

B

B

A

0%

50%

100%

Multi length (-23%)

Multi length/time

(-49%)

Model size CPU time

UC

rUC 0%

50%

100%

0%

50%

100%

rUC UC

UC rUC

-85% CPU Modelling time

-89% CPU Analysis time

Figure 6.1: Multiple length/time-scales framework for the virtual testing of largecomposite components. Such framework consists of a Mesh Superposi-tion Technique for coupling differently-discretized subdomains (A) andon the exploitation of symmetries in the solid-to-shell numerical homo-genization of periodic structures.

137

Chapter 6

(ii) using the MST, the size of the lower-scale models (the most computation-ally demanding) can be minimized without jeopardizing the response at theglobal/local transition. Therefore, compared to conventional solid-to-shell coup-ling techniques based on a sudden discretization-transition the MST enables ahigh-fidelity damage pattern prediction at a lower computational cost;

(iii) the progressive transition provided by the MST mitigates the spurious reflec-tions of stress waves at the interfaces between differently-discretized subdomains(see example in Appendix A);

(iv) the MST can be used in combination with an implicit/explicit co-simulationtechnique [127, 128], for a multiple time/length-scale analysis. The capabilityprovided by the MST to minimize the size of the lower-scale models allowsto maximize the computational efficiency of the implicit/explicit co-simulationtechnique; hence, in the context of a multiple length/time scale analysis, the useof the MST for the length-scale transitions (instead of a sudden discretization-transition) results into a significant computational cost reduction.

Regarding the exploitation of symmetries in the solid-to-shell homogenization ofperiodic structure, the following points should be highlighted:

(i) the framework proposed in [164] leads to the derivation of the exact periodicboundary conditions that apply to rUCs and enable the numerical solid-to-shellhomogenization of periodic structures, including their bending and twistingresponse.

(ii) no limitations on the deformed/undeformed shape of the rUCs, as well as tothe nonlinear behaviour at the lowest length-scale were made (provided thereis no localization [151]);

(iii) when comparing results obtained using conventional UCs with those obtainedwith rUCs, in the latter case time savings of about 90% can be achieved in theanalysis CPU time, as well as in the modelling/meshing CPU time.

6.3 Multiple length/time-scale simulation of a largeaeronautical component

6.3.1 Problem description

The multiple length/time-scale framework detailed in § 6.2 is applied to the analysisof a low-velocity impact on real-sized helicopter rotor blade. The latter is idealized asan hollow structure with dimensions provided in Figure 6.2. The profile of the rotor

138


blade corresponds to the NASA/Langley Whitcomb integral supercritical airfoil [165].The rotor blade is impacted with a 2.5 kg impactor of diameter equal to 40 mm anda 12 J impact energy was considered.

The profile of the rotor blade is assumed to be made of a composite sandwichstructure [158] with pin-reinforced foam core [166]; relevant geometrical and materialproperties of the pin-reinforced sandwich structure can be found in [164]. The materialconsidered for the pin-reinforcements is the carbon-epoxy T300/913 [69].

3500

200

Impactor 115

1750 Á 40

9 0 ±

0 ±

Figure 6.2: Schematic (not in scale) of the helicopter rotor blade considered in thisstudy; dimensions are in mm.

6.3.2 Implicit & Explicit FE submodels

The multiple time-scale connotation of the framework described in § 6.2 consistsin simulating the mechanical response of different portions of the structure usingdifferent solvers, depending on where they are expected to provide the most compu-tationally efficient solution. Generally, complex material failure, contact interactionsand highly nonlinear response are best analysed using FE solvers based on explicittime-integration schemes, e.g. Abaqus/Explicit, while the elastic behaviour of lightand stiff components can, more efficiently, be simulated with FE solvers using implicittime-integration schemes, e.g. Abaqus/Standard.

Using different FE solvers implies the definition of multiple time-scales at whichthe structural response is analysed: (i) an explicit/micro time-scale, characterizedby a high number of short and relatively inexpensive time-steps and (ii) an impli-cit/macro time-scale, characterized by a reduced number of larger time-steps, due tothe unconditional stability of implicit solvers.

Therefore, for the problem described in § 6.3.1, the FE model of the rotor blade(see Figure 6.3) consists of two separate FE submodels:

139

Chapter 6

(i) an explicit submodel for the area of the rotor blade surrounding the impactlocation (where contact interactions and complex failure modes are expected);within this submodel, different parts of the structure are modelled at differentlength-scales as described in Chapter 4;

(ii) an implicit submodel for the remaining of the structure, whose deformationremains elastic throughout the entire analysis.

Explicit multiscale FE submodel Implicit FE submodelExplicit multiscale FE submodel Implicit FE submodel

Multi time/length-scale FE model

Figure 6.3: Implicit/explicit submodels of a multiple length/time-scale FE modelof an helicopter rotor blade. The area surrounding the impact locationis analysed using an explicit solver, to exploit its capabilities to betterhandle contact interactions and complex failures modes; the responseof the remaining of the structure (the largest part) is more efficientlysimulated using an implicit solver.

In the implicit FE subdomain, the rotor blade is discretized using first-ordertriangular shell elements (S3). The equivalent homogenised 2D behaviour, i.e. theentire ABD matrix of the pin-reinforced sandwich structure is defined exploiting the

140


solid-to-shell homogenization described in Chapter 5. The meshing of the entire rotorblade required 63744 S3 elements.

In this example, the multiple length/time-scale analysis is performed exploitingAbaqus built-in Co-Simulation Engine (CSE) [85]; the simulation of the mechan-ical response of the implicit and explicit submodels is performed using, respectively,Abaqus/Explicit and Abaqus/Standard.

6.3.3 Multiscale explicit FE submodel

Within the explicit submodel, different parts of the structure are modelled at differentlength-scales and the coupling between them is performed using the MST [136]. Thesize of the differently-discretized subdomains needs to be established a priori, i.e.before the analysis, for example by exploiting analytical methods to estimate theexpected size of the damaged areas [167].

For the example described in this chapter, we considered three different regionswhere the structure is modelled at different length-scales (see Figure 6.4):

Meso-scale subdomainWithin this subdomain, the structure is modelled at the lowest length-scale.The plies in the composite facesheets, the interfaces between them and betweenthe facesheets and the core, as well as the pin-reinforcements within the core,are modelled individually.

Eight-noded reduced-integration solid elements (C3D8R) with enhanced hour-glass control were used for the discretization of the composite plies, pin-reinforcements and core layer: eight-noded cohesive elements (COH3D8) wereused to account for possible delamination (between composite plies) and de-bonding (between the facesheets and the core). To simplify the meshing of thepin-reinforced core, the pins’ finite elements were embedded within the foamcore exploiting the Embedded Element Method (EEM) available in Abaqus.This method has been already used for modelling pin-reinforced cores [159] andtextile composites [119]. The discretization of the meso-scale subdomain (in-cluding the MST transition regions) required 71896 C3D8R elements and 11520COH3D8 elements.

The complex shape of the pin-reinforcements (i.e. misalignments, bended ends,etc.) is also explicitly modelled. Within the meso-scale subdomain, severalfailure modes are accounted for (see Figure 6.4), i.e

(i) delamination between composite plies, as well as debonding between thecomposite facesheets and the core, using multiple layers of cohesive ele-

141

Chapter 6

Explicit multiscale FE submodel

Meso-scale subdomain Macro-scale subdomain

Debonding

Core Crushing Pins’ failureDelamination/Debonding

Meso/macro-scale subdomain

Figure 6.4: Multiscale explicit FE submodel. Different areas of the structure can bemodelled at different length-scales and their coupling be performed us-ing the MST. The mechanical response of the structure can be correctlycaptured at all the length-scales of interest (both in terms of geomet-rical details and failure modes), while keeping the computational costof the analysis to a minimum.

142


ments with a mixed-mode bilinear traction-separation law to simulate thesoftening and fracture response [125,126];

(ii) crushing of the foam core [166,168];

(iii) failure of the pin reinforcements [69];

Meso/Macro scale regionWithin this subdomain, the composites facesheets, the core, the interfacebetween them, as well as the pin-reinforcements within the core layer, are ex-plicitly modelled.

Eight-noded reduced-integration continuum shell elements (SC8R) with en-hanced hourglass control were used for the discretization of the compositefacesheets while eight-noded reduced-integration solid elements (C3D8R) wereused for the discretization of the pin-reinforcements and of the core. As in themeso-scale subdomain, the pins’ finite elements were embedded within the coreexploiting the EEM.

The discretization of the meso-scale subdomain (including the MST transitionregions) required 181952 C3D8R elements, 54784 SC8R elements and 11520COH3D8 elements.

Within the meso/macro-scale subdomain, the only failure model accounted foris the debonding between the core and the composite facesheets (see Figure 6.4).This can actually occur even at a significant distance from the impact location,as it might initiates as a result of the stress concentrations due to the pinreinforcements in the foam core layer.

Macro-scale regionWithin this subdomain, the pin-reinforced sandwich structure is discretized us-ing reduced-integration first-order quadrangular shell elements (S4R) with en-hanced hourglass control. The equivalent homogenised 2D behaviour, i.e. thecorresponding entire ABD matrix of the sandwich structure is computed ex-ploiting the solid-to-shell homogenization described in § 6.2. The discretizationof the meso-scale subdomain (including the MST transition regions) required1516 S4R elements.

6.4 Discussion

The key-aspects of the multiple length/time-scale framework for the virtual testing oflarge composite structures illustrated in this chapter can be summarized as follows:

143

Chapter 6

(i) coupling differently-discretized subdomains using the MST [136] presentedin Chapter 4 allows to minimize the areas of the structure modelled at thelowest (and therefore computationally costly) length-scales. Thus, the desiredlevel of accuracy is attained only in selected portions of the model, which allowsnot to jeopardize the computational cost of the analysis;

(ii) the response of different portions of the structure is simulated using differentsolvers, depending on where they are expected to provide the most efficientsolution;

(iii) the exploitation of symmetries in the solid-to shell homogenization [164] presen-ted in Chapter 5 allows to strongly reduce the time required for the computationof equivalent homogenised properties;

(iv) the size of the resulting multiscale FE models can be minimized without anyloss of accuracy at the lowest-scales; as a result, such FE models can be easilyhandled by high-performance computing resources, as well as modern laptopcomputers. As an example, the multiscale FE model of the helicopter rotorblade discussed in § 6.3 consisted of a total number 1291202 DOFs (410846nodes); simulations were performed using a laptop computer with an Intel i7-3610QM processor (2.3GHz Max Turbo Frequency), 8 CPUs and 8Gb RAM fora total analysis time of approximately 6 hours.

6.5 Conclusions

In this chapter, a multiple length/time-scale framework for the virtual testing of largecomposite structures is presented. Its primitives and key-aspects are detailed throughthe analysis of a real-sized helicopter rotor blade subjected to a low-velocity impact.The proposed multiple length/time-scale framework consists of the Mesh Superposi-tion Technique (MST) presented in Chapter 4 and the solid-to-shell homogenizationpresented in Chapter 5.

The multiple length/time-scale framework enables a significant reduction of theCPU time required to compute the homogenized properties used in the higher-scalesof multiscale FE models; furthermore, it allows to minimize the areas of the struc-ture that need to be modelled at the lowest length-scales, by opportunely couplingdifferently-discretized subdomains using the MST.

The multiple length/time-scale framework presented in this work represents aclear step towards the systematic and efficient virtual testing of large composite com-ponents.

144

Chapter 7

Conclusions

7.1 Introduction

The individual contributions presented in this work, i.e. the analytical model forpredicting the post-crushing response of polymeric foam cores of chapter 2, the meas-urement of the translaminar fracture toughness of Non-Crimp Fabric composites ofchapter 3, the Mesh Superposition Technique (MST) of chapter 4 and the MultiscalePeriodic Boundary Conditions of chapter 5 are extensively discussed in § 2.5, § 3.4,§ 4.5, and § 5.4, respectively. In the following, the above chapters’ contributions tothe state of the art and impact are summarized.

7.2 Novelty

The following contributions of this work to the state of the art are highlighted:

(i) the analytical model presented in chapter 2 is the first (to date) capable ofaccurately predicting the post-crushing compressive response of crushable foamsusing, as input, experimental measurements obtained exclusively from standardmonotonic compressive tests;

(ii) the experimental analysis in chapter 3 represents the first (to date) attempt inthe literature to measure the translaminar fracture toughness of NCF compos-ites with multiaxial blankets; furthermore, the concept of an homogenised NCFblanket-level translaminar fracture toughness is introduced for the first time;

(iii) chapter 3 demonstrates that the translaminar fracture toughness of off-axis fibretows/NCF blankets can be analytically related to that of axially-loaded fibretows/NCF blankets;

145

Chapter 7

(iv) in chapter 4, a novel Mesh Superposition Technique (MST) for coupling dif-ferent areas of the composite structure modelled at different length-scales andwhose associated discretizations consist of different element types is presented.The interface between local and global meshes is replaced by transition regionswhere the corresponding discretizations are superposed. Within these transitionregions, a model based on partition of unity is derived to achieve a progressivediscretization-transition;

(v) chapter 4 proposes the first (to date) attempt towards a multiple length/time-scale modelling of composite structures experiencing localized damage; and

(vi) chapter 5 presents the first set of Periodic Boundary Conditions for the solid-to-shell numerical homogenization of periodic structures, including their bendingand twisting response, using domains smaller than the Unit Cell, obtained byexploiting the internal symmetries of the latter.

7.3 Impact

The potential impact of the findings and conclusions presented in this work can besummarized in the following points:

(i) the calibration of the model presented in chapter 2 is performed using experi-mental measurements obtained exclusively from standard monotonic compress-ive tests; therefore, the need for performing time-consuming compressive testsincluding multiple unloading-reloading cycles is avoided and the effective testingtime significantly reduced;

(ii) the measured values of the translaminar fracture toughness of NCF composites(see chapter 3) can be used as input value of physically-based failure criteria topredict/analyse failure of NCF composites;

(iii) results presented in chapter 3 indicate that the translaminar fracture toughnessof laminates with complex layups (with several differently-oriented off-axis fibretows and off-axis NCF blankets) can be accurately estimated from the translam-inar fracture toughness of axially-loaded fibre tows/NCF blankets. This resultis highly relevant for the design of NCF composite laminates to be used inlarge-scale structural applications;

(iv) coupling differently-discretized subdomains using the MST presented in Chapter4 allows to minimize the areas of the structure modelled at the lowest (and there-fore computationally costly) length-scales. Thus, the desired level of accuracyis attained by refining only selected portions of the model, which allows not tojeopardize the computational cost of the analysis;

146

Conclusions

(v) when compared to PBCs which do not exploit symmetries, the PBCs presentedin chapter 5 are shown to achieve a significant reduction in analysis CPU time,as well as modelling/meshing CPU time. Therefore, homogenized properties ofvery complex periodic materials and structures can be computed with a muchlower use of computational resources and time; and

(vi) the multiple length/time-scale framework presented in chapter 6 represents aclear step towards the systematic and efficient virtual testing of large compos-ite components. The size of the resulting multiple length/time-scale FE modelscan be minimized without any loss of accuracy at the lowest-scales (chapter4), the response of different portions of the structure is simulated using dif-ferent solvers - depending on where they are expected to provide the mostefficient solution (chapter 4) - and the exploitation of symmetries allows tostrongly reduce the time required for the computation of equivalent homogen-ised properties (chapter 5). As a result, such FE models can be easily handledby high-performance computing resources, as well as modern laptop computers.

7.4 Future work

The work presented in this thesis offers several opportunities for further development.These include, but are not limited to, the points highlighted in the following:

(i) the analytical model presented in chapter 2 predicts the post-crushing com-pressive response of polymeric foam materials based on the computation of thethickness of crushed and uncrushed layers of material. It would be relevant forpractical purposes, to extended it to predict the post-crushing tensile stiffnessand strength of polymeric foam materials;

(ii) the analytical model presented in chapter 2 can serve as basis to develop apredictive model capable of predicting the effect of prior crushing on the en-tire constitutive response of polymeric foam materials, including multiaxial andshear loadings;

(iii) at this stage of development, chapter 2 requires as input experimental datafrom standard monotonic compressive tests. For more complicated scenarios,e.g. pin-reinforced foam cores, stress vs. strain curves obtained through FEanalyses could be used as input data;

(iv) in chapter 3, CT specimens are used to measure the translaminar fracture tough-ness of a carbon NCF composite. Different specimen’s geometry could be usedto obtained meaningful propagation translaminar fracture toughness values and,therefore, a complete R-curve;

147

Chapter 7

(v) although it would lead to greater challenges than those associated to the meas-urement of the tensile translaminar fracture toughness, the measurement ofthe compressive translaminar fracture toughness is of paramount importancefor damage-tolerant design of composite structures. Following an experimentalprocedure similar to that described in chapter 3, the compressive translaminarfracture toughness of NCF composites could be measured;

(vi) in chapter 3, the translaminar fracture toughness of off-axis fibre tows/NCFblankets is analytically related to that of axially-loaded fibre tows/NCFblankets. However, further verification would be required for other values ofthe orientation angle α;

(vii) the MST presented in chapter 4 requires at this stage of development that allnodes belonging to the higher-scale discretization coincide with a subset of thenodes belonging to the lower-scale discretization. To render the methodologymore flexible, it would be of interest to eliminate this constraint and let thetwo discretizations be mutually independent. This would simplify the practicalimplementation of the proposed MST and improve its usability;

(viii) the MST presented in chapter 4 could be coupled with adaptive meshing meth-ods. In this manner, the size of the lower-scale subdomains can be minimizedat any moment in time during the FE analysis, their shape would not need tobe defined prior to the analysis and, more importantly, their shape and size canadapt to evolving features in the model, e.g. propagating cracks;

(ix) to implement the Multiscale Periodic Boundary Conditions presented in chapter5 in a FE code, it is required that the discretizations on the boundaries of therUC are periodic; in many practical applications, e.g. textile composites, thisentails an enormous amount of user modelling time or it can even be practicallyimpossible to achieve. Therefore, it would be beneficial to develop a robustnumerical methodology to eliminate this constraint; and

(x) the proposed formulation of the Multiscale Periodic Boundary Conditions inchapter 5 could be extended so that also the computation of the equivalenttransverse shear stiffness is enabled.

148

Bibliography

[1] J. Tomblin, T. Lacy, B. Smith, S. Hooper, A. Vizzini and S. Lee, “Review ofdamage tolerance for composite sandwich airframe structures,” DOT/FAA/AR-99/49 Federal Aviation Authority, USA, 1999.

[2] J. Vinson, “Sandwich structures,” Applied Mechanics Review, vol. 54, no. 3,pp. 201–214, 2001.

[3] S.-E. Hellbratt, “Experiences from design and production of the 72 m CFRP-sandwich corvette Visby.,” in Proceedings of the 6th International Conferenceon Sandwich Structures (J.R. Vinson, Y.D.S. Rajapakse and L.A. Carlsson,ed.), pp. 15–24, CRC Press, 2003.

[4] G. Belingardi, M.P. Cavatorta and R. Duella, “Material characterization of acomposite-foam sandwich for the front structure of a high-speed train,” Com-posite Structures, vol. 61, no. 1, pp. 13–25, 2003.

[5] J-S. Kim, S-J. Lee and K-B. Shin, “Manufacturing and structural safety eval-uation of a composite train carbody,” Composite Structures, vol. 78, no. 4,pp. 468–476, 2007.

[6] P.M. Schubel, J-J. Luo and I.M. Daniel, “Impact and post impact behavior ofcomposite sandwich panels,” Composites Part A: Applied Science and Manu-facturing, vol. 38, no. 3, pp. 1051–1057, 2007.

[7] D. Zenkert, A. Shipsha, P. Bull and B. Hayman, “Damage tolerance assessmentof composite sandwich panels with localised damage,” Composites Science andTechnology, vol. 65, no. 15, pp. 2597–2611, 2005.

[8] P.M. Schubel, J-J. Luo and I.M. Daniel, “Low velocity impact behavior ofcomposite sandwich panels,” Composites Part A: Applied Science and Manu-facturing, vol. 36, no. 10, pp. 1389–1396, 2005.

149

Bibliography

[9] V.L. Tagarielli, V.S. Deshpande and N.A. Fleck, “The dynamic response of com-posite sandwich beams to transverse impact,” International Journal of Solidsand Structures, vol. 44, no. 7, pp. 2442–2457, 2007.

[10] R. Olsson, “Engineering method for prediction of impact response and damagein sandwich panels,” Journal of Sandwich Structures and Materials, vol. 4, no. 1,pp. 3–29, 2002.

[11] R. Olsson, “Methodology for predicting the residual strength of impacted sand-wich panels,” tech. rep., FFA TN 1999-08, The Aeronautical Research Instituteof Sweden, 1999.

[12] R. Avila-Dominguez, “NCF processability and performance - A literature re-view,” DERA/MSS/MSMA2/WP003564, 2000. Farnborough, UK.

[13] C. Kjelgaard, “Challenges in composites,” Aircraft Technology, vol. 116, pp. 52–57, 2012.

[14] A. Quilter, “Composites in aerospace applications,” IHS - The Source for Crit-ical Information and Insight, 2010.

[15] G. Marsh, “Composites lift off in primary aerostructures,” Reinforced Plastics,vol. 48, no. 4, pp. 22–27, 2004.

[16] E. Greenhalgh, Failure analysis and fractography of polymer composites. El-sevier, 2009.

[17] S. Lomov, Non-crimp fabric composites: manufacturing, properties and applic-ations. Elsevier, 2011.

[18] G.A. Bibo, P.J. Hogg, R. Backhouse and A. Mills, “Carbon-fibre non-crimp fab-ric laminates for cost-effective damage-tolerant structures,” Composites Scienceand Technology, vol. 58, no. 1, pp. 129–143, 1998.

[19] V. Cecen, M. Sarikanat, H. Yildiz and I.H. Tavman, “Comparison of mechanicalproperties of epoxy composites reinforced with stitched glass and carbon fabrics:Characterization of mechanical anisotropy in composites and investigation onthe interaction between fiber and epoxy matrix,” Polymer Composites, vol. 29,no. 8, pp. 840–853, 2008.

[20] F. Edgren, L.E. Asp and P.H. Bull, “Compressive failure of impacted ncf com-posite sandwich panels-characterisation of the failure process,” Journal of Com-posite Materials, vol. 38, no. 6, pp. 495–514, 2004.

[21] F. Edgren, Physically based engineering models for NCF composites. PhD thesis,KTH, 2006.

150

Bibliography

[22] J. Filsinger, R.-P. Dittmann and T. Bischoff, “Textile preforming technolo-gies for fuselage applications,” in Proceedings of the 25th Jubilee InternationalSAMPE European Conference 2004 of the Society for the Advancment of Ma-terials and Process Engineering, Materials & Process Technnology - The driverfor tomorrow’s improved performance, 2004.

[23] “BMW formally launches i3 manufacture and assembly.” http://www.compositesworld.com/news/bmw-formally-launches-i3-manufacture-and-assembly, Sept. 2013.

[24] P. Martensson, D. Zenkert and M. Akermo, “Method for the cost-efficient andweight-efficient material diversity and partitioning of a carbon fibre compositebody structure,” Proceedings of the Institute of Mechanical Engineering, PartD: Journal of Automotive Engineering, 2015.

[25] L. Asp, “Local models for NCF composite materials mechanical preformanceprediction,” in Proceedings of the 16th International Conference on CompositeMaterials, 2007.

[26] FALCOM, “Failure, performance and processing prediction for enhanced designwith non-crimp-fabric composites,” EU-project, contract No. G4RD-CT-2001-00694, 2001.

[27] TECABS, “Technologies for carbon fibre modular automotive body structures,”Framework V project G6RD-CT2001-00636, 2000.

[28] I. Verpoest, T.C. Thanh and S. Lomov, “The TECABS Project: Developmentof Manufacturing, Simulation and Design Technologies for a Carbon Fiber Com-posite Car,” in Proceedings of the 9th Japan International SAMPE Symposium,pp. 56–61, 2005.

[29] H. Molker, D. Wilhelmsson, R. Gutkin and L.E. Asp, “Orthotropic criteriafor transverse failure of non-crimp fabric reinforced composites,” Journal ofComposite Materials, 2015.

[30] H. Molker, R. Gutkin and L.E. Asp, “Implementation of failure criteria fororthotropic Non Crimp Fabric composite structures,” in Proceedings of the 5thECCOMAS Thematic Conference on the Mechanical Response of Composites -COMPOSITES 2015, 2015.

[31] J. Pora, “Composite Materials in the Airbus A380 - From History to Future,”in Proceedings of the 13th International Conference on Composite Materials,2001.

151

Bibliography

[32] M.G. Ostergaard, A.R. Ibbotson, O. Le Roux and A.M. Prior, “Virtual testingof aircraft structures,” CEAS Aeronautical Journal, vol. 1, no. 1-4, pp. 83–103,2011.

[33] J. Imbert, “Airbus, challenges in aircraft structure analysis,” in ESA/NAFEMSSeminar on Engineering Quality, verification and validation. Noordwijk, 6December 2007.

[34] A. Prior, “Dassault Systemes, simulating damage and failure in aircraft struc-tures,” in RAeS Conference: challenges for the next generation-concept to dis-posal, October 14âĂŞ16, 2008.

[35] A. Prior, “Dassault Systemes, nonlinear simulation of large scale aircraft struc-turesâĂŤimplications for certification methodology and high performance com-puting infrastructure,” in NAFEMS World Congress, 2009.

[36] T. Brown, “Airbus, Working to meet the challenges of ext generation com-posite wing structural design,” in RAeS Conference: challenges for the nextgeneration-concept to disposal, October 14âĂŞ16, 2008.

[37] European Aviation Safety Agency, “Certification Specifications for Large Aero-planes CS-25 Amendement 5,” 2008.

[38] V. Koissin, A. Shipsa and V. Skvortsov, “Compression strength of sandwichpanels with sub-interface damage in the foam core,” Composites Science andTechnology, vol. 69, pp. 2231–2240, 2009.

[39] S. Abrate, “Localized impact on sandwich structures with laminated facings,”Applied Mechanics Review, vol. 50, no. 2, pp. 69–82, 1997.

[40] S. Minakuchi, T. Uezono and N. Takeda, “Formation and relaxation of residualfacesheet dent on foam-core sandwich structures by localized transverse load-ing,” Journal of Sandwich Structures and Materials, vol. 15, no. 21, pp. 7–91,2013.

[41] V. Koissin and A. Shipsa, “Residual in-plane mechanical properties of trans-versely crushed structural foams,” Journal of Sandwich Structures and Mater-ials, vol. 11, pp. 199–211, 2009.

[42] Q. Liu and G. Subhash, “A phenomenological constitutive model for foamsunder large deformations,” Polymer Engineering and Science, vol. 44, no. 3,pp. 463–473, 2004.

[43] Q. Liu, G. Subhash and X.L. Gao, “A parametric study on crushability of open-cell structural polymeric foams,” Journal of Porous Materials, vol. 12, no. 3,pp. 233–248, 2005.

152

Bibliography

[44] M. Avalle, G. Belingardi and A. Ibba, “A parametric study on crushability ofopen-cell structural polymeric foams,” International Journal of Impact Engin-eering, vol. 34, no. 1, pp. 3–27, 2007.

[45] E.A. Flores-Johnson, Q.M. Li and R.A.W. Mines, “Degradation of elastic mod-ulus of progressively-crushable foams in uniaxial compression,” Journal of Cel-lular Plastics, vol. 44, no. 5, pp. 415–434, 2008.

[46] D. Zenkert, A. Shipsha and K. Persson, “Static indentation and unloading re-sponse of sandwich beams,” Composites Part B: Engineering, vol. 35, pp. 511–522, 2004.

[47] V. Rizov, A. Shipsha and D. Zenkert, “Indentation study of foam core sandwichcomposite panels,” Composite Structures, vol. 69, pp. 95–102, 2005.

[48] J. Zhang, Z. Lin, A. Wong, N. Kikuch, V.C. Li, A.F. Yee and G.S. Nusholtz,“Constitutive modeling and material characterization of polymeric foams,”ASME Journal of Engineering Materials and Technology, vol. 119, no. 3,pp. 284–291, 1997.

[49] J. Zhang, N. Kikuchi, V.C. Li, A.F. Yee and G.S. Nusholtz, “Constitutivemodeling of polymeric foam material subjected to dynamic crash loading,” In-ternational Journal of Impact Engineering, vol. 21, no. 5, pp. 369–386, 1998.

[50] L. Gibson and M. Ashby, Cellular solids: structure and properties. Cambridgeuniversity press, 1997.

[51] Q.M. Li, R.A.W. Mines and R.S. Birch, “The crush behaviour of rohacell-51WF structural foam,” International Journal of Solids and Structures, vol. 37,pp. 6321–6341, 2000.

[52] S. Arezoo, V.L. Tagarielli, N. Petrinic and J. M. Reed, “The mechanical re-sponse of rohacell foams at different length scales,” Journal of materials science,vol. 46, pp. 6863–6870, 2011.

[53] C. Kelley, Iterative methods for optimization, vol. 18. Siam, 1999.

[54] M.A. Rodrıguez-Perez, S. Dıez-Gutierrez and J.A. De Saja, “The recovery be-havior of crosslinked closed cell polyolefin foams,” Polymer Engineering andScience, vol. 38, no. 5, pp. 831–837, 1998.

[55] M.A. Rodrıguez-Perez, J.I. Velasco, D. Arencon, O. Almanza and J.A. De Saja,“Mechanical characterization of closed-cell polyolefin foams,” Journal of AppliedPolymer Science, vol. 75, no. 1, pp. 156–166, 2000.

153

Bibliography

[56] V. Koissin and A. Shipsa, “Deformation of foam cores in uniaxial compression-tension cycle,” Journal of Sandwich Structures and Materials, vol. 8, no. 5,pp. 395–406, 2006.

[57] M.F. Caliri Junior, G.P. Soares, R.A. Angelico, R.B. Canto and V. Tita, “Studyof an anisotropic polymeric cellular material under compression loading,” Ma-terials Research, vol. 15, no. 3, pp. 359–364, 2012.

[58] J.D. Ferry, Viscoelastic properties of polymers. John Wiley & Sons, 1980.

[59] ASTM International, ASTM D5229, Standard test methods for moisture ab-sorption properties and equilibrium conditioning of polymer matrix compositematerials.

[60] ASTM International, ASTM C365-57, Standard test methods for flatwise com-pressive strength of sandwich cores.

[61] E. Marklund, L.E. Asp and R. Olsson, “Transverse strength of unidirectionalnon-crimp fabric composites: Multiscale modelling,” Composites Part B: En-gineering, vol. 65, pp. 47–56, 2014.

[62] S. Drapier and M.R. Wisnom, “Finite-element investigation of the compressivestrength of non-crimp-fabric-based composites,” Composites Science and Tech-nology, vol. 59, no. 8, pp. 1287–1297, 1999.

[63] L.G. Zhao, N.A. Warrior and A.C. Long, “Finite element modelling of damageprogression in non-crimp fabric reinforced composites,” Composites Science andTechnology, vol. 66, no. 1, pp. 36–50, 2006.

[64] A. Gonzalez, Enrique Graciani and Federico Parıs, “Prediction of in-plane stiff-ness properties of non-crimp fabric laminates by means of 3D finite elementanalysis,” Composites Science and Techonology, vol. 68, no. 1, pp. 121–131,2008.

[65] J. Juhasz, R. Rolfes and K. Rohwer, “A new strength model for application ofa physically based failure criterion to orthogonal 3D fiber reinforced plastics,”Composites Science and Technology, vol. 61, no. 13, pp. 1821–1832, 2001.

[66] P.P. Camanho, C.G. Davila, S.T. Pinho, L. Iannucci and P. Robinson, “Pre-diction of in situ strengths and matrix cracking in composites under transversetension and in-plane shear,” Composites Part A: Applied Science and Manufac-turing, vol. 37, no. 2, pp. 165–176, 2006.

[67] S.T. Pinho, C.G. Davila, P.P. Camanho, L. Iannucci and P. Robinson, “Failuremodels and criteria for FRP under in-plane or three-dimensional stress statesincluding shear non-linearity,” NASA/TM-2005-213530, 2005.

154

Bibliography

[68] S.T. Pinho, L. Iannucci and P. Robinson, “Physically-based failure models andcriteria for laminated fibre-reinforced composites with emphasis on fibre kink-ing: Part I: Development,” Composites Part A: Applied Science and Manufac-turing, vol. 37, no. 1, pp. 63–73, 2006.

[69] S.T. Pinho, L. Iannucci and P. Robinson, “Physically based failure models andcriteria for laminated fibre-reinforced composites with emphasis on fibre kink-ing. Part II: FE implementation,” Composites Part A: Applied Science andManufacturing, vol. 37, no. 5, pp. 766–777, 2006.

[70] R. Olsson, “Experimental observations on the orthotropic transverse strength ofNon-Crimp Fabric composites,” Tech. Rep. TN15001, 2015. Swerea SICOMP.

[71] S.T. Pinho, R. Darvizeh, P. Robinson, C. Schuecker and P.P. Camanho, “Ma-terial and structural response of polymer-matrix fibre-reinforced composites,”Journal of Composite Materials, vol. 46, no. 19-20, pp. 2313–2341, 2012.

[72] E.K Tschegg, K. Humer and H.W. Weber, “Fracture tests in Mode I on fibre-reinforced plastics,” Journal of Materials Science, vol. 28, no. 9, pp. 2471–2480,1993.

[73] K. Humer, E.K. Tschegg, S.M. Spiessberger and H.W. Weber, “Fracture mech-anical characterization of fibre reinforced plastics in the intralaminar crackopening mode (mode I) and in the shear mode (mode II),” Cryogenics, vol. 38,no. 1, pp. 69–74, 1998.

[74] S. Kamiya and H. Sekine, “Prediction of the fracture strength of notched con-tinuous fiber-reinforced laminates by interlaminar crack extension analysis,”Composites Sci.ence and Technology, vol. 56, no. 1, pp. 11–21, 1996.

[75] S.T. Pinho, P. Robinson and L. Iannucci, “Fracture toughness of the tensile andcompressive fibre failure modes in laminated composites,” Composites Scienceand Technology, vol. 66, no. 13, pp. 2069–2079, 2006.

[76] G. Catalanotti, P.P. Camanho, J. Xavier, C.G. Davila and A.T. Marques,“Measurement of resistance curves in the longitudinal failure of compositesusing digital image correlation,” Composites Sci.ence and Technology, vol. 70,no. 13, pp. 1986–1993, 2010.

[77] J. Masters, “Translaminar Fracture Toughness of a Composite Wing Skin Madeof Stitched Warp-knit Fabric,” NASA/CR-201728, 1997.

[78] N. Blanco, D. Trias, S.T. Pinho and P. Robinson, “Intralaminar fracture tough-ness characterisation of woven composite laminates. Part II: Experimental char-acterisation,” Engineering Fracture Mechanics, vol. 131, pp. 361–370, 2014.

155

Bibliography

[79] R. Teixeira, Translaminar fracture toughness of CFRP: from the toughness ofindividual plies to the toughness of the laminate. PhD thesis, Imperial CollegeLondon, 2015.

[80] R. Olsson, “Analytical models for impact response and damage growth in sand-wich panels,” Tech. Rep. CR13-011, 2013. Swerea SICOMP.

[81] T.B Block, C. Brauner, M.I Zuardy and A.S. Herrmann, “Advanced numericalinvestigation of the impact behaviour of CFRP foam core sandwich structures,”in Proceedings of the 3rd ECCOMAS Thematic Conference on Mechanical Re-sponse of Composites, 2011.

[82] M.J. Laffan, S.T. Pinho, P. Robinson and A.J. McMillan, “Translaminar frac-ture toughness testing of composites: A review ,” Polymer Testing, vol. 31,no. 3, pp. 481–489, 2012.

[83] M.J. Laffan, S.T. Pinho, P. Robinson, A.J. McMillan, “Translaminar fracturetoughness: The critical notch tip radius of 0◦ plies in CFRP,” CompositesScience and Technology, vol. 72, no. 1, pp. 97–102, 2011.

[84] M.J. Laffan, S.T. Pinho, P. Robinson and L. Iannucci, “Measurement of the insitu ply fracture toughness associated with mode I fibre tensile failure in FRP.Part I: Data reduction,” Composites Science and Technology, vol. 70, no. 4,pp. 606–613, 2010.

[85] Dassault Systemes Simulia Corp., ABAQUS 6.12, User’s Manual. 2012.

[86] W.C. Jackson and J.G. Ratcliffe, “Measurement of fracture energy for kink-band growth in sandwich specimens,” in Proceedings of the 2nd InternationalConference on Composites Testing and Model Identification (CompTest 2004),2004.

[87] C.G. Davila, C.A. Rose and P.P. Camanho, “A procedure for superposing lin-ear cohesive laws to represent multiple damage mechanisms in the fracture ofcomposites,” International Jorunal of Fracture, vol. 158, no. 2, pp. 211–223,2009.

[88] M.J. Laffan, S.T. Pinho, P. Robinson and L. Iannucci, “Measurement of the insitu ply fracture toughness associated with mode I fibre tensile failure in FRP.Part II: Size and lay-up effects,” Composites Science and Technology, vol. 70,no. 4, pp. 614–621, 2010.

[89] R.W. McCune, C.G. Armstrong and D.J. Robinson, “Mixed-dimensional coup-ling in finite element models,” International Journal of Numerical Methods inEngineering, vol. 49, no. 6, pp. 725–750, 2000.

156

Bibliography

[90] C.G. Davila, “Solid-to-shell transition elements for the computation of inter-laminar stresses,” Computing Systems in Engineering, vol. 5, no. 2, pp. 193–202,1994.

[91] Z.P. Bazant, “Spurious reflection of elastic waves in non-uniform finite elementgrids,” Computer Methods in Applied Mechanics and Engineering, vol. 16, no. 1,pp. 91–100, 1978.

[92] Z. Celep and D. Turhan, “Transient wave propagation in constant and linearstrain finite elements,” Journal of Sound and Vibration, vol. 116, no. 1, pp. 15–23, 1987.

[93] J.N. Reddy, Mechanics of laminated composite plates and shells: theory andanalysis. CRC press, 2004.

[94] I. Kreja, “A literature review on computational models for laminated compositeand sandwich panels,” Open Engineering, vol. 1, no. 1, pp. 59–80, 2011.

[95] J.N. Reddy, “An evaluation of equivalent-single-layer and layerwise theories ofcomposite laminates,” Composite Structures, vol. 25, no. 1, pp. 21–35, 1993.

[96] W.C. Chao and N.J. Reddy, “Analysis of laminated composite shells usinga degenerated 3D element,” International Journals of Numerical Methods inEngineering, vol. 20, no. 11, pp. 1991–2007, 1984.

[97] E.J. Barbero, “A 3-D finite element for laminated composites with 2-D kin-ematic constraints,” Computers & Structures, vol. 45, no. 2, pp. 263–271, 1992.

[98] Dassault Systemes Simulia Corp., ABAQUS 6.12, Theory Manual. 2012.

[99] K.R. Shim, D.J. Monaghan and C.G. Armstrong, “Mixed dimensional couplingin finite element stress analysis,” Engineering with Computers, vol. 18, no. 3,pp. 241–252, 2002.

[100] N. Osawa, K. Hashimoto, J. Sawamura, T. Nakai and S. Suzuki, “Study onshell-solid coupling FE analysis for fatigue assessment of ship structure,” MarineStructures, vol. 20, no. 3, pp. 143–163, 2007.

[101] R.D. Mindlin, “Influence of rotary inertia and shear on flexural motions ofisotropic, elastic plates,” Journal of Applied Mechanics, vol. 18, pp. 31–38,1951.

[102] E. Reissner, “The effect of transverse shear deformation on the bending of elasticplates,” Journal of Applied Mechanics, vol. 12, pp. 69–77, 1945.

157

Bibliography

[103] A. Tessler, “A higher-order plate theory with ideal finite element suitabil-ity,” Computer Methods in Applied Mechanics and Engineering, vol. 85, no. 2,pp. 183–205, 1991.

[104] J. Reinoso, A. Blazquez, A. Estefan, F. Parıs, J. Canas, E. Arevalo and F.Cruz, “Experimental and three-dimensional global-local finite element analysisof a composite component including degradation process at the interfaces,”Composites Part B: Engineering, vol. 43, no. 4, pp. 1929–1942, 2012.

[105] J. Reinoso, A. Blazquez, A. Estefani, F. Paris and J. Canas, “A compositerunout specimen subjected to tension-compression loading conditions: Experi-mental and global–local finite element analysis,” Composite Structures, vol. 101,pp. 274–289, 2013.

[106] R.Krueger and T.K. O’Brien, “A shell/3D modeling technique for the analysisof delaminated composite laminates,” Composites Part A: Applied Science andManufacturing, vol. 32, no. 1, pp. 25–44, 2001.

[107] R.Krueger, J. Ratcliffe and P.J. Minguet, “Analysis of composite panel-stiffenerdebonding using a shell/3D modeling technique, NASA/CR-2007-214879 NIAReport No. 2007-07,” 2007.

[108] R. Krueger and P.J. Minguet, “Analysis of composite skin-stiffener debondspecimens using a shell/3D modeling technique,” Composite Structures, vol. 81,no. 1, pp. 41–59, 2007.

[109] R. Krueger, J.C. Ratcliffe and P.J. Minguet, “Panel stiffener debonding ana-lysis using a shell/3D modeling technique,” Composites Science and Technology,vol. 69, no. 14, pp. 2352–2362, 2009.

[110] T. Belytschko and S.P. Xiao, “Coupling methods for continuum model withmolecular model,” International Journal for Multiscale Computational Engin-eering, vol. 1, no. 1, pp. 115–126, 2003.

[111] S.P. Xiao and T. Belytschko, “A bridging domain method for coupling con-tinua with molecular dynamics,” Computer Methods in Applied Mechanics andEngineering, vol. 193, no. 17, pp. 1645–1669, 2004.

[112] N. Bicanic, Discrete Element Methods. Encyclopedia of Computational Mech-anics, John Wiley & Sons (Ltd), 2004.

[113] J. Rojek and E. Onate, “Multiscale analysis using a coupled discrete/finiteelement model,” Interaction and Multiscale Mechanics, vol. 1, no. 1, pp. 1–31,2007.

158

Bibliography

[114] Y. Lee and C. Basaran, “A multiscale modeling technique for bridging molecu-lar dynamics with finite element method,” Journal of Computational Physics,vol. 253, pp. 64–85, 2013.

[115] F. Aymerich, C. Pani and P. Priolo, “Damage response of stitched cross-plylaminates under impact loadings,” Engineering Fracture Mechanics, vol. 74,no. 4, pp. 500–514, 2007.

[116] F. Aymerich, F. Dore and P. Priolo, “Prediction of impact-induced delaminationin cross-ply composite laminates using cohesive interface elements,” CompositesScience and Technology, vol. 68, no. 12, pp. 2383–2390, 2008.

[117] W-G. Jiang, S.R. Hallett and M.R. Wisnom, “Development of domain super-position technique for the modelling of woven fabric composites,” in Mechanicalresponse of composites, pp. 281–291, Springer, 2008.

[118] E.V. Iarve, D.H. Mollenhauer, E.G. Zhou, T. Breitzman and T.J. Whitney,“Independent mesh method-based prediction of local and volume average fieldsin textile composites,” Composites Part A: Applied Science and Manufacturing,vol. 40, no. 12, pp. 1880–1890, 2009.

[119] S.A. Tabatabaei, S.V. Lomov and I. Verpoest, “Assessment of embedded ele-ment technique in meso-FE modelling of fibre reinforced composites,” Compos-ite Structures, vol. 107, no. 0, pp. 436–446, 2014.

[120] S.V. Lomov, D.S. Ivanov, I. Verpoest, M. Zako, T. Kurashiki, H. Nakai and S.Hirosawa, “Meso-FE modelling of textile composites: Road map, data flow andalgorithms,” Composites Science and Technology, vol. 67, no. 9, pp. 1870–1891,2007.

[121] S. Osher and J.A. Sethian, “Fronts propagating with curvature-dependentspeed: algorithms based on Hamilton-Jacobi formulations,” Journal of Compu-tational Physics, vol. 79, no. 1, pp. 12–49, 1988.

[122] S. Abrate, Impact on composite structures. Cambridge University Press, 2005.

[123] J. LLorca, C. Gonzalez, J.M. Molina-Aldareguıa, J. Segurado, R. Seltzer, F.Sket, M. Rodrıguez, S. Sadaba, R. Munoz and L.P. Canal, “Multiscale mod-eling of composite materials: a roadmap towards virtual testing,” AdvancedMaterials, vol. 23, no. 44, pp. 5130–5147, 2011.

[124] Dassault Systemes Simulia Corp., ABAQUS 6.12, Scripting Reference Manual.2012.

159

Bibliography

[125] G. Alfano and M.A. Crisfield, “Finite element interface models for the delamin-ation analysis of laminated composites: mechanical and computational issues,”International Journal of Numerical Methods in Engineering, vol. 50, no. 7,pp. 1701–1736, 2001.

[126] Q. Yang and B. Cox, “Cohesive models for damage evolution in laminatedcomposites,” International Journal of Fracture, vol. 133, no. 2, pp. 107–137,2005.

[127] A. Gravouil and A. Combescure, “Multi-time-step explicit–implicit method fornon-linear structural dynamics,” International Journal of Numerical Methodsin Engineering, vol. 50, no. 1, pp. 199–225, 2001.

[128] A. Combescure and A. Gravouil, “A numerical scheme to couple subdomainswith different time-steps for predominantly linear transient analysis,” ComputerMethods in Applied Mechanics and Engineering, vol. 191, no. 11, pp. 1129–1157,2002.

[129] C. Farhat and F.X. Roux, “A method of finite element tearing and intercon-necting and its parallel solution algorithm,” International Journal of NumericalMethods in Engineering, vol. 32, no. 6, pp. 1205–1227, 1991.

[130] C. Farhat, K. Pierson and M. Lesoinne, “The second generation FETI methodsand their application to the parallel solution of large-scale linear and geomet-rically non-linear structural analysis problems,” Computer Methods in AppliedMechanics and Engineering, vol. 184, no. 2-4, pp. 333–374, 2000.

[131] M. Brun, A. Batti, A. Limam and A. Gravouil, “Explicit/implicit multi-timestep co-computations for blast analyses on a reinforced concrete frame struc-ture,” Finite Elements in Analysis and Design, vol. 52, pp. 41–59, 2012.

[132] M. Brun, A. Batti, A. Limam and A. Combescure, “Implicit/explicit multi-time step co-computations for predicting reinforced concrete structure responseunder earthquake loading,” Soil Dynamics and Earthquake Engineering, vol. 33,no. 1, pp. 19–37, 2012.

[133] E. Duni, G. Toniato, R. Saponaro, P. Smeriglio and V. Puleo, “Vehicle dynamicsolution based on finite element tire/road interaction implemented through Im-plicit/Explicit sequential and Co-Simulation approach,” SAE Technical Paper2010-01-1138, 2010.

[134] H. Surendranath,V. Oancea and S. Subbarayalu, “Full vehicle durability predic-tion using Co-Simulation between Implicit and Explicit finite element solvers,”tech. rep., Dassault Systemes Simulia Corp., 2011.

160

Bibliography

[135] M. Brun, A. Batti, A. Combescure and A. Gravouil, “External coupling softwarebased on macro-and micro-time scales for explicit/implicit multi-time-step co-computations in structural dynamics,” Finite Elements in Analysis and Design,vol. 86, pp. 101–119, 2014.

[136] L. Gigliotti and S.T. Pinho, “Multiple length/time-scale simulation of local-ized damage in composite structures using a Mesh Superposition Technique,”Composite Structures, vol. 121, pp. 395–405, 2015.

[137] A. Gusev, “Representative volume element size for elastic composites: a nu-merical study,” Journal of Mechanics of Physics and Solids, vol. 45, no. 9,pp. 1449–1459, 1997.

[138] I.M. Gitman, H. Askes and L.J. Sluys, “Representative volume: existence andsize determination,” Engineering Fracture Mechanics, vol. 74, no. 16, pp. 2518–2534, 2007.

[139] C. Sun and R. Vaidya, “Prediction of composite properties from a representativevolume element,” Composites Science Technology, vol. 56, no. 2, pp. 171–179,1996.

[140] Z. Xia, Y. Zhang and F. Ellyin, “A unified periodical boundary conditions forrepresentative volume elements of composites and applications,” InternationalJournal of Solids and Structures, vol. 40, no. 8, pp. 1907–1921, 2003.

[141] Z. Xia, C. Zhou, Q. Yong and X. Wang, “On selection of repeated unit cell modeland application of unified periodic boundary conditions in micro-mechanicalanalysis of composites,” International Journal of Solids and Structures, vol. 43,no. 2, pp. 266–278, 2006.

[142] S. Li, “Boundary conditions for unit cells from periodic microstructures andtheir implications,” Composites Science and Technology, vol. 68, no. 9, pp. 1962–1974, 2008.

[143] S. Li, “General unit cells for micromechanical analyses of unidirectional com-posites,” Composites Part A: Applied Science and Manufacturing, vol. 32, no. 6,pp. 815–826, 2001.

[144] S. Li and A. Wongsto, “Unit cells for micromechanical analyses of particle-reinforced composites,” Mechanics of Materials, vol. 36, no. 7, pp. 543–572,2004.

[145] J.D. Whitcomb, C.D. Chapman and X. Tang, “Derivation of boundary condi-tions for micromechanics analyses of plain and satin weave composites,” Journalof Composite Materials, vol. 34, no. 9, pp. 724–747, 2000.

161

Bibliography

[146] N.V. de Carvalho, S.T. Pinho and P. Robinson, “Reducing the domain in themechanical analysis of periodic structures, with application to woven compos-ites,” Composites Science and Technology, vol. 71, no. 7, pp. 969–979, 2011.

[147] N.V. de Carvalho, S.T. Pinho and P. Robinson, “A mathematical frameworkfor reducing the domain in the mechanical analysis of periodic structures,”arXiv:1012.3133v1, pp. 1–18, 2010.

[148] X. Peng and J. Cao, “A dual homogenization and finite element approach formaterial characterization of textile composites,” Composites Part B: Engineer-ing, vol. 33, no. 1, pp. 45–56, 2002.

[149] J. Llorca, A. Needleman and S. Suresh, “An analysis of the effects of matrixvoid growth on deformation and ductility in metal-ceramic composites,” ActaMetallica Materialia, vol. 39, no. 10, pp. 2317–2335, 1991.

[150] J. Yeh, “The effect of interface on the transverse properties of composites,”International Journal of Solids and Structures, vol. 29, no. 20, pp. 2493–2502,1992.

[151] Z. P. Bazant, “Can multiscale-multiphysics methods predict softening damageand structural failure?,” International Journal of Multiscale Computational En-gineering, vol. 8, no. 1, 2010.

[152] R. M. Jones, Mechanics of composite materials. CRC Press, 1998.

[153] HexFlow, “RTM6-2 - Product Information..” http://www.hexcel.com/Resources/DataSheets/RTM-Data-Sheets/RTM6-2_eu.pdf, 2010.

[154] Evonik Industries AG, “ROHACELL HERO - Technical Information.” http://www.rohacell.com/sites/dc/Downloadcenter/Evonik/Product/ROHACELL/product-information/ROHACELL%20HERO%20Product%20Information.pdf,2014.

[155] R. Olsson, “Private communications.” Swerea SICOMP AB, 2014.

[156] T. Bru, “Shear tests of carbon non-crimp fabric composites used in sandwichskins,” Tech. Rep. CR14-069, Swerea SICOMP, 2014.

[157] HexFlow, “IM7 Carbon Fiber - Product Information..” http://www.hexcel.com/resources/datasheets/carbon-fiber-data-sheets/im7.pdf, 2014.

[158] G. Endres, “Innovative sandwich constructions for aircraft application,” in Pro-ceedings of the 9th International Conference on Sandwich Structures, Pasadena(CA), 2010.

162

Bibliography

[159] T. Liu, Z.C. Deng and T.J. Lu, “Analytical modeling and finite element simula-tion of the plastic collpase of sandwich beams with pin-reinforced foam cores,”International Journal of Solids and Structures, vol. 45, pp. 5127–5151, 2008.

[160] S. Cron, “Finite Element Analysis of a Coil Spring RFID Antenna Embeddedin a Rubber Patch,” in 2007 ABAQUS UsersâĂŹ Conference, 2007.

[161] J. Hohe, “A direct homogenisation approach for determination of the stiff-ness matrix for microheterogeneous plates with applcation to sandwich panels,”Composites Part B: Engineering, vol. 34, pp. 615–626, 2003.

[162] J. Gager and H.E. Pettermann, “Numerical homogenization of textile compos-ites based on shell element discretization,” Composites Science and Technology,vol. 72, no. 7, pp. 806–812, 2012.

[163] D.M. Thompson and O.H. Griffin, “Two-dimensional/three-dimensionalglobal/local finite element analysis of composite laminates,” in Proceedings ofthe 4th Annual American Society for Composites Conference, Lancaster, PA(T. Publishing, ed.), pp. 234–243, 1989.

[164] L. Gigliotti and S.T. Pinho, “Exploiting symmetries in solid-to-shell homo-genization, with application to periodic pin-reinforced sandwich structures,”Composite Structures, vol. 132, pp. 995–1005, 2015.

[165] R.T. Whitcomb, “Review of NASA supercritical airfoils,” in Proceediings of the9th International Council of the Aeronautical Sciences Congress, Haifa, Israel,1974.

[166] L. Gigliotti and S.T. Pinho, “Prediction of the post-crushing compressive re-sponse of progressively crushable sandwich foam cores,” Composites Part A:Applied Science and Manufacturing, vol. 80, pp. 148–158, 2016.

[167] R. Olsson and T.B. Block, “Criteria for skin rupture and core shear crackinginduced by impact on sandwich panels,” Composite Structures, vol. 125, pp. 81–87, 2015.

[168] V. Deshpande and N. Fleck, “Isotropic constitutive models for metallic foams,”Journal of the Mechanics and Physics of Solids, vol. 48, no. 6, pp. 1253–1283,2000.

[169] R. Courant, K. Friedrichs and H. Lewy, “On the partial difference equationsof mathematical physics,” IBM Journal of Research and Development, vol. 11,no. 2, pp. 215–234, 1967.

163

Appendix A

Propagation of stress waves innon-uniform FE meshes usingthe MST

A.1 Motivation

The benefits provided by the MST in avoiding the artificial stress disturbances atthe interfaces between differently-discretized subdomains are demonstrated in [136].In this example, we aim to investigate the capability of the MST to mitigate thespurious reflections of stress waves at the interfaces between differently-discretizedsubdomains.

A.2 FE models

An infinite bar of uniform cross-sectional area (4×4 mm2) and made of an isotropicelastic material (E =60 MPa, ν =0.25 and ρ=1600 kg/m3) was modelled in Abaqus.The initial 300 mm of the bar were discretized using conventional finite elements;the infinitely long idealization is achieved using infinite elements (CIN3D8). Forthe region were conventional FE models were used, three different FE models of theinfinite bar were created (see Figure A.1):

Fully Solid (FS) modelthe bar is discretized entirely using eight-noded reduced-integration solid ele-ments (C3D8R) with enhanced hourglass control. Results obtained with thismodel are used as reference solution;

165

Appendix A

Sudden Transition (ST) modelthe bar is discretized using eight-noded reduced-integration solid elements(C3D8R) with enhanced hourglass control for the initial 150 mm; the remain-ing portion, eight-noded reduced-integration continuum shell elements (SC8R),with hourglass control and three integration points along the element’s thick-ness, were used;

MST modelthe sudden discretization-transition between solid (C3D8R) and continuum shell(SC8R) elements is replaced by a region where the MST is exploited. The lengthof the MST region is equal to 100 mm.

Fully Solid Fully Solid

MST model

Sudden Transition model

Sudden discretization-transition

MST region

Region 1 Region 2

1

1

1x

x

x

Figure A.1: Three different FE models of the infinite bar: Fully Solid model, SuddenTransition model and MST model.

In Figure A.1, Region 1 corresponds to the portion of the bar between x=100 mmand x = 200 mm, while Region 2 corresponds to the portion of the bar betweenx = 200 mm and x = 300 mm. Region 1 includes the sudden distretization-transitionfor the ST model, and the entire MST region for the MST model.

A.3 Stress-wave propagation analysis

The applied displacement and boundary conditions considered are shown in Fig-ure A.2. The applied displacement consists of a finite discrete impulse defined usingFourier series as

u(z, t) =u(z)

κ

κ∑i=1

sin(2πifmint) , (A.1)

where κ =fmax

fminis the number of terms in the Fourier series, u(z) is the amplitude of

the applied displacement as a function of the out-of-plane coordinate, and fmin andfmax are, respectively, the minimum and maximum frequency of interest.

166

Propagation of stress waves in non-uniform FE meshes using the MST

4 mm

4 mm

y

z

u (z ; t)

y

x

300 mm , 300 elem ents

1

C onventiona l F in ite E lem ents

Figure A.2: Applied displacement and boundary conditions for the stress-wavepropagation study.

In Figure A.3, the evolution of the applied displacement as a function of time(t) Figure A.3a and of the z-coordinate Figure A.3b is shown. Simulations wereperformed using Abaqus/Standard. The analysis time-step was chosen so that theCourant-Friedrichs-Lewy (CFL) condition [169] is fulfilled. The essence of the CFLcondition is that, given a certain wave travelling across a discrete spatial grid, itsamplitude can be correctly computed at fixed time steps of constant duration, onlyif the latter is smaller than the time needed for the wave to travel between two con-secutive nodes of the spatial grid. In the simplified scenario of 1D wave propagation,

0

0.25

0.5

0.75

1

0E+0 4E-3 8E-3

0

0.05

0.1

0E+0 5E-4 1E-3

·Xi= 1

s in (2¼ ifm in t )

·[¡ ]

t [m s]

(a) Normalized amplitude of the applied dis-placement as a function of time.

0

0.25

0.5

0.75

1

0 1 2 3 4

z [m m ]

u (z ) [mm ]

(b) Amplitude of the applied displacement asa function of the z-coordinate.

Figure A.3: Time and space evolution of the applied displacement u(z, t). Theresulting applied displacement is a finite discrete impulse with variableamplitude along the through-the-thickness direction.

167

Appendix A

the CFL condition can be expressed in terms of the Courant number C as

C =cΔtmax

≤ Cmax , (A.2)

where c is the dilatational wave speed and is the characteristic length of the spa-tial grid. A value of Cmax = 1 is usually tolerated when implicit solvers are used.Nevertheless, a maximum time-step of 20 ns (C≈ 0.122) was chosen for the analysispresented here.

A.4 Results and discussion

Figure A.4 shows the internal energy U of Region 1 and Region 2 (see Figure A.1)as a function of time. Such internal energy corresponds to the energy carried by thestress waves travelling along the bar.

The stress waves enter into Region 1 at t ≈ 0.014 ms when the energy contentof the latter starts increasing Figure A.4a. At t ≈ 0.029 ms, the stress waves travelfrom Region 1 into Region 2 (the energy content of Region 1 starts decreasing and,correspondingly, that of Region 2 starts increasing).

However, after t≈ 0.029 ms the energy content of Region 1 computed with the STmodel is higher than that computed with both the FS and MST model; a specular

0

10

20

30

0 0.02 0.04

SOL

ST

MST

U [J ]

t [m s]

F S

ST

M ST

(a) Region 1.

0

10

20

30

0 0.02 0.04

SOL

ST

MST

U [J ]

t [m s]

F S

ST

M ST

(b) Region 2.

Figure A.4: Time history of the internal energy of Region 1 (a) and Region 2 (b)associated with the propagation of the travelling stress waves along theinfinite bar.

168

Propagation of stress waves in non-uniform FE meshes using the MST

behaviour is associated to the energy content of Region 2. This evidence indicatesthat part of the stress waves travelling along the bar has been reflected by the suddendiscretization-transition included in the ST model.

When compared to the reference solution obtained with the FS model, an energyloss of approximately 17 % due to stress-wave reflections is observed when using theST model. On the other hand, when using the MST for the progressive transitionbetween the two differently-discretized subdomains, the energy loss associated tostress-wave reflections is almost completely eliminated (approximately 99 % reductionif compared to the ST model).

The results presented in this section demonstrate that the MST [136] presentedin Chapter 4 allows to avoid, in addition to the artificial stress disturbances, also thespurious reflections of stress waves at the interfaces between differently-discretizedsubdomains.

169

Gigliotti-L-2016-PhD-Thesis (2)

Documents

Transcript of Gigliotti-L-2016-PhD-Thesis (2)