DF Development of modelling and homogenisation procedures for stochastic tape-based discontinuous composites Master’s thesis in Applied Mechanics Olle Haglund Nilsson Jacob Sjöberg DEPARTMENT OF INDUSTRIAL AND MATERIALS SCIENCE Chalmers University of Technology Gothenburg, Sweden 2022 Master’s thesis in applied mechanics 2022 Development of modelling and homogenisation procedures for stochastic tape-based discontinuous composites OLLE HAGLUND NILSSON JACOB SJÖBERG DF Department of Industrial and Materials Science Division of Material and Computational Mechanics Chalmers University of Technology Gothenburg, Sweden 2022 Development of modelling and homogenisation procedures for stochastic tape-based discontinuous composites © Olle Haglund Nilsson, Jacob Sjöberg, 2022. Supervisor: Prof. Martin Fagerström, Division of Material and Computational Mechanics. Examiner: Prof. Martin Fagerström, Division of Material and Computational Mechanics. Master’s Thesis 2022:07 Department of Industrial and Materials Science Division of of Material and Computational Mechanics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Generated geometrical model of an stochastic tape-based discontinuous composite. Typeset in LATEX Printed by Chalmers Reproservice Gothenburg, Sweden 2022 Development of modelling and homogenisation procedures for stochastic tape-based discontinuous com- posites OLLE HAGLUND NILSSON JACOB SJÖBERG Department of Industrial and Materials Science Division of Material and Computational Mechanics Chalmers University of Technology Abstract Conventional laminated composites have, due to their high stiffness- and strength-to-weight ratios, shown to be a viable alternative to metals in applications where lightweight structures are of essence. However, the current reach of conventional laminated composites is limited to certain industries and high-end products as a result of their high manufacturing cost. As a way to reduce cost while maintaining high performance, stochastic tape-based discontinuous composites (STBDCs) have emerged, offering a middle ground between the easy manufacturing of short fibre composites and high performance of continuous fibres composites. The irregular mesostructure and high interest from industries have increased the de- mand for efficient and predictive models (analytical and numerical) describing the mechanical response of this class of material. As a response to the increased interest, this project aimed to develop a method to predict the elastic properties of STBDCs, based on the mesostructural configuration and constituent properties (impregnated tapes and pure matrix). To replicate the complex material structure, a modular MATLAB code was developed as the core of the project. The MATLAB code takes the tape and plate dimensions, material properties, voxel resolution, etc, as inputs, then builds the material layer-by-layer by placing tapes at random positions and angles until the desired plate thickness has been reached. The developed algorithm allows tapes to form around each other (drape) to replicate the compression part of the real manufacturing process. The outputs of the code are, among other visualisations, a 3D-image of the material structure and Abaqus input files for the coming analysis. Based on the generated material structures, smaller statistical volume elements (SVEs) were extracted and used for the analysis. Nearly 400 SVEs of different dimensions from different samples were extracted and analysed to have a large enough sample size to account for the variability of the material. Computational homogenisation was used to determine the volume averaged in-plane elastic parameters. The homogenisation process was carried out by applying periodic boundary conditions (PBCs) to the SVEs using the Abaqus plug-in EasyPBC. The computational homogenisation was initially focused on a full 3D-model but as a subsequent step, the full 3D-model was reduced to a 2D-model by an intermediate analytical homogenisation process using classical laminate theory. A full 3D- and reduced 2D-model were generated for each SVE sample to allow a comparative analysis. Finally, experimental data of manufactured STBDCs plates was used as a comparison to verify the results of the numerical model. The study found that the voxel-based 3D-model and the developed methodology can be used to ac- curately evaluate the elastic properties of STBDCs. More specifically, the generated material samples and methodology used provided reliable results for all studied SVE dimensions, converging to elastic properties within one standard deviation compared to experimental tests. Furthermore, the reduced 2D-model showed a similar accuracy compared to the 3D-model while also requiring significantly less computational power, indicating that this is a reasonable simplification to make. It should be noted that the developed model is not a perfect representation of reality, simplifications had to be made to stay within the limitations of the project. As a consequence, the fibre volume fraction (FVF) of experimental data could not be reached, making the model slightly under-predict the elas- tic properties. With further improvements to raise the FVF, the data indicates that the model would produce more accurate results compared to experimental tests, thus being a reliable source material for future industrial use. Keywords: Stochastic Tape Based Discontinuous Composite, Periodic Boundary Condition, Statistical Volume Element, Voxel, Chopped Fibre Composite, Tow, Homogenisation, i ii Preface The master’s thesis was conducted as a 30 credits project at Chalmers University of Technology during spring of 2022. Both authors, Olle Haglund Nilsson and Jacob Sjöberg were students at the master’s programme Applied Mechanics with previous Bachelor of science degrees in Mechanical Engineering. Acknowledgements This thesis work was conducted at the Division of Material and Computational Mechanics at Chalmers University of Technology. We would first and foremost like to thank our supervisor and examiner Professor Martin Fagerström, without you this project would have been only a fraction of what it became. Your expertise, experience and guidance has been an incredible help throughout the project. Furthermore, we want to express our gratitude to everyone at the division. Thank you for opening your arms to us, making us feel like a part of the division. Thanks to all of you, the time during our thesis has not only been a great learning experience, but also a genuinely enjoyable time that we will cherish moving forward. iii Table of Contents Nomenclature vii 1 Introduction 1 1.1 Aim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Previous work 4 2.1 Experimental testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Digital modelling and finite element analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 Theory 7 3.1 Stochastic tape-based discontinuous composites . . . . . . . . . . . . . . . . . . . . . . . . 7 3.1.1 Manufacturing process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.1.2 Microstructure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.1.3 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2 Geometrical modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Finite element analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3.1 Discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3.2 Homogenisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3.3 Periodic boundary conditions (PBCs) . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3.4 Evaluation of effective material properties . . . . . . . . . . . . . . . . . . . . . . . 17 3.4 Classical laminate theory - analytical homogenisation . . . . . . . . . . . . . . . . . . . . . 18 4 Methodology 19 4.1 Geometrical modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.1.1 Definition of domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.1.2 Tape generation and placement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.1.3 Draping algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1.4 Geometric material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.2 Finite element analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2.1 Discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2.2 Implementation of periodic boundary conditions . . . . . . . . . . . . . . . . . . . 26 4.2.3 Evaluation of elastic properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.2.4 Abaqus input file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.3 2D-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4 Design of experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5 Results 32 5.1 Geometrical properties of the numerical STBDC model . . . . . . . . . . . . . . . . . . . 32 5.1.1 Results of domain modifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 5.1.2 Mesostructure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.1.3 Microstructure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.1.4 Distribution of fibre volume fraction and average tape angle . . . . . . . . . . . . . 36 5.2 Homogenised material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 5.2.1 Convergence study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.2.2 Stress and strain response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.3 Material properties versus geometrical properties . . . . . . . . . . . . . . . . . . . . . . . 44 5.3.1 Young’s modulus versus global average tape angle . . . . . . . . . . . . . . . . . . 44 5.3.2 Shear modulus versus Young’s modulus . . . . . . . . . . . . . . . . . . . . . . . . 46 5.3.3 Young’s modulus versus fibre volume fraction . . . . . . . . . . . . . . . . . . . . . 46 5.4 Application to a practical structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 iv 6 Discussion 51 6.1 Geometrical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.1.1 Deforming and sliding tapes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 6.1.2 Tape-resolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.1.3 Selective tape placement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.1.4 Removal of top layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 6.1.5 Under-stiff draping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.2 Homogenised material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.2.1 Scaling of input material parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.2.2 Variation between SVE sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.2.3 Reduced integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.2.4 Convergence of Young’s moduli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6.2.5 Appearance of tape bands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 6.3 Practical implementation of homogenised material properties . . . . . . . . . . . . . . . . 55 7 Conclusions 56 8 Future work 57 v Nomenclature Standard variables ν Poisson’s ratio σ Stress ε Strain φ In-plane tape angle D Plate volume domain E Young’s modulus G Shear modulus l Length N t vox,h Through thickness number of voxels (tape) N t vox,l Number of length voxels (tape) N t vox,w Number of width voxels (tape) Rvox Voxel aspect ratio t Thickness u Displacement V Volume Vf Fibre volume fraction w Width Xc Tape centre position (x-direction) Yc Tape centre position (y-direction) Subscripts vox Voxel c Centre L Longitudinal T Transverse Superscripts CA Cumulative average e Element p Plate t Tape Abbreviations FEA Finite Element Analysis vi FVF Fibre Volume Fraction PBC Periodic Boundary Condition RP Reference Point RSA Random Sequential Adsorption RVE Representative Volume Element STBDC Stochastic Tape Based Discontinuous Composite SVE Statistical Volume Element (aperiodic RVE) vii 1 Introduction Conventional laminated composites are an advantageous alternative to metals in numerous applications, especially because compared to metals these materials offer better stiffness- and strength-to-weight ra- tios, corrosion- and fatigue resistance [1]. Given their properties and material structure, continuous composites are particularly suitable for long continuous structures from a manufacturability perspective. Currently these materials are more commonly used in the aircraft industry than in the automotive in- dustry, although a similar carbon fibre reinforced polymer (CFRP) usage in the automotive industry as the aircraft industry could lead to a 54 Mt reduction of CO2 emissions per year due to weight reductions [1, 2]. For reference, a 10% reduction in weight of a personal vehicle leads to a 6-8% reduction in fuel consumption, or for electric vehicles, an increased range [3]. The use of carbon fibre composites in the automotive sector is low and limited to luxury vehicles due to their high manufacturing cost [2]. This project aims to investigate a material that due to its ease of manufacturing and mouldability introduce a viable alternative to metals that can extend the use of CFRP across all industries; stochastic tape-based discontinuous composites (STBDCs) [1]. STBDCs is a class of high performance lightweight discontinuous composites composed of randomly dis- tributed, polymer matrix impregnated, carbon fibre tapes, schematically shown in Figure 1(b) [2, 4]. The term carbon fibre tapes is used to describe flat strands of unidirectional carbon fibres, pre-impregnated with either a thermoplastic or thermoset polymer matrix [5]. The carbon fibre tapes may be sourced as residual scrap from the production of continuous fibre composites which leads to a significant reduction of waste [1] or directly from unidirectional carbon fibre sheets. Due to the discontinuity in the layout of the tapes, STBDCs can be compression moulded to complex geometries, a process where the tapes are placed and press-formed in a mould and then cured under pressure [4, 6]. Compression moulding gives STBDCs great potential for high volume production due to the short cycle time and low cost compared to conventional composites [2, 7]. STBDCs offer a middle ground between the easy manufacturing of short fibre composites and high per- formance of continuous fibre composites [8]. The key benefits of short fibre composites are that they are easy and fast to produce by either injection- or compression moulding and with that carry a low cost per unit [9]. The downside, however, is that they typically have only a fraction of the stiffness and strength of continuous fibre composites [9]. Continuous fibre composites are stiffer and stronger than short fibre composites but also more costly and difficult to manufacture into as complex shapes [1]. Discontinuous tape-based composites can due to their drapability and discontinuity be used to manufacture complex shapes while maintaining a low cost and high performance [1, 2, 4]. (a) Short fibre composites. effect on the performance of TBDCs, but this has not been quanti- tatively examined yet. A noticeable feature of TBDCs is that they are highly heteroge- neous, due to the large tow dimensions (up to 50 mm long and 10 mm wide) and random tow orientations. This results in loca- lised anisotropic properties, and in the presence of inherent weak spots (e.g. clusters of tow-ends, resin pockets or voids) in TBDCs, which leads to notch insensitivity [16]. These characteristics are also shared by other classes of Discontinuous Fibre-Composites (DFCs) [17–19], although the macroscopic effect of the heteroge- neous microstructure is reduced in the latter as the fibres are more dispersed in DFCs than in TBDCs. Models to predict the mechanical properties of TBDCs and of other DFCs commonly use an equivalent laminate analogy (Fig. 1), to account for the randomly-oriented architecture of the fibres or tows [2,20–23]. The equivalent laminate assumption rep- resents the random orientations of discontinuous composites into a ply-by-ply equivalent laminate, which contains unidirectional discontinuous plies at different orientations, usually forming a quasi-isotropic lay-up (see Fig. 1). Other models using orientation homogenisation methods may not explicitly mention the equiva- lent laminate to represent the microstructure of TBDCs and DFCs, but the mathematical formulation used to homogenise the modu- lus of the material is similar to the models directly considering an equivalent laminate [24–28]. Although the equivalent laminate assumption has been widely adopted, there is a striking lack of experimental validation compar- ing the mechanical performance of ply-by-ply equivalent lami- nates and that of actual discontinuous composites with a random architecture. Consequently, the applicability and accuracy of the equivalent laminate assumption to predict the mechanical perfor- mance of DFCs and TBDCs remains an open question. Therefore, this study aims to investigate the validity of the equivalent laminate assumption to model TBDCs, by comparing the mechanical properties of (i) actual TBDCs with a random microstructure and (ii) their equivalent laminates, for two different tow thicknesses. This is achieved by manufacturing and testing random TBDC plates and equivalent laminates, as described in Sec- tion 2. The testing results are presented in Section 3, and they are further discussed in Section 4; conclusions are presented in Section 5. 2. Experimental procedures 2.1. Materials HexPly–M77 [29] carbon-fibre/epoxy prepreg was used to man- ufacture all specimens analysed in this study. Two types of discon- tinuous composites were manufactured: 1. Random Composite (RC), made of chopped prepreg tows ran- domly oriented and distributed in a plate (see Fig. 2a). These RC plates are aimed to have similar microstructure and constituents to those of HexMC-M77 [30], a commercial version of TBDC. 2. Equivalent Laminate (EL), in which each lamina is a prepreg ply with discontinuous cuts defining aligned tows (see Fig. 2b). Two RC plates and two EL plates were made, using prepregs of either tow thickness tply ¼ 0:164 mm (later referred to as ‘thin- prepreg’) or tply ¼ 0:285 mm (later referred to as ‘thick-prepreg’). In all cases, the tow length (lt) and the tow width (wt) are lt ¼ 50 mm and wt ¼ 8 mm respectively. 2.2. Manufacturing of the random composite and equivalent laminate plates 2.2.1. Lay-up of random composites The number of tows in the thin-prepreg and thick-prepreg RC plates (Np tows, see Table 1) were selected to be exactly the same as in their corresponding EL plates. The tows at the edges of the RC plates were trimmed, and therefore the total volume of tows was slightly lower in the RC than in the EL; consequently, the nominal laminate thicknesses (tn) of the RC plates are slightly lower than their EL counterparts (see Table 1). To randomise the position and orientation of the tows, the following procedure was carried out: 1. Realising the position (xi; yi) and the orientation (hi) of each tow from uniform statistical distributions, using a randomisation algorithm (written in Matlab [31], see Fig. 2a and Appendix A). 2. Printing on paper several layers of the pattern generated by the randomisation algorithm on a one-to-one scale (see Fig. 2a), where each layer consisted of approximately 105 tows (which corresponds to half coverage of a uni-directional EL ply). 3. Cutting the prepreg to 50 mm � 8 mm tows using an auto- mated ply cutter, and manually placing the tows, one-by-one, on each printed layer according to the positions and orienta- tions defined in the random pattern generated in Step 2. The natural tackiness of the tows allowed them to stay in place rel- atively to each other. 4. Laying-up the layers (each consisting of approximately 105 ran- domly placed tows), so as to make RC plates with the required number of tows (defined in Table 1). The plates were cured as described in Section 2.2.3. 2.2.2. Lay-up of equivalent laminates The EL plates consist of uni-directional (UD) plies in a quasi- isotropic lay-up, where each ply consists of perfectly-aligned 50 mm � 8 mm tows (see Fig. 2b). The specifications of the EL plates are detailed in Table 1; classical laminate theory [28] was used to confirm that the lay-ups used had no coupling between tension and bending or torsion (the extensional stiffness matrix A and the bending matrix D of the lay-ups were ensured to be Fig. 1. Equivalent laminate assumption for modelling TBDCs. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Y. Li et al. / Composites: Part A 102 (2017) 64–75 65 (b) Layup of STBDCs, [4]. (c) Layup of con- ventional composite laminates. Figure 1: Schematic layups of (a) short fibre composite, (b) STBDCs, and (c) conventional composite laminates, fibre orientation visualised with black stripes. 1 Numerous companies globally are actively exploring the potential of implementing STBDCs into the design of new products, a selection of which are presented in Figure 2. Complex metal structures can easily be converted to use STBDC material with little design changes required, whereas large design modifications are typically needed to convert a complex metal part to use continuous laminate composites [10]. As a practical example, the suspension control arm for the Lamborghini concept car Sesto Elemento, shown in Figure 2(a), is a component under high loads. By replacing the aluminium with STBDC, the component became on average 27 % lighter and had a competitive cost compared to its aluminium counterpart [11]. Converting the suspension control arm from aluminium to STBDC only required one major modification; the addition of a third cross member to increase lateral stiffness [11]. A second example is the window frames of the Boeing 787 Dreamliner, shown in Figure 2(b). They used to be made of aluminium but when converted to STBDC, the weight was reduced by nearly 50 % and they gained a superior damage tolerance [4, 12]. (a) Lamborghini Sesto Elemento concept front upper and lower sus- pension control arms. (b) Boeing 787 window frame. Figure 2: Current usage of STBDC materials, [13]. When working with STBDCs it is important to consider that the stochastic nature of the material creates a highly heterogeneous structure containing significant variability of local properties and mesostructure [2, 4]. This results in a material that is locally anisotropic but globally in-plane isotropic [14, 15]. Locally anisotropic refers in this case to that the material properties are direction-dependant on a local scale, even though the average macroscopic response is in-plane isotropic [14]. A drawback of the randomised position and orientation of the tapes is that internal weak spots, e.g. in the form of matrix filled spaces, are distributed throughout the material. These matrix filled spaces occur at two places: i) when tapes lay next to each other without contact, as shown in Figure 1(b), and ii) when one tape drape over another as shown in Figures 7 and 8. Stress concentrations may arise at these matrix rich regions which could cause critical or sub-critical damage to the structure [1]. Due to the irregularities in the meso-scale of STBDCs, no two finished components will ever be identical. This in combination with moulding process sensitivity makes it difficult to characterise general mechanical properties [14]. Thus, increasing the need for predictive analytical as well as numerical models. 1.1 Aim The main aim of this project has been to develop a method to predict elastic properties of STBDCs. This has then been separated in two main goals: i) to develop a method to generate a voxel-based geo- metrical representations of material samples, and ii) to evaluate homogenisation strategies for predicting the volume averaged macroscopic properties. More specifically, for the first goal, the aim has been to develop a modular and multi-functional MATLAB code that based on a set of input parameters outputs a voxel-based numerical model of a flat STBDC plate, along with visualisations of geometrical properties and 3D-models of said plate. For the second goal, the aim has been to investigate, implement, and compare different methods of homogenisation on models of different scale and dimension, namely full 3D models and reduced 2D models of different sizes. 2 1.2 Method The method used to reach the goals of the project can be described in three steps. The first step was a literature study to learn more about the real material and capture the state of the art of modelling and analysing STBDCs. This step gave an important understanding of the micro (fibre), meso (tape), and macro (component) scale properties and effects of the STBDC material, as well as procedures to analyse it numerically. The second step was to use the knowledge gained in the literature study and implement it to create a representative voxel-based numerical model of the material. The final step was to transfer the model to the finite element solver Abaqus, evaluate the in-plane material properties, and compare them to experimental data. The last step was done for both a full-scale 3D-model and a reduced 2D-model. The accuracy of the homogenised material properties was then verified by comparing results to experimental data. The method of homogenisation was validated by comparing results of using homogenised material properties to heterogeneous material properties in a practical example. 1.3 Limitations This project has been conducted as a 30 credit master’s thesis, which amounts to roughly 800 hours divided equally between the two authors. Time was not the only limiting factor but also available com- putational power which limited the size and resolution of the numerical models. The first limitation is that this project has been conducted purely digitally. No prototypes were built, and no experimental tests were made. For experimental reference data, the project relied solely on pub- lished literature. Regarding the numerical model of the material structure. It was set to be a three-dimensional voxel- based structure with cuboid voxels, no tetrahedral or prismatic representations have been used. The macroscopical structure was set to represent a flat plate, to enable the comparison of results to experi- mental test data. The tapes were assumed not to deform in the in-plane direction during the compression moulding phase but only in the out-of-plane direction. They did therefore not get wider nor longer than the pre-determined size. From a top view perspective, the tapes had a constant shape independent of draping or other effects of out-of-plane compression. As a result of this, the cross-sectional area of the tapes remained constant. For the finite element modelling, a few limitations and assumptions were put in place to meet the time deadline. Firstly, the tapes and the matrix were assumed to be perfectly bonded. This meant that adjacent elements with different material properties shared nodes without any additional interface conditions like cohesive contacts. Secondly, only periodic boundary conditions were applied to the model for the homogenisation procedure as they were assumed to give the best response. There was no comparison of periodic boundary conditions to Dirichlet nor Neumann boundary conditions. Thirdly, thermal effects and the resulting residual stresses were not accounted for as they were expected to not affect the prediction of the elastic properties of the material. Fourthly, the model was defined as fully linear elastic, thus without damage or fracture modelling included. Finally, related to the numerical material structure, only cuboid brick elements were used in the 3D finite element analysis, and quadrilateral shell elements for the 2D analysis. This means that the possible implementation of prismatic, tetrahedral, or triangular elements was not explored in this project. 3 2 Previous work The following subsections highlights previous work in the field of stochastic tape-based discontinuous composites. The first part focuses on experimental studies, summarising results and briefly explaining the method used in each study. The results gathered from these studies will later be used to verify the accuracy of the FE models. Then, a number of alternative methods for creating and analysing numerical models with different application capabilities will be summarised to provide an overview of the state of the art. 2.1 Experimental testing In this subsection, three different studies that include experimental testing and manufacturing of STB- DCs will be summarised. Test results and physical dimensions of the plates produced in the different studies are presented in Table 1. Mechanical properties of the tapes used in each study are presented in Table 2. The first study was conducted by Wan and Takahashi in 2016 [14] and was aimed to investigate how the aspect ratio of the tapes affects properties of the final material. Carbon fibre reinforced thermoplastic tapes with a constant width wt = 5 mm and thickness tt = 50 µm were used in the study. The tape length (lt) for the different plates was varied between 12 mm and 30 mm in 6 mm increments. The tape randomisation method used in the study was essentially to distribute tapes in water and use a strainer to pick them up, thus creating thin sheets. The sheets were then stacked to create STBDC plates. The tensile properties were tested following the JIS K 7073 standard but with an increased specimen width to account for the tape length. The plate fibre volume fraction was evaluated using the ash test where the initial mass and volume of a sample was measured, the thermoplastic matrix was then burnt off in a furnace at a temperature of 500 ◦C for three hours. Then, with only carbon fibre strands remaining, the mass was re-measured and the fibre volume fraction calculated. The second study was conducted by Li et al. in 2017 [4], where they manufactured and tested STB- DCs consisting of two different tape thicknesses tt = 164 µm and tt = 285 µm. MATLAB was used to generate a blueprint of the stochastic tape layout. The tapes were then placed by hand according to the blueprint to form the plates. The systematic tape placement makes the position and angle of every single tape known which allows for geometrical cause and effect studies. For example, how a low local fibre content or tape alignment affect the failure mechanisms. Tensile tests were conducted according to the ASTM-D3039 standard but with increased coupon width to account for the tape length. The fibre volume fractions of the plates were evaluated by digital microstructure analysis. They further noted a presence of voids in the thick tape laminate. A third study was conducted as a master’s thesis by Johansen in 2019 [15]. In this study, the randomised layup was achieved by cutting tapes 1.1 m above a substrate and letting them fall freely through a square channel onto the substrate. The substrate was also rotated every 30 seconds (≈ 300 tapes) to avoid any tendency in the landing position. The weight, and later on volume fractions, of the samples were calculated by matrix digestion using nitric acid. The tensile tests specimens were taken from different parts of the plates and cut at different angles which show some variation in the stiffness and strength results. 4 Table 1: Data for STBDC plates manufactured and analysed in [14, 4, 15]. † - values read from plot. Study Tape size [mm] tt [µm] Plate size [mm] E [GPa] σp,fail [MPa] Vf,p [%] Wan [14] 12x5 50 - 42.5 ± 2.5† 440 ± 40† 52.2 Wan [14] 18x5 50 - 49 ± 4† 505 ± 15† 55.1 Wan [14] 24x5 50 - 46 ± 4† 480 ± 60† 52.8 Wan [14] 30x5 50 - 47 ± 2† 520 ± 20† 53.1 Li [4] 50x8 164 290x290x2.444 43.41 ± 4.14 290.1 ± 34.9 57.85 ± 3.84 Li [4] 50x8 285 290x290x2.106 39.66 ± 5.06 186 ± 34.8 53.5 ± 5.39 Johansen [15] 40x20 22.5 275x275x0.5 - - 43.47-49.68 Table 2: Material data for tapes used in [14, 4, 15]. Study Material Et L [GPa] Et T [GPa] σt L [MPa] σt T [MPa] V t f [%] Wan [14] TR50s - PA6 240 - 4900 - - Li [4] HexPly-M77 140 9 2400 73 57 Johansen [15] HS40 - BI018 425 - 4610 - 43.17 It should be noted that in both [4] and [15], the fibre volume fraction of the plate (V p f ) is higher than the original fibre volume fraction of the tape (V t f ). Pre-impregnated tapes were used in both these studies and no extra resin was added to the mould. This means that some resin must have bled through the mould during the curing process. Considering that a part of the resin originally bound to the tapes has bled through the mould, and that a part of the resin has migrated to fill resin pockets, it is very likely that the fibre content of individual tapes in the final plate is higher than when originally placed. 2.2 Digital modelling and finite element analysis Digital models have been created and analysed in terms of elastic properties, failure strength, and failure mechanisms in previous projects. These projects have had different goals and therefore also different methods of implementation. The following paragraphs briefly explain the goals and methodology in five different projects and report the relevant findings. Li et al. [2] wanted to develop a method to analyse large structures while still capturing the local variabil- ity of the STBDC material. The proposed methodology utilises an equivalent laminate based analytical model for prediction of stiffness and strength properties of STBDCs. By varying input parameters, a varying set of material properties are computed. These analytically predicted properties are then assigned to different areas (sets of elements) of the larger structure in order to replicate the stochastic variability of the material. The authors used this method to analyse an automotive monocoque and found that the point of initial failure shifts depending on where weak spots are located. This monocoque example shows that the goal of analysing large structures and capturing local effects of the STBDC material was reached. Selezneva et al. [16] developed a 2D model to evaluate in-plane material properties of STBDCs. The authors used a fine grid with square 2D elements to represent a plate or specimen and randomly placed tapes onto the grid. Multiple tapes were allowed to occupy the same element and classical laminate the- ory was used to evaluate the effective in-plane material properties in each element. The thickness of the plate was pre-defined and constant. In elements where the number of tapes did not amount to the exact plate thickness, the tape thickness was adjusted to compensate, thereby ignoring the presence of resin pockets. The model was used to create and compare specimens with different tape lengths. Abaqus was used for the finite element analysis with pure displacement boundary conditions. The developed model over-predicted the stiffness with roughly 25% for short tape (12 and 25 mm) STBDCs and roughly 6% for long tape (50 mm) STBDCs compared to experimental results. One reason for this over-prediction of stiffness may be the neglection of resin pockets. Another reason why the model over-predicted the modu- lus for short tape STBDCs was that these are more prone to have tapes oriented out-of-plane, which this model neglected. Any out-of-plane tape variation results in a lower in-plane modulus which the 2D model 5 was not capable of capturing. Capturing this out-of-plane effect would require a more complex 3D model. Harper et al. [7] developed an architectural model in which individual fibre bundles were modelled as spring systems. The spring system representation of tapes enabled a more advanced tape interaction model. This model could accurately replicate the compression moulding process and the reaction of individual fibre bundles as they contact other fibre bundles. The fibre bundles were therefore allowed to both bend out-of-plane and twist around their longitudinal axis, with bending and shear stiffness repre- sented by different spring stiffnesses. The authors achieved a reasonable fibre volume fraction (>50%) under the assumption that individual fibre bundles had a fibre volume fraction of 60%. They used the embedded element technique in Abaqus where the fibre bundles were discretised using shell elements (STRI65) embedded in solid brick elements (C3D8) representing the matrix. Pure displacement bound- ary conditions were used in the analysis. The tensile strength and stiffness predictions were both within 5% of experimental results for STBDC fibre volume fractions in the range 40-45%. The difference in stiffness compared to experimental results increased to 30% for fibre volume fractions in the 50-55% range. Another downside to this method was the long model generation time. It took nearly four hours to generate a square plate RVE with side length slightly larger than one tape length and the generation time increased exponentially for larger models. A generation time of 24 hours was reported for a large sized square plate RVE with side length 3.9 times the tape length. Shah et al. [17] wanted to develop a more computationally efficient way of predicting strength and stiffness of STBDC materials. A lot of emphasis was put on developing both a fast model generation procedure and an efficient solution procedure compared to other methods. The tapes were given a prescribed curved shape using spline paths to replicate the out-of-plane distortion, and the tapes were sequentially added to the domain without allowing interactions. This resulted in a tape volume fraction of 61.5% which equates to a fibre volume fraction of roughly 37%. The tapes were then discretised using shell elements embedded in a solid element host mesh representing the matrix material, similar to Harper et al. [7]. Also in this case, pure displacement boundary conditions were used in the finite element analy- sis. The stiffness was underestimated by 12-29% depending on tape length and tape thickness compared to experimental data. Ryatt and Ramulu [1] set out to develop a model to predict elastic and tensile failure properties of STBDC materials. They added voxelised tapes sequentially to a three-dimensional domain and cre- ated digital plates of different thicknesses and widths. I this case, the out-of-plane distortion of tapes was not pre-defined but depended solely on where the tapes were placed. From these digital plates, test coupons of fixed length but varying width were taken out and analysed in Abaqus using reduced integration continuum shell elements (SC8R). They used a twice as coarse mesh for the analysis than the created voxel structure to reduce computational cost and mapped material properties from multiple voxels to a single element. They used pure displacement boundary conditions and found that the com- puted stiffness was slightly higher (4.2%) than the experimental reference value. The authors found no clear trend in modulus as a function of plate thickness nor width. The evaluated strength of the material was found to be dependent on both plate thickness and width but in agreement with experimental results. 6 3 Theory The following section initially covers STBDCs in further detail with regards to manufacturing, mi- crostructure and material properties. Secondly, a thorough introduction on how to geometrically model STBDCs and conduct a FEA is covered in detail in order to ensure a full theoretical background of how the project was conducted, more importantly, why specific choices were made. Finally, a brief overview of the classical laminate theory is covered, introducing key concepts and formulas used in this project. 3.1 Stochastic tape-based discontinuous composites When it comes to carbon fibre composites, the STBDC class of materials is unique. The part manufac- turing process can be described as a mix of the production of short fibre composite parts and continuous fibre composite parts. The resulting microstructure can be structured in layers with subtle hints of the stochastic nature, almost like a continuous fibre composite, but it can also be highly stochastic, mainly depending on the dimensions of the tapes. The effective elastic material properties of the material are typically in-plane isotropic, as long as biases in the randomisation process can be avoided. The stiffness and strength of STBDCs depend not only on the material properties of the constituents but also the mesostructural composition. 3.1.1 Manufacturing process The manufacturing process of STBDCs is visualised in Figure 3, where pre-impregnated carbon fibre tapes are compression moulded to form a geometrically complex part. The process starts with gen- erating the stochastic material structure. Some methods of generating a random layup of tapes were discussed in Subsection 2.1, including tape dropping [15] and a method similar to a wet-type paper making process [14]. Another method is to place small batches of tapes in the mould and shuffle it back and forth between each batch to evenly distribute the tapes [8]. Further, the tapes can either be placed directly in the mould or made into sheets before being placed in the mould [4]. The moulding material can also be bought pre-made in the form of sheet moulding compounds (SMCs) [6]. The second part of the process is the compression moulding, which itself consists of four steps [6]. First, the moulding material is placed in a hot mould, either directly or as a sheet. Second, the mould is closed, the pressure increases, and the material is allowed to flow. How much the material flows depends on the material. However, for STBDCs, it has been observed that the material flow is low, at least for flat plates [4]. Third, the material is cured while the pressure and the temperature is maintained. The curing time depends on the resin, part thickness and temperature [6]. Finally, the mould is opened, the part is extracted and left to cool down. The compression moulding process can be completed in down to two minutes if a fast curing matrix material is used [2]. The moulds may include multiple features such as stiffening ribs, holes, and inserts that reduce post-production machining [6]. a) b) c) Figure 3: The production process of STBDCs, a) pre-impregnated carbon fibre tapes, b) compression moulding phase, and c) complex STBDC part with multiple features [16]. 7 3.1.2 Microstructure Depending on the tapes and the manufacturing process, the microstructure may look very different for different STBDCs [8]. Two examples are shown in Figure 4 where short tapes (12x6 mm) and long tapes (50x6 mm) are compared. The microstructure of the plate using short tapes is more irregular, containing more resin pockets and large out-of-plane oriented tapes. The long tape microstructure, on the other hand, is more resembling a continuous fibre laminate in that the tapes are more structured in layers with less out-of-plane variations. Figure 6. Microstructure of (a) short and (b) long strand ROS specimens. Article Copyright © 2016 Authors, Source DOI: 10.1177/0021998316654748. See content reuse guidelines at: sagepub.com/journals-permissions Figure 4: Microstructure images for STBDC materials with (a) 12x6x0.14 mm tapes, and (b) 50x6x0.14 mm tapes, [16]. STBDC plates can have a varying plate thickness despite a fixed number of stacked tapes in a region, an example of this is shown in Figure 5(a). The varying plate thickness is a result of the amount of resin between the layers and how much the tapes are deformed, which stems from the compression phase [15]. Resin pockets are a common characteristic feature in STBDC materials, appearing either where one tape drape over the end of another tape or where two tapes lie next to each other without contact. Figure 5(b) displays a resin pocket in the microstructure. 4. Results Figure 4.5: A cross-section with varying thickness from 630 µm to 775 µm, but with fixed amount of layers. The white lines inside the market region indicate each layer in that region. Figure 4.6: A cross-section with 18 layers. 20 (a) STBDC microstructure image showing thickness difference in a plate despite a fixed number of layers. White lines indicating different layers. 4. Results Figure 4.7: A cross-section with wavy pattern of layers. The white lines indicate layers. The red line indicates a specifically curved layer with fibre direction parallel to the cross-section. White dots indicate tape edges Figure 4.8: A cross-section with the region around a tape edge marked. The tape edge is at the fibre ends where the tape has been cut. Figure 4.9: A cross-section with the region around a tape edge marked. The tape edge is parallel to the fibres. 21 (b) STBDC microstructure image showing a resin pocket. Figure 5: Microstructure images showing (a) thickness variation, and (b) a resin pocket [15]. 8 Figure 6 shows a comparison of a thin tape (tt = 0.164 mm) and a thick tape (tt = 0.285 mm) STBDC microstructure where multiple characteristic features can be observed. Resin pockets appear in both the thin and the thick tape material, visible in regions 2 and 3, respectively. Further, the varying tape thickness, showcased in Figure 5(a), is also visible in region 1 in Figure 6(a). The final notable characteristic feature is the out-of-plane tape rotation (draping), shown in region 2 in Figure 6(a), where one tape has been pressed over the edge of another tape. Although the microstructure images in Figure 6 only show a small section of the materials, one major difference can be noted between the thin and thick tape STBDCs. The thin tape material shows, in general, a more uniform distribution of tapes with distinct layers and little out-of-plane variations compared to the thick tape material. Region 4 in Figure 6(b) shows this well, where there is no clear boundary between tapes and the layers nor out-of-plane misalignment can be distinguished. (a) Thin (164 µm) tape STBDCs. (b) Thick (285 µm) tape STBDCs. Figure 6: Microstructure images of (a) thin, and (b) thick tape STBDCs [4]. Note that what in this report is called tape is in this image denoted tow. 9 3.1.3 Material properties The macroscale material properties of STBDCs are a result of the micro- and mesostructure, which in turn is a result of the tapes and the manufacturing process. Besides the material properties of the tapes and the polymer matrix, there are multiple factors affecting the macroscopical material properties. These factors are mainly the mesostructure of the material and the dimensions of the tapes, but also e.g. the manufacturing process. One effect of the randomised layup of tapes is that the material is locally anisotropic [2, 4]. Globally, on the other hand, the material is in-plane isotropic as long as there is no bias in the tape distribution [14, 15]. The local variability in the material creates both strong and weak spots where the fibre volume fraction is either high or low, where fibres have a preferred direction, and in the presence of resin pockets [4, 8]. Li et al. [2] analysed the impact of weak spots in their finite element model for an automotive monocoque, discussed in Subsection 2.2. The authors found that the point of initial failure shift depending on where the local weak spots are located, even in the presence of notable geometrical stress concentrations. A consequence of this is that the region at which the material will start to fail also varies, due to the failure of STBDCs typically occurring at material weak spots. A positive aspect of the inherent weak spots is that they aid in making the material notch insensitive [4], which allows holes and inserts to be used in the final product without compromising its structural integrity. Another factor that has been shown to affect the properties of the material are the dimensions of the tapes. The mesostructural dependence on tape length and aspect ratio (lt/tt) is clearly visible in Fig- ure 4. This effect is not only visual, both the stiffness and the strength of the material increase with increasing aspect ratio of the tapes up to a certain point where the effect plateaus [14]. Li et al. [4] experimentally tested STBDCs made of thin tapes with tape thickness tt = 0.164 mm and thick tapes (tt = 0.285 mm) and found that the thick tape plates failed at a lower stress level than the thin tape plates. The earlier failure of the thick tape plates was determined to be partly due to earlier debonding and partly due to a larger microstructural variability than in the thin tape plates. Both the lower aspect ratio and the thickness of the thicker tapes makes them more prone to debonding, i.e. matrix failure between tapes. Also the larger microstructural variability, shown in Figure 6, of the thick tape plates play a role, with a larger presence of out-of-plane oriented tapes and resin pockets. Li et al. [4] found that tape pull-out is more dominant for thick tape STBDCs whereas tape fracture is dominant for thin tapes STBDCs. Tape pull-out and tape (fibre) fracture are the two most common failure modes in STBDC materials [4, 8, 15]. As there are multiple factors affecting the STBDC material properties, it is very difficult to state a typical value for e.g. modulus and stiffness. In examples found in literature, the modulus ranges from 20 GPa to 60 GPa, and the strength ranges from 100 MPa to 500 MPa, depending on the constituents and manufacturing process [4]. 3.2 Geometrical modelling The accuracy of the finite element analysis results, in terms of predicted effective material properties, depend greatly on the accuracy of the geometrical model [17]. Therefore, it is important to capture the major mesostructural effects (random tape distribution, tape draping, fibre volume fraction, etc) accurately. On the other end of the spectrum is computational cost, which increases with the level of detail and size of the model. Therefore, simplifications must be made where it is possible and where the impact on the predicted properties is minimal. The following paragraphs describe some geometrical effects needed to be accounted for and methods of doing so in a numerical model. To create a numerical STBDC model from the ground up, without a physical reference material, there are two general methods to use; the random sequential adsorption technique (RSA) and the Monte-Carlo procedures [18]. In the random sequential adsorption technique, tapes with a random position and angle are sequentially added to a pre-defined domain. When a tape has been accepted in the domain, it is fixed and following tapes are not allowed to intersect previously added tapes. The Monte-Carlo procedures instead add all tapes to a large domain at the beginning. The domain then shrinks to obtain a desired fibre volume fraction by rearranging the tapes. A problem with both these techniques is that without modifications or constraints, the maximum achievable fibre volume fraction is much lower than the real 10 material [18]. For this project, RSA was chosen due to familiarity and relative ease of implementation. Therefore, the Monte-Carlo procedures will not be covered in more detail. Tape interactions Depending on the purpose of the numerical model, different tape interaction conditions may be benefi- cial to consider. In broad terms there are three different tape interaction categories [17], which are all presented in Figure 7. For a reduced 2D model, without bending effects considered, the method shown in Figure 7(a) provides an easy implementation. This method allows tapes to intersect and thereby discards the effect of indi- vidually modelled tapes, presence of resin pockets, and all out-of-plane variations. Analysis wise for this tape interaction method, one option is to use classical laminate theory (CLT) to evaluate effective mate- rial properties where multiple tapes intersect, as done by Selezneva et al. [16], described in Subsection 2.2. A second tape interaction method, shown in Figure 7(b), does not allow intersections between tapes and also discards out-of-plane rotations (draping). This method would result in an inefficiently packed material and a low fibre volume fraction [17]. Therefore, it is not suitable for modelling of STBDCs but may be used as an intermediate step (followed by a Monte-Carlo procedure and a compression phase) or to include bending effects in a reduced 2D model. Finally, Figure 7(c) shows the method with the highest level of fidelity. It prohibits tape intersections and includes out-of-plane draping. Therefore, the tapes may be more efficiently packed, and a higher fibre volume fraction can be achieved. At the same time, the effect of individually modelled tapes is preserved, and the out-of-plane effects are considered. Composites Part B 183 (2020) 107676 2 of fibres [17–20]. Thus, a large number of experiments are required for generating design data, making mechanical property characterization expensive and challenging. This factor has contributed to relatively lesser use of DFC composites. As a response, several researchers pre sented analytical [18,19,21,22] and numerical models [23,24] to pre dict mechanical properties of DFC material. Earlier models like Halpin and Pagano [21,25] were based on CLT. These were originally proposed to model random fiber composites and not just prepreg based DFC’s, and can be used to estimate the in-plane stiffness. The accuracy is limited however, because the size of strands is not considered. Feraboli et al. [19] proposed modified Halpin and Pagano [21] model to predict tensile modulus of DFC material by considering different strand sizes. They, however, considered only the in-plane orientation of strands. More recently, Nakashima et al. [18] proposed an analytical model to predict the flexure modulus of ultra-thin DFC. The proposed model captured the variation in the flexure modulus and accurately predicts the flexure strength. Although the model was effective it was only validated for ultra-thin laminates where through thickness effects are not important. In terms of strength prediction, Selezneva et al. [22] proposed an analytical model to predict strength and failure modes. This model does not account for the out of plane orientation of strands as well as the resin-rich areas. Consequently, it overpredicts the strength of DFC. A more accurate but semi-empirical model was proposed recently by Visweswaraiah et al. [26]. This is based on interlaminar shear and bending stress distribution to predict the strength of hybrid DFC. The authors found that Skin/DFC/Skin configuration gives the highest shear strength (25% higher). More recently, Li and Pimenta [27] have advanced the stochastic equivalent laminate analogy of Feraboli et al. [19] and proposed a multiscale model that uses analytical physically based failure criteria to predict strength and stiffness of tow-based DFCs while accounting for discontinuities in material architecture. In terms of analytical models this is perhaps one of the most comprehensive so far. It still, however, needs a larger data set for fuller validation and does not account for features such as out of plane bending (kinking) of tows (strands) and change of tow thickness and resin rich pockets in different regions of DFC. In addition to analytical models, few authors have also used the FE method to predict the properties of DFC materials. While analytical models are faster and more suitable for producing general trends in property prediction they may be limited by the accuracy of underlying assumptions and requirement of parametrs which are hard to measure physically. It is also not easy to apply these to complex structures unless used in a multi-scale modelling framework. The FE models on the other hand offer the advantage that the methodology once validated can be easily automated and extended to any size of component (given that computational resources are available), even if the user is not an expert of the underlying theory. The accuracy of the FE models however, de pends strongly on the accuracy of the geometric model of internal ar chitecture. A detailed literature review on strength/stiffness prediction of DFC presented by Visweswaraiah et al. [28] identified that the ho mogenized methodologies have been mostly adapted due to their ease in the implementation. However, these approaches overlooked geomet rical effects (waviness) and process parameters (strands dimensions), which significantly influence the characteristics of DFC. The use of FE models that explicitly model the strand geometry thus opens up exciting opportunities for a more realistic predictions by better representing the internal architecture. Generating the strand geometry and its distribu tion accurately and computationally efficiently is challenging, however. Similarly, defining the connection between these hundreds of strands and subsequently generating and solving a 3D continuous mesh is also computationally expensive for any reasonably sized part. Thus, in the next section we have discussed these and other challenges, in the light of recent studies that utilize an explicit representation of internal archi tecture for FEA based prediction of properties of DFC composites. 2. Geometric modelling strategies of DFC material As pointed out in the last section, the primary requirement for a good predictive model for characterization of DFC materials is a reasonably accurate geometric representation of the DFC material architecture. This is particularly important for high fibre volume fraction composites where the architecture contributes as significantly to the material properties as the volume fraction. Generating an accurate geometric representation which is representative of the in-plane as well as out-of- plane distribution is challenging, and simplified models based on con ventional shell approaches often fail at that because these ignore the out- of-plane orientation of strands, defects due to voids, resin-rich areas, adhesion and interlaminar interaction between strands [19,22]. The geometric models of DFC are broadly divided into three cate gories i.e. over-lapped geometric models, non-over-lapped in-plane geometric models, and non-over-lapped out-of-plane geometric models, as shown in Fig. 1. In overlapped geometric models each strand is shared between different cells and this does not represent the continuity of strands. On the other hand, non-over-lapped in-plane geometric models do not allow for representation of the true features of the high-volume fraction DFCs as observed in experiments. Thus, non-over-lapped out- of-plane geometric models have the potential to be most realistic. In this regards the work of Luchoo et al. [29] was a good initial effort. He proposed a non-contact algorithm and spline interpolation to generate 3D architecture of fibre bundles. The main limitations of Luchoo’s work were that it was firstly computationally expensive in terms of solution time. Secondly since 1D truss elements were used the input material behaviour for each tow had to be one-dimensional as well, which is not what is sufficient to model a prepreg strand based DFC where each strand can have a significant width requiring an orthotropic material definition and a non-circular cross-section. Thirdly although the spines of the fibre remained spaced there were possibilities of intersection elsewhere due to the use of circular cross-section profile. The work of LT Harper et al. [30] significantly improved the earlier model and [29] used a force-directed algorithm to generate 3D archi tecture of DFC material. The author’s considered out-of-plane orienta tion, bending and twisting of fibres bundles, and successfully generated RVE of small fibre bundles (3 k having 1.7 mm width) with volume fraction close to 50%. However, for large fibre bundles (12 k having 5.2 mm width) the author found significant intersection/penetration be tween strands with over 40% fibre volume fraction. This model was not specifically validated for the prepreg based DFC, although in theory it seems that one should be able to simply replace tows with pre-preg Fig. 1. Types of DFC geometric models. (a) Over-lapped geometric models, (b) Non over-lapped in-plane geometric models and (c) Non over-lapped out-of- plane geometric models. S.Z.H. Shah et al. Figure 7: Different tape interaction methods: (a) overlapping strands, (b) no overlapping strands, (c) out-of-plane non-overlapped strands (draping), [17]. 11 Draping When one tape is placed over the edge of another tape and then cured under pressure, the upper tape will drape over the lower tape until it reaches an obstacle (another tape or a mould surface). This draping is illustrated in Figure 7(c), and shown in microstructure images in Figures 4, 5(b), and 6(a). Depending on the thickness and bending stiffness of the tapes, the draping ratio varies. For a voxel-based structure, the draping ratio (out-of-plane misalignment angle) is based on the height difference between adjacent voxels [1]. The simplest case, a 1:1 (45°) draping ratio without modifications is shown in Figure 8 a). Gentler slopes, e.g. a 2:1 ratio, may be achieved either by implementing a different draping condition, as shown in Figure 8 b). Or by increasing the in-plane dimensions of the voxels while maintaining the same height, as shown in Figure 8 c). An added benefit of the latter is that it also reduces the number of elements. a) b) c) Figure 8: Voxel-based tape draping, light grey tape draping over dark grey tape, a) 1:1 draping, b) 2:1 draping, c) 2:1 draping with increased element size. Additional constraints for consistent model generation To avoid problems like non-uniform tape distributions and unrealistically low fibre volume fractions, a number of constraints can be used in the model generation [17]. Although the real STBDC material is in fact produced by a random tape placement, there are some effects that a numerical model cannot capture. Two effects that occur in the real material but are often neglected in computational models are tape sliding and tape deformation. In-plane sliding of tapes was observed, although to a limited extent, in [4] and implies that tapes move from their original position during the curing process. Tape deformation means that the cross-sectional area of the tapes does not remain constant during the curing process. Individual tapes get compressed and laterally expand in certain high-pressure areas resulting in thinner but wider tapes. This compression and lateral expansion was observed in [4, 15] and is one of the causes of the thickness variation in Figures 5(a) and 6(a). Both these effects aid in packing the tapes more efficiently and if neglected in a numerical model, the fibre volume fraction will likely not be representative of the real material. It is possible that a completely stochastic model generation procedure can generate a realistic material structure, but not consistently as it relies on probability. To increase the consistency of the model gen- eration, additional constraints can be added to the model generation procedure. To uniformly distribute the tapes over the entire plate, the algorithm created in [17] partitions the plate in a number of bins, like a grid. The algorithm then places tapes in the bins systematically, ensuring that no part of the plate has a significantly higher or lower number of tapes. Without this additional constraint, there is just a statistical probability that the tape distribution will be uniform, no guarantee. To ensure that the fibre volume fraction is representative, there are different constraints that can be implemented. The algorithm developed in [1] has the fibre volume fraction as a target variable and iteratively adds tapes to the domain until the desired fibre volume fraction has been reached. The use of fibre volume fraction as a target value inspired the method used in this project, where the plate is built layer-by-layer, without allowing the algorithm to start a new layer before the first layer has a high enough fibre volume fraction. The implementation of the layer-by-layer method is covered in further detail in Subsection 4.1.2 12 3.3 Finite element analysis The following subsections cover the fundamental theory on how to conduct an FEA with the aim to evaluate homogenised material properties. Initially, the basic concepts that goes into setting up an FEA are detailed, including choice of elements, number of nodes and rule of integration. Secondly, the concept of homogenisation and the general theory of periodic boundary conditions are explained as well as how to evaluate the effective material properties. 3.3.1 Discretisation With a voxel-based structure, the obvious choice is to use cuboid brick (hexahedral) elements. It is also possible to use tetrahedral or prismatic elements, but at a higher risk of bad aspect ratios as the cuboid voxel shape would be cut in at least two segments. In Abaqus, there are a plethora of cuboid elements to choose from, suitable for different problems. For incompressible materials for example, so called hybrid elements (suffix H) are recommended. In this application, regular stress - displacement elements are sufficient but three major choices must be made regarding the discretisation: i) shell or solid elements, ii) number of nodes, and iii) number of integration points and integration method. Solid or shell elements There are two options for hexahedral shell elements in Abaqus: continuum shell elements (SC8R), and continuum solid shell elements (CSS8) [19]. These elements look like solid hexahedral elements but be- have like, and have the same limitations as, regular shell elements. The main limitation of these elements being that they should only be used for thin structures [19]. For the solid hexahedral elements, there is only one choice for the base configuration (C3D• • •) but it may be altered in number of nodes, number of integration points, and method of integration [19]. The solid hexahedral elements are generally not recommended in structures subjected to bending as they are either too stiff or too flexible, due to shear locking and hourglassing, respectively [19]. Number of nodes Regarding the number of nodes for a brick element, Abaqus offers elements with 8, 20, and 27 nodes [19]. An 8-node brick element has one node in each corner, a 20-node element has 12 additional nodes - one on the centre of each edge, and a 27-node element has 6 additional nodes - one on each face and one additional node in the centre of the element. If the minimal number of nodes is used, the edges and surfaces remain straight during deformation, and the displacement approximation is trilinear. For elements with additional nodes, the edges and faces may obtain a curvature as they deform, and the approximation of the displacement is at least quadratic in all directions. Integration rule For the solid hexahedral elements, there are generally three choices regarding the numerical integration for solid elements: full integration, reduced integration, and incompatible modes. Fully integrated el- ements have just enough integration points for numerical integration to achieve exact integration. For an 8-node hexahedral element, this is 2 × 2 × 2 = 8 integration points. Reduced integration has the suffix R in the element description code and means that the number of integration points is reduced. Reduced integration saves computational power and can, in practice, be used to avoid the over-stiff response of fully integrated elements [20]. Typically, the reduced integration cuboid elements give a worse approximation in bending and are not as capable of capturing surface stresses compared to fully integrated elements [19]. Reduced integration elements are also prone to hourglassing that can make them too flexible, which is why hourglass control is included as a default for these elements. Finally, the incompatible modes elements (suffix I) have been implemented to solve some of the problems related to the regular, fully integrated versions of the same element. The C3D8I element for example uses different shape functions than C3D8 and has additional internal degrees of freedom in order to reduce the shear locking phenomena and thereby increase the accuracy in bending. The incompatible mode elements are more computationally expensive than the fully integrated versions. 13 Quadrilateral shell elements For a shell model, the selection of elements is not as wide as for the 3D elements. If the selection is limited to quadrilateral four node elements, excluding triangular elements, it is only one to chose, the S4• • • element with various alterations [19]. In the standard configuration, it is a general purpose shell element with a finite strain formulation, meaning it is suitable for both thick and thin plate structures with large deformations. There is also a reduced integration (suffix R) version with one in-plane integration point instead of four, also including hourglass control [19]. In addition, there are more alterations like the small strain formulations (suffix S) and elements with reduced number of degrees of freedom (suffix 5). 3.3.2 Homogenisation The macroscale mechanical properties of any material depends on the properties of its subscales [21]. In a practical implementation, it is not possible to model every single grain in a steel or every fibre in a carbon fibre composite. Therefore, the material properties are homogenised to effective properties at a higher scale. The rule of mixtures is an example of analytical homogenisation in which the constituent material properties and volume fractions are used to calculate effective material properties. The rule of mixtures can be used in classical laminate theory to calculate some of the effective lamina properties such as the longitudinal modulus [9]. For more complex and aperiodic structures, analytical models may not be accurate enough and computational homogenisation can be used instead. The first step of computational homogenisation is to extract or define a representative volume element (RVE) based on the full microstructure as a sampling space. For a periodic microstructure, the RVE is the smallest volume element that by periodicity represents the material [22]. An example of a periodic microstructure and its RVE is shown in Figure 9(a), where the RVE has the domain Ω and boundaries δΩ+ and δΩ−. For completely irregular (aperiodic) microstructures, the RVE should in theory be in- finitely large [21]. In practice, infinitely large RVEs cannot be used. Instead, the materials are assumed to be statistically homogeneous, meaning that if the RVE is sufficiently large, the position from where it is extracted will not affect the homogenised results. Since an RVE of an aperiodic microstructure is not representative, they are commonly called statistical volume elements (SVEs) [21]. SVEs do not necessarily have be large enough to give the same material properties irrespective of how they are sam- pled, instead multiple smaller SVEs can be homogenised to have the same statistical significance as fewer larger SVEs. An example of an aperiodic microstructure and its SVE is shown in Figure 9(b). The second step of the homogenisation procedure is to carry out the numerical analysis. Different boundary conditions can be used in the homogenisation process, namely: displacement based (Dirich- let), traction based (Neumann), and periodic boundary conditions (PBCs) [23]. In a displacement based homogenisation procedure, the displacement is prescribed to the domain such that it represents a pre- defined macroscopic strain state, and the resulting macroscopic stress state is evaluated as a volume average [21]. The effective macroscopic elasticity parameters are then evaluated from the resulting re- lation between macroscopic strain (prescribed) and macroscopic stress (computed). The evaluation of material parameters is further detailed in Subsection 4.2.3. Depending on if the material response is linear or non-linear, the homogenisation process is different [21]. For a linear material response, the homogenisation can be done once and for all, a priori, and elastic engineering parameters can be determined. For a nonlinear material response, an a priori homogenisation is not sufficient. Instead, one either has to device a two-way coupled system where the homogenisation is nested with the solver [21], or perform off-line so-called virtual testing to generate data to which a macroscale surrogate model can be calibrated. 14 Ω 𝛿Ωା 𝛿Ωା 𝛿Ωି 𝛿Ωି Ω 𝛿Ωା 𝛿Ωା 𝛿Ωି 𝛿Ωି x 𝑦 x 𝑦 (a) Periodic microstructure and its RVE. Ω 𝛿Ωା 𝛿Ωା 𝛿Ωି 𝛿Ωି Ω 𝛿Ωା 𝛿Ωା 𝛿Ωି 𝛿Ωି x 𝑦 x 𝑦 (b) Aperiodic microstructure and its SVEs. Figure 9: Schematic examples of (a) a periodic microstructure and its RVE, and (b) an aperiodic mi- crostructure and its SVEs. 3.3.3 Periodic boundary conditions (PBCs) Periodic boundary conditions (PBCs) are a type of boundary condition that allows for the approximation of a large structure by using a smaller representative volume, typically these would be used on structures with a periodic phase arrangement. However, studies such as [23, 24, 25] indicate that for the case of a heterogeneous material with an aperiodic structure, PBCs provide a better prediction of the effective material properties than both Dirichlet and Neumann boundary conditions, assuming that the SVE is not large enough for all three to converge fully. It has also been shown that the use of Neumann boundary conditions provides the lower bound of the mechanical properties, Dirichlet boundary conditions give the upper bound, and periodic boundary conditions generally fall within these bounds, providing better predictions of the effective mechanical properties for aperiodic structures [24]. A visual representation of this relation between the three is shown in Figure 10. It is important to note that the better prediction of PBCs is the general case, there are special cases where PBCs might produce a less correct estimation of the effective mechanical properties for small domains, as shown in [26]. Similar studies have also come to the conclusion that in terms of convergence, in general, periodic bound- ary conditions are the most efficient compared to the traction (Neumann) and displacement (Dirichlet) boundary conditions [23, 25], this can also be seen in Figure 10 where the PBC reaches the effective value significantly faster. aaaH � " aH11 aH12 sym: aH22 # , �29� the former of which can be recognized as either ÃaH, ÄaH or aH per. Fig. 6 shows the trend of convergence of two-norm of eqn (28). Here, the characteristic length in the ®gure indicates the length scale normalized by that of the model of 32 � 32 pixels (which was used in the above). It can be seen from the ®gure that the two-norms of homogenized elasticity matrices obtained by using periodic boundary conditions are bounded by those by the other two boundary conditions. This result agrees with those by using actually periodic microstructures in Section 2. Furthermore, all the values become close to each other as the RVE size increases irrespective of boundary conditions. In other words, there seems a limit that corresponds to the limit in eqn (21) and therefore we could simulate the converging trend of the e�ective elasticity properties. Since the medium appears to be statistically homogeneous and possibly be isotropic, let us check whether or not the medium has such a feature from the computed homogenized material constants. We ®rst de®ne the parameters such that x � �����ÿ aH1111 aH2222 ����� and X � �����1ÿ aH11 aH22 ����� �30� which would be reasonable measures. For an isotropic nature of the heterogeneous media, the following measures for elasticity matrix and CTE are introduced: Z � max fjaH1112j, jaH2212jg and Y � jaH12j: �31� Figs. 7 and 8 show the convergence trends associated with these parameters by taking the characteristic length normalized by the model of 128 � 128 pixels. As can be seen from these ®gures, the considered medium shows orthotropy or isotropy when the RVE sizes are taken larger and larger. It should also be noted that since the material con®guration is incompatible with the periodicity at the boundary of the Fig. 6. Convergence associated with orthotropy and isotropy for homogenized elasticity matrix. K. Terada et al. / International Journal of Solids and Structures 37 (2000) 2285±2311 2299 Figure 10: Trend of convergence for homogenised elasticity matrix using different boundary conditions, [25]. 15 As a starting point of implementing periodic boundary conditions on an SVE, consider the domain Ω with boundaries δΩ. The external boundaries δΩ are divided into their respective opposing parts δΩ+ and δΩ−, as shown in Figure 9. This boundary division should satisfy the following conditions: (i) δΩ+ ∪ δΩ− = δΩ (ii) δΩ+ ∩ δΩ− = 0 (iii) n− = −n+. (1) The above conditions say that (i) all pairs of opposing boundary sides compose the total boundary δΩ and (ii) the intersection of two opposing sides should be equal to zero, meaning no intersection. The final condition (iii) says that the normal of two opposing boundaries are equal in opposite directions. Therefore, a cuboid domain is suitable in the application of periodic boundary conditions, preferably with a periodic mesh meaning that the mesh distribution is the same on the parallel faces. The displacement field u of an SVE submitted to a macroscopic strain can be divided in to two parts: the displacement due to the applied macroscopic strain denoted ū and a zero-mean fluctuations field ũ, as seen in Equation (2). The periodicity is enforced such that the fluctuating part of two parallel boundaries are equal. This implies that when setting up a displacement constraint between two sides, the only remaining part is the contribution of the linear displacement field (derived from the applied macroscopic strain). This results in the general form of the constraint equation, presented in Equation (3) ui = ūi + ũi = ε̄ijxj + ũi, thus (2) u+ i − u− i = ε̄ijx + j − ε̄ijx − j + ũ+ i − ũ− i︸ ︷︷ ︸ 0 = ε̄ij(x + j − x− j ), (3) where ε̄ij is the macroscopic strain and x+ j and x− j are the respective coordinate of the material point on opposing boundaries. Thus, (x+ j − x− j ) represents the distance between paired external boundaries. Equation (3) ensures the continuity of the displacement field across the domain. A similar constraint is enforced on the traction as: t+ = −t−, defined as an anti-periodicity condition. This comes from the definition of the traction ti = σ̄ijnj which should satisfy condition (iii) in Equation (1) for all pairs δΩ+ and δΩ−. These in turn, also fulfil the Hill-Mandel principle, ensuring that the sub-scale modelling is energetically consistent, and therefore valid. The full theoretical background of the Hill-Mandel principle will not be covered in this thesis, for further reading see [23]. 16 3.3.4 Evaluation of effective material properties The effective material properties are evaluated using the known/assumed constitutive relation (in Voigt notation) between the macroscopic stress (σ̄) and strain (ε̄) via the compliance matrix (S) as ε̄ = Sσ̄. The material at hand is anisotropic on the lower scale, but in the ideal case quasi-isotropic, or in-plane isotropic, on the macroscale [4]. Therefore, in order to evaluate the effective material properties, an orthotropic constitutive relation is used. This is presented in Equation (4), where S is the orthotropic compliance matrix. As orthotropic materials have material properties that differ along three orthogonal axes, quasi-isotropy is a special case of orthotropic materials. In theory, with the material being ideally quasi-isotropic, the two in-plane components of the elastic stiffness (Ex and Ey), the in plane Poisson’s ratios (νxy and νyx), and the out-of-plane shear moduli (Gxz and Gyz) should be identical.  εxx εyy εzz εyz εzx εxy  =  1 Ex −νyx Ey −νzx Ez 0 0 0 −νxy Ex 1 Ey −νzy Ez 0 0 0 −νxz Ex −νyz Ey 1 Ez 0 0 0 0 0 0 1 2Gyz 0 0 0 0 0 0 1 2Gzx 0 0 0 0 0 0 1 2Gxy  ︸ ︷︷ ︸ S  σxx σyy σzz σyz σzx σxy  (4) Extracting the full set of equations from Equation (4) results in the following relations between the different material properties: εxx = σxx Ex − νyxσyy Ey − νzxσzz Ez , εyy = σyy Ey − νxyσxx Ex − νzyσzz Ez , εyy = σyy Ey − νxyσxx Ex − νzyσzz Ez , εyz = σyz 2Gyz , εzx = σzx 2Gzx , εxy = σxy 2Gxy . (5) Using periodic boundary conditions, the applied strain is known for the different load cases. The macro- scopic stress state, however, is unknown but may be evaluated as a volume average using Gauss integra- tion [24] as: σ̄ij = 1 Ω ∫ Ω σijdΩ = 1 Ω Ne∑ e=1 Ne,int∑ I=1 σij(yI)J(yI)W (yI), (6) where σ̄ij is the average stress, e denotes element, yI is the integration point, J is the determinant of the Jacobian, and W is the integration weight. Looping through all elements and integration points of a large RVE is not computationally efficient and the volume average stress may instead be obtained by evaluating the traction at the boundary. Applying the Gauss theorem on Equation (6) yields: σ̄ij = 1 Ω ∫ Ω σijdΩ = 1 Ω ∫ ∂Ω σijnjxjd∂Ω = 1 Ω ∫ ∂Ω tixjd∂Ω ≈ 1 Ω Nnode∑ n=1 Fixj , (7) where ti is the traction, nj the outwards facing normal, xj the position, ∂Ω the boundary, and Fi the nodal force. In practice, the average stress can be evaluated as the total reaction force of a surface divided by the cross-sectional area, as shown in Equation (8), where F tot i is the total force and li is the SVE side length. σ̄xx = F tot x lylz σ̄yy = F tot y lxlz σ̄zz = F tot z lxly σ̄xy = F tot y lylz σ̄xz = F tot z lylz σ̄yz = F tot z lxlz (8) 17 3.4 Classical laminate theory - analytical homogenisation Analysing a laminated composite with layered plies of known angles allows for the use of classical laminate theory (CLT) in order to evaluate the effective in-plane properties: Ē11, Ē22, Ḡ12, and ν12, under the assumption that the out-of-plane stress is zero. CLT does not apply to STBDCs on the macroscale, since the material is neither continuous nor structured in layers. But on a sufficiently small scale, the material is actually both structured in layers and continuous. How CLT was used for homogenisation on a smaller scale is presented further in Subsection 4.3. The smaller scale laminate structure of the STBDC material does not only consist of tapes (plies) but could also consist of pure matrix. Therefore, the ply/lamina/pure matrix will from this point be referred to as layers. Knowing the longitudinal, transverse, and shear material properties for each transverse isotropic layer, the stress strain relation can be expressed with the reduced stiffness matrix Q, defined through its four independent components. The definition of Q in Voigt notation for a transverse isotropic layer and its respective components are stated in Equations (9) and (10). The subscripts L and T in Equation (10) represent the longitudinal and transverse directions, respectively. Furthermore, since the layers can have different angles and the material properties are direction dependent, the reduced stiffness matrix needs to be transformed to comply with the global coordinate system. The transformation between the local and global coordinate system is done through Equation (9), where T 1 is the stress-transformation matrix and T 2 is the strain transformation matrix, the components of these will not be stated in this report but can be found in [9]. Q̄ = T−1 1 QT 2 with Q = Q11 Q12 0 Q12 Q22 0 0 0 Q66  (9) Q11 = EL 1− νLT νTL , Q22 = ET 1− νLT νTL , Q12 = νTLEL 1− νLT νTL , Q66 = GLT (10) Once the transformed stiffness matrix Q̄ has been evaluated for the different possible angles within the system, the laminate stiffness matrix K can be defined as the relation between the moment and stress resultants (M, N), and laminate mid-plane strains and curvatures (ε0,k). The relation has been defined in Equation (11): N · · · M  = [K] ε0 · · · k  = [ A B B D ]ε0 · · · k  , (11) where A, B, and D are the extensional stiffness matrix, coupling stiffness matrix, and bending stiffness matrix, respectively [9]. The evaluation of the A, B, and D matrices has been stated in Equation (12). The factor (hk − hk−1) denotes the distance between the top and bottom surface of each layer relative to the mid-plane, resulting in the layer thickness. A = n∑ k=1 (Q̄)k(hk − hk−1), B = 1 2 n∑ k=1 (Q̄)k(h 2 k − h2 k−1), D = 1 3 n∑ k=1 (Q̄)k(h 3 k − h3 k−1) (12) Once all component matrices of the stiffness matrix K have been evaluated, the effective engineering constants can be calculated using the compliance matrix S [16], defined as: S = [ a b b d ] = [ A B B D ]−1 , where a = a11 a12 a16 a21 a22 a26 a61 a26 a66  . (13) From this, the in-plane material parameters can be evaluated from the inverse components of the com- pliance matrix as: Ē11 = 1 a11h , Ē22 = 1 a22h , Ḡ12 = 1 a66h , ν12 = −a12 a11 , (14) where h is the total thickness of the laminate. 18 4 Methodology In this section, the methodology used to realise the aim throughout the project are detailed. It starts with an in-depth explanation of the geometrical modelling of the STBDC material, including both an overview of the entire procedure and in-depth details on all the major topics. Then covering the finite element analysis, in a more practical implementation manner. This is followed by a description of how the reduced 2D-model was developed and ends with a design of experiment subsection. 4.1 Geometrical modelling The first part of the project was to create a voxel-based numerical model, replicating the microstruc- ture of a compression moulded STBDC. A modified RSA technique was used, partly due to the ease of implementation and partly since it is the most commonly used in literature. This section covers the methodology used in order to create the modular MATLAB code with the functionality of generating a realistic material structure. A flowchart outlining the general functionality of the developed MATLAB code is shown in Figure 11. In short, the flowchart shown in Figure 11 describes the following processes: a main file is used to manage settings such as tape dimensions and toggles such as which plots to plot and if an input file is wanted or not. The geometry is then created using multiple different functions, detailed further in the following subsections. Based on the toggles and settings defined in the main file, the geometrical post-processing is completed and wanted plots produced. The relevant data is then saved in a separate directory for later access. Based on toggles, the sides and top of the domain are cut to a predetermined SVE size to remove edge effects and resin rich top layers. If wanted, an input file is written and saved. Also, a through thickness homogenisation procedure can be done and a reduced 2D-model input file written. Finally, a data check can be run using Abaqus as well as the full simulation. Abaqus Viewer may be automatically opened once the simulation is completed to visualise the results. The entire process can be looped to run multiple simulations in series, without using the Abaqus CAE graphical user interface. Figure 11: Schematic flowchart of the algorithm used to produce and analyse the digital STBDC material structure. 19 4.1.1 Definition of domain For the ease of manipulation and data storage, a 3D array was chosen to virtually depict the full domain of the STBDC plate. For future reference, the indices i, j, k are used to represent the row, column, and layer or y, x, and z axes of the 3D array, respectively. Every tape placed in the domain is stored with an identifiable integer in sequential order (tape index) and the rest of the array is filled with zeros representing the polymer matrix. In terms of size, the length and width of the domain were set to be equal and defined as a function of the tape length. Similarly, the size of the SVEs was based approximately on the tape length. Doing this, a relation between the size of the SVE and full domain could be maintained, ensuring that it was sufficiently large. The height of the domain depends on the maximum number of stacked tapes, which is a user input parameter. The maximum size of the domain was limited by the available memory of the computer. This inhibits the creation of large structures with high resolution, for a computer with 16 Gb RAM memory, the size of the 3D array was limited to roughly 7000×7000×36 elements. For reference, a tape modelled in high resolution could require 4000 elements along the length of the tape. The full domain was partitioned twice into two different subdomains, shown in Figure 12. Firstly, the full domain was partitioned by the safety margin ls to create the safe space subdomain, shown in Figure 12(a). The placement of each tape is governed by its centre position, the centre position may not exceed the boundaries of the internal safe space subdomain. The safety margin was set equal to one tape length, ls = lt. The rationale for this was to avoid any parts of the tape exceeding the outer boundary which would complicate the programming as MATLAB does not allow calling positions outside the array limits. The second partitioning was to divide the safe space subdomain into a number of so-called bins, which were used to force a more uniform tape distribution (further detailed in Subsection 4.1.2). The layout of these bins are shown in Figure 12(b). 𝑙௦ 𝑙௦ Safe space Full domain (a) 3D array with subdomain division. 𝑙௦ 𝑙௦ (b) 3D array with subdomain and bin division. Figure 12: Schematic visualisation of geometric domain and subdomains. 4.1.2 Tape generation and placement The individual tapes making up the domain, with pre-defined width wt, length lt, and thickness tt, were constructed using cuboid volume pixels (voxels). The amount of voxels used to build up each tape was a result of the user-defined tape dimensions, voxel aspect ratio Rvox, and number of voxels in the height direction per tape, N t vox,h (all input-values to the code). The voxel aspect ratio dictates the in-plane dimensions to thickness relation of the voxels. The number of length voxels (N t vox,l), and number of width voxels (N t vox,w) for each tape were calculated as: N t vox,l = round ( ltN t vox,h ttRvox ) , N t vox,w = round ( wtN t vox,h ttRvox ) , (15) where the operator round(•) refers to the MATLAB function with the same name which rounds (•) to the nearest integer value. Increasing the amount of voxels defining each tape results in an increased resolution, providing a more accurate model at the cost of an increased need for computational power. 20 The placement of each tape inside the domain was decided strictly based on its randomised centre posi- tion (Xc, Yc) and tape angle (φ), defined in Figure 14(b). The centre position was randomised, however, the algorithm has a condition implemented to ensure that the search for an open position was more efficient, and to ensure a uniform tape distribution. For each new iteration, the tape content of each bin (see Figure 12(b)) was evaluated and the tape was randomly placed in the bin with the lowest tape content. This ensured not only that the search area was smaller but that the randomised tape placement was more uniform across the different layers. Another restriction set on the modelling was that the placement of tapes were layer constrained (layer condition in Figure 11), this implies that initially tapes may only be placed in the first layer. Once a user-specified desired volume fraction was achieved for the active layer, tapes may be placed one layer higher. This was then repeated until the user-specified number of layers (plies) had been reached. The tape angle was pseudo randomised, meaning it was randomly chosen from a set of angles in the range 0◦-179◦ in equal increments determined by the user. An angle increment of 1◦, for example, yields the following set of angles to choose from: [0◦, 1◦, . . . , 178◦, 179◦]. The reason for using 179◦ instead of 180◦ as an upper limit was to prevent the average tape angle of all tapes to be skew towards the x-direction. If 180◦ was used as an upper limit, both 0◦ and 180◦ would represent the x-direction whereas only 90◦ would represent the y-direction. The tape angle and centre position was sufficient to represent a tape with straight edges but to fit the voxel-based architecture, this linear representation was required to be converted to voxel format, as exemplified in Figure 13. (a) Visualisation of tape discretised with straight edges - linear represen- tation. (b) Visualisation of tapes with ragged edges - voxel format. Figure 13: Schematic visualisation of two different tape representations - linear or voxel. To do the linear to voxel conversion, a Bresenham algorithm created by A. Wetzler was implemented in MATLAB [27]. This is a line drawing algorithm that given two points in a 2-dimensional space provide a pixelated approximation of the line between the two points [28]. Bresenham’s algorithm therefore allows for the transformation of angled lines within the domain into its corresponding pixelated representation, this process is visualised in Figure 14(a). The tapes were created by initially using the corner coordinates in combination with Bresenham’s algorithm to draw all the edges of one tape, as seen in Figure 14(b). The rest of the tape was then filled in by looping through all the x-positions of the tape and filling it from the minimum to the maximum y-position. One drawback to the voxelisation of the tape was that if the resolution was not high enough, multiple given tape angles could result in identical structures. For example, there was no visual difference between tapes with 0◦ and 1◦ tape angles, both were perfectly aligned with the x-axis. 21 (a) Illustration of the functionality of Bresenham’s algorithm. (Xc, Yc) ϕ X Y (b) Visual representation of centre posi- tion and tape angle. Figure 14: Definitions used to define a voxel-based tape Once the tape position within the domain and the voxel-based tape representation was complete, the entire voxelised tape was placed at the lowest available height, without penetrating another tape already placed and not breaching the height restriction. If this could not be achieved at the given position, the centre position and angle were re-selected until the tape could be placed while not penetrating and staying within the active layers. This step was repeated for a user-defined number of iterations. If the maximum number of iterations was reached, the fibre volume fraction target was lowered. This final step was repeated, either until a possible tape position was found or until the fibre volume fraction target was lowered enough to let a new layer be created. 4.1.3 Draping algorithm The process of creating a realistic geometric model emulating a STBDC involved allowing the tapes to adapt as they fall on top of each other (form around each other). In order to capture this during the geometry generation, a draping algorithm was implemented, simulating how the real material acts when cured under pressure. The draping algorithm does not allow for any discontinuity of tapes and the draping ratio (out-of-plane angle) is controlled by the thickness to length ratio of the voxels (Rvox) as shown in Figure 8(c). If the in-plane dimensions of each voxel is increased, the draping angle is decreased. The draping algorithm can be divided into three main parts: i) when a suitable in-plane position has been located for a voxelised tape, all voxels are placed at their lowest available position in the z-direction. ii) Once the tape has been placed, a topology matrix is established, describing the highest position of each voxel position of the current tape from a top-view perspective, as seen in the leftmost image in Figure 15. The purpose of this is to locate any discontinuities of the tape, i.e. where the height difference between adjacent voxels is larger than one. iii) The discontinuities are sequentially removed using a search algo- rithm, looping over all elements in the topology map and calculating the height differences between the current element and its respective neighbours. If any height difference between two neighbouring ele- ments is larger than one, the current element voxel position is moved up one position in the out-of-plane direction (see yellow-marked voxels from left to right in Figure 15). Changing the out-of-plane position of a voxel alters the topology matrix and it has to be updated. Therefore, steps ii) and iii) are repeated until no more changes are made to the out-of-plane position of the tape. A simplified visualisation of how this draping algorithm works is presented in Figure 15. The leftmost illustration shows how a light grey tape lands on a dark grey tape and takes the lowest possible position. The remaining illustrations in Figure 15 show how the nearby voxels move upward to satisfy the one-to- one draping condition. The evolution of the topology matrix for the light grey tape is shown underneath. 22 Figure 15: Draping algorithm example, yellow boxes indicating moved voxels. Corresponding topology matrix shown underneath. The draping ratio was set to 5:1, and achieved by using voxels with a length and width to height ratio (Rvox) of 5:1, as schematically shown in Figure 8 (c). This draping ratio was determined by examining microscope images from three different studies and measuring the height and length of the matrix filled crevices. These measurements are shown in Figure 16. Using a 5:1 voxel aspect ratio did not only achiev