Computational micromechanics of the effect of fibre misalignment on the longitudinal compression and shear properties of UD fibre-reinforced plastics

Computational micromechanics of the effect of fibre misalignment on the longitudinal compression and shear properties of UD fibre-reinforced plastics

Composite Structures 248 (2020) 112487 Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/comp...

2MB Sizes 0 Downloads 52 Views

Composite Structures 248 (2020) 112487

Contents lists available at ScienceDirect

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

Computational micromechanics of the effect of fibre misalignment on the longitudinal compression and shear properties of UD fibre-reinforced plastics T.A. Sebaey a,b,c,⇑, G. Catalanotti d, C.S. Lopes e,f, N. O’Dowd a a

Irish Composites Centre (IComp), Bernal Institute, University of Limerick, Limerick, Ireland Engineering Management Department, College of Engineering, Prince Sultan University, Riyadh, Saudi Arabia Mechanical Design and Production Department, Faculty of Engineering, Zagazig University, P.O. Box 44519, Zagazig, Sharkia, Egypt d Advanced Composites Research Group (ACRG), School of Mechanical and Aerospace Engineering, Queen’s University Belfast, Belfast, UK e IMDEA Materials Institute, C/Eric Kandel 2, 28906 Getafe, Madrid, Spain f Luxembourg Institute of Science and Technology, 5, avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg b c

A R T I C L E

I N F O

Keywords: Micromechanics Fibre misalignment 3D representative volume elements Finite element analysis Longitudinal compression

A B S T R A C T In this paper, 3D representative volume elements (RVEs) of UD carbon/epoxy composites are generated, taking into consideration a realistic misalignment of the fibres. To construct the RVE, a 2D distribution of the fibre sections is considered and then extruded in the longitudinal direction. Experimental measurements to the misalignments are modelled by perturbing the positions of the control points defined on the centreline of each fibre, representing the fibre path as a Bézier curve. The individual fibres are considered as linear elastic orthotropic whereas, the matrix is modelled as isotropic using a damage‐plasticity model. The fibre/matrix interface is modelled using a cohesive formulation law. The IM7/8552 material was simulated using three loading conditions: longitudinal compression, longitudinal shear, and transverse shear. The results show clear correlation between fibre misalignment and compression strength, and stiffness. For shear loading, no effect is recorded. Under the three loading conditions, the predicted material properties and damage propagation are in agreement with the data available in the literature.

1. Introduction Fibre‐reinforced plastics (FRPs) are currently used in a wide range of areas, including aerospace, automotive, marine, and petrochemical. For design of components using composite materials, prototyping is extremely expensive in most of these applications. Consequently, reliable virtual testing (computational tools) is often used to fill this gap. Computational modelling of composites can be performed at different length scales (micro, meso, or macroscale) either separately or in an integrated multiscale approach [1,2]. Typical length scales for describing the response of engineering materials range from atomic/molecular scales (< 109 m) to structural length scales (P 100 m) [3]. Micromechanical modelling predominantly involves numerical modelling of a reduced volume of material [4]. It gives the resulting properties of the material at the mesoscale/macroscale from known characteristics of the constituents, i.e., the matrix and the reinforcement, through the analysis of representative volume element (RVE) [5]. ⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (T.A. Sebaey).

https://doi.org/10.1016/j.compstruct.2020.112487 Received 18 March 2020; Revised 8 May 2020; Accepted 14 May 2020 Available online 30 May 2020 0263-8223/© 2020 Elsevier Ltd. All rights reserved.

Fibre misalignment can be classified, based on the size scale, as either ply misalignment (usually referred to in the literature as fibre/ply waviness or wrinkling) or individual fibre misalignment [6]. Misalignment is one of the characteristics that is mainly affecting the compression properties of the composites. For a multidirectional laminate under general loading conditions, Fedulov et al. [7] showed using mesoscale finite element analysis that the presence of ply misalignment can lead to a reduction of up to 49% in the laminate compressive strength. Using a similar computational approach, Yuan et al. [8] reached a reduction of 40% in the compressive strength by introducing the ply waviness. This effect was confirmed experimentally, by Wilhelmsson et al. [9], who demonstrated that for carbon/ epoxy laminates with a misalignment (on the ply level) of 6°, the compressive strength is reduced by 50%. Davidson et al. [10] studied the effect of a controlled ply waviness on compression strength with misalignment angles ranging from 10° to 35°. The results showed a reduction of 50% of the compressive strength within the range of studied misalignment angles. The effect of introducing artificially induced

T.A. Sebaey et al.

Composite Structures 248 (2020) 112487

represent fibre misalignments in a 3D RVE. In the current paper, the generated RVEs, with the desired fibre misalignment, are numerically analysed using the well‐established material models for the fibres, matrix and interface with the aim of investigating the effect of fibre misalignment on the mechanical response. To the best of the authors’ knowledge, no similar study has appeared in the literature, to date. Numerical simulations are conducted for three different loading case studies using the finite‐element method. In order, will be presented a summary of the method used for geometrical modelling of fibre misalignment, the materials’ constitutive models, and the approach used in setting up the finite element models used in the analysis. Finally, the results will be reported and the effect of fibre misalignment on the investigated mechanical properties will be discussed.

fibre‐wrinkling on the compression strength was investigated experimentally and numerically by Mukhopadhyay et al. [11]. The analysis showed a high stress concentration at the wrinkling area, which initiated early damage and resulted in strength reduction. Larranaga‐ Valsero et al. [12] measured controlled wrinkling in glass‐carbon hybrid laminates using ultrasound techniques. High reductions (up to 63%) in the compression and tensile strength were recorded. All these analyses adopted ply misalignment of kinking. In contrast, there is a lack of studies that introduced misalignment to the individual fibre [13,14], that is where this study fits. The simplest form of micromechanical simulations of FRP is to model a single fibre embedded in matrix, in a plane parallel to the fiber cross‐section, to form a periodic square or hexagonal fibre packing [15]. In most unidirectionally oriented composites, the fibres are distributed randomly in the transverse cross‐section, which is required to provide an accurate RVE [16]. The majority of research on such models considered in‐plane 2D models [17] examining such topics as: methods to represent fibre dispersion in the matrix [18–22], influence of void content (porosity) [23–25]; influence of the resin rich area [26]; material models/properties for matrix [27–31,13], fibre [32] and/or the fibre/matrix interface [33,13]; optimisation the RVE dimensions [34,17,30]; the development of cracks in ply surrounded by two other plies [35–37] and investigation of boundary conditions [38–40]. The outputs of such models are generally homogenized sets of material properties and/or damage identifications. A three‐dimensional RVE was adopted in [17,41] without any consideration of fibre misalignment/waviness. The models assumed that the matrix and fibres were isotropic linear elastic materials. In addition, a perfect bond was implemented between the fibres and the matrix. More details on the matrix behaviour and/or interface characteristics were introduced by Bai et al. [42], Yang et al. [43], and Soni et al. [44] for a 3D RVE with perfectly aligned fibres. Fibre misalignment/waviness was considered in a 3D model by examining a single periodic fibre, with a predefined wavelength and amplitude, embedded in a polymeric matrix in [45–47]. Due to simplifications, this approach is still far from reality and limited on the prediction of the real material behaviour under general loading conditions. Paluch [48] constructed a geometrical model that takes into account the waviness length, amplitude, fibre’s rotation on its own axis and fibre diameter of individual fibre. As per the review conducted by Bishara et al. [49], the problem of this model was that the wavelength was limited to 2 lm to avoid overlap between adjacent fires. In previous publications [14,6], the authors of the current study presented an integrated approach to measure, model and graphically

2. Generation of the RVE with fibre misalignment Deviations from the desired fibre direction occur in FRPs as a result of the manufacturing process and/or the thermal stresses during curing. In an earlier study, [6], the authors presented a methodology based on image processing of 3D computed tomography (CT) scans to assess the misalignment in an FRP composite based on the von Mises distribution. The von Mises distribution uses the (dimensionless) concentration parameter j to characterize the misalignment—the larger the value of j the more aligned the fibres are, and for perfectly aligned fibres (the ideal case) j ¼ 1. For the material systems and manufacturing techniques examined in [6], j ranges from 1000 to 2000. Fig. 1a shows the distribution of the fibre misalignments of a typical carbon fibre/epoxy (IM7/8552) specimen with the best fit von Mises distribution, for j ¼ 1583. Fig. 1b shows the effect of changing the value of j on the von Mises distribution. The process of constructing the 3D RVE is discussed in detail in [6,14]. The modelling part uses the concentration parameter j, the diameter of each fibre df , the fibre volume fraction v f , and the number of fibres nf to construct a 2D RVE using the modules developed by Catalanotti [50] and Varandas et al. [51]. The reader can refer to these two references for detailed information as the current paper does not has any contribution to the 2D construction of the RVE. The centre of each individual fibre in the 2D RVE is extended in the z‐direction and control points are defined on the centres at equal spaces along the RVE depth H. The control points are then perturbed to mimic fibre misalignments using a radius q and angle h from its original position. The value of q is designed to be smaller than df þ D to ensure that

20 = 800 = 1200 = 1600 = 2000 = 2400

15

10

5

0 -0.15

-0.1

-0.05

0

Fig. 1. Characterization of fibre misalignment using von Mises distribution. 2

0.05

0.1

0.15

Composite Structures 248 (2020) 112487

T.A. Sebaey et al.

Fig. 2. Construction of the generated and meshed RVE for the IM7/8552.

Table 1 Constituent materials and interface properties of the IM7/8552 carbon/epoxy composites. Fibre properties [53,54]: Diameter = 5.2 lm; E 11 = 276.0 GPa; E22 = E33 = 15.0 GPa; Matrix properties [13]: E m = 5.07 GPa; mm = 0.35;

m12 = 0.2; G12 = G13 = 15.0 GPa; G23 = 7.0 GPa

t 2 C0 CU rt0 m = 121 MPa; Gm = 90 J/m ; rm = 176 MPa; rm = 180 MPa

Interface properties [55,28,54]:

r0c = 50 MPa; s0c = 70 MPa; gBK = 1.45; lc = 0.52 GIC = 2 J/m2 ; GIIC = 6 J/m2

there is no intersection between the adjacent fibres, with D is the minimum distance between the adjacent fibres. The perturbed control points of each fibre are used to define a Bézier curve for each individual fibre and the angles associated with each segment in the curve is measured. For the whole RVE, these angles are represented by a von Mises distribution and compared to the experimentally measured ones (e.g. as in Fig. 1a). The fibre misalignments for the constructed RVE are optimized to match the  for the constructed RVE with the experimentally measured value of j j. The process continues, by adjusting the perturbations iteratively,  and j are matched within a margin smaller than the maximum until j possible error E M . The control points were then used to construct the whole geometry of the fibres and the matrix of the RVE with holes that match the fibre geometry. The output of the algorithm is shown in Fig. 2.

experiences a typical elastoplastic behaviour. The tensile response is linear and elastic with elastic modulus and Poison’s ratio, E m and mm , respectively, until the tensile failure stress, rt0 m , is reached (see Fig. 3a). Beyond this point, a quasi‐brittle softening is induced in the material, with Gtm being the matrix fracture energy. Under compression, the response is linear up to the initial yield limit, rC0 m . Then, is the matrix strain hardens until the ultimate stress value, rCU m reached. The 8552 matrix properties have been characterised in situ by Naya et al. [13] employing nano indentation techniques and are listed in Table 1. The fibre‐matrix interface is modelled using the cohesive surfaces available in Abaqus [56]. A cohesive interaction between fibre and matrix surfaces is governed by a mixed‐mode traction‐separation law where damage onset is ruled by a quadratic stress criterion that considers the interface strengths in both opening r0c and shear modes s0c [51]. The propagation of damage is governed by linear dissipation of energy with displacement jump for any mode mixity. The corresponding mixed‐mode fracture energy is given by the Benzeggagh‐Kenane law that is governed by a power exponent gBK and the pure mode critical energy release rates (GIC and GIIC ) [57]. Isotropic coulomb friction, with coefficient lc , is coupled with the cohesive behaviour, and enabled at cohesive damage initiation ((see Fig. 3b). Further details on the fibre/matrix interface modelling are presented in [13]. The fibre/matrix interface parameters used in this work were obtained from [55,28,54] and are presented in Table 1.

3. Materials and constitutive behaviour For a micromechanical model of fibre‐reinforced plastic composites, three different regions can be defined as fibres, matrix and fibre/matrix interface. Here, the fibres are considered to be orthotropic, and behave linearly elastically. Fibre damage is neglected as for the present configurations (material and loading cases) it is not expected to be important [52]. The elastic properties of the IM7 carbon fibres have been obtained from [53,54] and are listed in Table 1. To simulate the epoxy matrix the concrete damage‐plasticity model included in Abaqus [56] was used. Details on the matrix constitutive modelling were presented in [13]. This constitutive equation implements a modification to the Drucker‐Prager plasticity yield surface to allow the material to behave as quasi‐brittle when subjected to dominant tensile stress. While under compressive loading, the material

4. Finite element model and boundary conditions To determine the effect of fibre misalignment on the mechanical properties, RVEs are constructed using the method discussed in Section 3. A range of values for the parameter j; 800 6 j  2400, is considered which covers the most relevant materials and manufacturing 3

T.A. Sebaey et al.

Composite Structures 248 (2020) 112487

Fig. 3. a) Schematic of the uniaxial tension-compression response of the epoxy matrix according to the damage-plasticity model for quasi-brittle materials, b) Schematics of the shear response of the damage-friction model for fibre/matrix interfaces. Adapted from [13].

systems, [6]. Most of the studies in the current paper adopt this range of j and compares the results with straight fibres (j ¼ 1). The main features of the RVE are shown in Fig. 4. Based on results from 2D studies, Trias et al. [34] recommended the ratio between the RVE side length L and the fibre diameter df be at least 20 (L=df P 20). Other researches showed that this value is over estimated and significantly affected the computational time, especially for 3D analyses. Gitman et al. [58] recommended L=df ¼ 3 : 7, Okereke et al. [17] recommended L=df ¼ 6, and Pulungan et al. [30] used L=df ¼ 15. To investigate this issue four RVE geometries were examined, with two values of L=df ¼ 8:4 and 14:6 and three values of H=L ¼ 1; 2 and 3, where H is the out of plane dimension of the RVE. The corresponding number of elements for each RVE are listed in Table 2. Considering the random distribution of the fibres in the matrix and to ensure confidence in the results, at least three RVEs are generated and tested for each condition in the sensitivity analysis as well as during the parametric study. The fibres are modelled using C3D8R 8‐node linear brick reduced integration elements with hourglass control. A combination of C3D8R elements and C3D4 4‐node linear tetrahedron elements are used to model the matrix. The element size is optimized through a mesh sensitivity analysis. The sensitivity analysis began with a coarse mesh of 479,163 elements, which corresponds to an element size of 0.75 lm, and continues to refine, down to convergence (convergence of the load–displacement profile regarding the mesh size, i.e., no significant change in the load displacement profile by changing the element size). The model with the finest mesh contained 8,929,419 elements (1,646,326 brick and 7,283,093 tetrahedron elements), which corresponds to an element size of 0.35 lm. A non‐linear geometrical analysis (available in Abaqus) is carried out which takes into account the non‐linear kinematics of the deformation. The analysis used the Abaqus default values of the bulk viscosity parameters as b1 ¼ 0:06 and b2 ¼ 1:2 [56].

Fig. 4. Features and face identification of the developed RVE.

The problem is solved for three different loading conditions to study longitudinal compression, longitudinal shear, and transverse shear. The displacement boundary conditions associated with each case are [39,28]: • Longitudinal F5 uF6 x ¼ uy ¼ 0:0.

Table 2 Dimensions of the RVE used in this study (df ¼ 5:2lm).

compression:

F6 uF5 z ¼ dLC ; uz ¼ dLC

and

F4 F4 F2 • Longitudinal shear: uF2 z ¼ dLS ; uz ¼ dLS , and ux ¼ uy ¼ 0:0.

Size

L (lm)

H (lm)

# of elements

F4 F2 F3 RP • Transverse shear: uF1 y ¼ ux ¼ dTS ; ux ¼ uy ¼ dTS and uz ¼ 0:0.

Size_1 Size_2 Size_3 Size_4

43.845 43.845 43.845 75.942

43.845 87.691 131.548 75.942

479,163 960,293 1,126,898 2,523,780

where ui is the degree of freedom with i ¼ x; y and z; d is the applied displacement at each loading condition, and the superscripts F1‐F6 refer to the faces of the RVE. 4

Composite Structures 248 (2020) 112487

T.A. Sebaey et al.

1600

4000

1400

3000

1200

2000

1000

1000

800

1

2

3

4

5

Fig. 5. Effect of mesh size on the ultimate stresses,

6

1200 900 600 300

0

0

0

0.2

0.4

0.6

0.8

1

1.2

rU , and CPU time T C . Fig. 6. Effect of RVE size on the stress strain behaviour under longitudinal compression.

In finite‐element studies of RVEs, periodic boundary conditions (PBCs) are generally used to represent the condition of a small volume of material (the RVE) surrounded by identical material, representing the overall material volume, as the RVE size is typically a very small fraction of the material being considered. This boundary condition is applied through a set of constraints to correlate the displacement/ strain of opposite faces [39,28,17]. Gitman et al. [58] and Kereke et al. [17] studied the effect of boundary conditions on the predicted elastic properties on composites. Their results showed that the error between using and eliminating periodic boundary conditions ranges from 5 to 30%. In the literature, several authors [59–61] concluded that microstructural RVEs can be solved without PBCs for big enough volume of the RVE. In the current analysis, there are two difficulties associated with the application of the boundary condition. The first is related to the size of the model and the number of constraints resulting from constraining the nodes of the opposite surfaces, which significantly increases the computational time to limits beyond the available computational capacity. The second issue is related to the topology of mesh in the opposite faces i.e., to apply (PBCs) the meshes of the two opposite faces should be identical. In our case, fibres at the boundaries experience misalignment and, it is very difficult to control the number and/or position of nodes in the subsequent tetrahedral meshing [41]. For these reasons, the current analysis does not consider the periodic boundary conditions and apart from the surfaces where the displacement boundary conditions are applied, other boundaries remain traction free and un‐restrained. This assumption may affect the quantitative predictions, but the predicted trends are expected to be representative of the composite behaviour, even without including a periodic boundary condition. For the sensitivity of the model to the mesh size, six different mesh sizes are adopted, with the mesh size previously defined and tested under compression. The values of compression modulus, E 11 , obtained for the meshes are 160.6, 162.8, 163.3, 164.3, 165.0, 165.6 GPa from coarse to fine mesh, respectively, demonstrating that the predicted modulus is insensitive to mesh size. The values of ultimate stresses, rU demonstrating that some sensitivity to mesh size. This sensitivity is further examined in Fig. 5, where rU and computational time T C are plotted as a function of mesh size. It may be noted that the finite‐element simulations are carried out at the Irish Centre for High End Computing (ICHEC) using an eight node machine with 40 central processing unit (CPU) cores per node. The CPU time, in Fig. 5, is calculated as the actual time the model is running  the number of nodes used  the number of CPU cores per node. Moving from the coarsest to the finest mesh increases the computational time exponentially while the value of rU converges at mesh Size_4. For the rest of this article, the same element size will be used as in mesh Size_4. The RVE for this mesh size composed of 2,734,645 elements with a typical run time of 1160 CPU hours.

3000 2500 2000 1500 1000 500 0

0

0.4

0.8

1.2

1.6

Fig. 7. Stress–strain behaviour of the IM7/8552 RVEs at different value of the von Mises concentration parameter j.

3500

172 170

3000

168 2500

166

2000

164 162

1500 1000

160 0

2

4

6

8

10

158 14

12 10

-4

Fig. 8. Effect of fibre misalignment on the longitudinal compression strength and stiffness.

5. Results and discussion 5.1. Influence of RVE size Sensitivity to RVE size is studied for the four RVEs listed in Table 2, at the same value of j ¼ 1600 under longitudinal compression. The same mesh size is used for the four cases. Fig. 6 shows the stress–strain behaviour of the material obtained from the four RVEs considered. The peak stress in each case is controlled by decohesion of the fibre/matrix

5

T.A. Sebaey et al.

Composite Structures 248 (2020) 112487

computational effort, the smallest RVE (Size_1) will be adopted during the subsequent analyses. 5.2. Effect of fibre misalignment on response under longitudinal compression The effect of fibre misalignment on the longitudinal compression properties is examined by constructing RVEs with different values of j. At least three RVEs are constructed and analysed for each j value. Fig. 7 shows a typical stress–strain response for each value of j. Results for the full test matrix are shown in Fig. 8. The results are presented as a function of 1=j for ease of representation (1=j is analogous to the variance in the normal distribution). The strong effect of fibre misalignment on the maximum stress (rU ) is noted in Fig. 8. As j increases (less misalignment), the value of rU increases. This is due to the fact that FRPs, subjected to compressive loads aligned with the fibre direction, fail due to fibre kinking, as a combination of fibre micro‐buckling and the transfer of shear stresses by the interface between the matrix and fibre [42,62–64]. In the current analysis, fibre kinking is predicted by the model, which promotes matrix yielding and fibre matrix de‐bonding, Fig. 9. The same mechanism occurs for all the RVEs in the test matrix under compressive loading. The decrease of rU by increasing the misalignment (reducing j) is explained by the initial micro‐buckling introduced to the individual fibres by increasing the misalignment, that promotes early kinking. It is worth remarking that at the onset of fibre kinking, excessive matrix plasticity and fibre matrix debonding begin to occur simultaneously. Note also from Fig. 9 that fibre kinking occurs over a very small strain range, from 0.93% to 0.96% strain for the case examined in Fig. 9. The composite stiffness (quantified by E 11 ) is also reduced by reducing j, although the effect is much less significant compared to rU . The difference between the highly aligned case (j ¼ 1) and the highly dispersed case (j ¼ 800) is 6.5%. Published values of E 11 of IM7/8552 include 150 GPa [65], 165 GPa [66], and 172 GPa [54]. Our predictions of E 11 (Fig. 8) are in good agreement with the experimentally published values i.e., the current model predicts 1600:75 6 E 11 6 1700:26 GPa, which agrees well with the prescribed experimental published data.

Fig. 9. Fibre kinking phenomenon under longitudinal compression and its effect on the matrix yielding. The blue colour in the figure represents matrix material which has not yielded whereas, the other colours represent matrix material which has yielded. Reader are advised to consult the electronic version for colour identification.

as no matrix softening occurs under compressive loading. In the elastic region, the stiffness E 11 does not show any dependency upon RVE size i.e., the maximum difference between the predicted stiffness is less than 2%, which can be justified by the stochastic nature of the RVE. For the ultimate stress rU , the two RVEs with H=L ¼ 1 (Size_1 and Size_4) show the same value of the maximum stress, despite the significant difference in the volume between the two cases (V 4 ¼ 5:2V 1 ). However, as the ratio H=L increases, the predicted maximum stress also decreases somewhat (compare Size_3, H=L ¼ 3 with Size_1 and Size_4, H=L ¼ 1). However, the difference in the maximum stress between the specimen is less than 5%. Thus, within the scope of this study, the smallest RVE (Size_1) is representative of the material within a margin of 5% of strength and at no deviation in the stiffness, compared with the largest RVE (Size_3). In order to optimize the

5.3. Transverse shear loading Fig. 10 shows typical stress–strain curves for different values of j under transverse shear loading. It may be noted that for misaligned fibre, failure occurs under shear loading at much higher strains than

100 80 60 40 20 0

0

1

2

3

4

5

Fig. 10. Shear stress-strain behaviour of RVEs at different value of 6

j.

Composite Structures 248 (2020) 112487

T.A. Sebaey et al.

Table 3 Predicted values of the strength stiffness under transverse shear loading. Property

s (MPa) U 23

G23 (GPa)

j ¼ 800

j ¼ 1600

3:2

4:3

92:8 4:320:10

92:7 4:340:14

j ¼ 2400

Table 4 Predicted values of the strength and stiffness under longitudinal shear loading.

j¼1

2:6

Property

3:6

91:5 4:320:02

s = s (MPa) U 13

91:8 4:280:05

U 12

G13 = G12 (GPa)

j ¼ 800 2:5

87:0 5:430:15

j ¼ 1600 3:2

89:6 5:430:03

j ¼ 2400 1:3

90:4 5:430:04

j¼1 88:92:6 5:370:14

Damage under longitudinal shear starts by fibre matrix debonding at c ¼ 0:8%. Debonding continues in a plane parallel to the loading direction up to the onset of matrix plasticity and damage, which occurs at c ¼ 2:3%. The combination of the two damage mechanisms is the cause for the drop in the load carrying capacity. This damage profile/mechanism is in agreement with the predictions of Yang et al. [43] and Tan et al. [2] for a similar loading configuration.

those predicted under compression loading. The predicted values of the shear modulus G23 and the ultimate shear stress sU23 are shown in Table 3. It is clear that neither the strength nor the stiffness are affected by the fibre misalignment. The average values obtained from all simulations are, sU23 = 92.2 MPa and G23 = 4.31 GPa. Compared to the published experimental data of the transverse shear modulus (G23 = 3.9 GPa [53,67]), the current predictions are 10% higher, which given the approximations in the model, is an acceptable difference. For the transverse shear strength, the range found in the literature is large (ST ¼ 57 MPa [66]; ST = 85 MPa [68]), which makes any comparison difficult. However, it should be noted that our predicted value falls outside the range of measured values in the literature. The initiation of damage under transverse shear in the model follows a well established propagation mechanism that starts with interface debonding followed by matrix plasticity [43,2]. In Fig. 10, two damage profiles are presented. The red lines in the contours correspond to fiber matrix debonding profile. In both cases damage initiates by debonding and then matrix damage. The final damage profile is not influenced by the value of j (and for the same value of j, different damage profiles may occur). Despite the fact that the principal plane is nearly at 45° in all cases, the crack direction follows the weak inter‐fibre matrix layer.

6. Conclusions The current paper examines the effect of including fibre misalignment in micromechanical models of FRPs. The RVE was constructed to have the same statistical representation (von Mises concentration parameter j) as experimental measurements for this material, [6]. Elastic orthotropic behaviour was assumed for the fibre whereas, for the matrix, the concrete damage plasticity model was adopted. The fibre/matrix interface was simulated using a cohesive law. The effect of fibre misalignments was studied under three loading cases namely longitudinal compression, longitudinal shear and transverse shear. Under longitudinal compression, the fibre misalignment has an effect on both the strength and the stiffness. By reducing j from 1 to 5000, the strength is predicted to reduce by 60%. This drop can be justified by the initial kinking applied to the individual fibres in the form of misalignment. Following this sudden reduction, the strength is almost linearly proportion to j.For the stiffness, the reduction of the strength by reducing j is also observed and can be justified by the area of the fibre in the loading direction. However, the sensitivity of the stiffness (measured by E 11 ) to j is significantly less than for the strength, less than 7% by reducing j from 1 to 800. The analysis did not recognize any effect of the misalignment on the failure mode i.e. all the tested RVEs failed by the onset of fibre kinking, that promotes matrix damage. For the shear loading (either longitudinal or shear) the analysis did not identify any effect of the misalignment on the strength, stiffness

5.4. Longitudinal shear Fig. 11 shows the predicted response under longitudinal shear. A summary of the maximum stress (s13 ) and the associated shear modulus (G13 ) is introduced in Table 4. Again, under longitudinal shear the predicted sensitivity to j is weak. The average values for sU12 and G12 are 89.0 MPa and 5.42 GPa, respectively. The published data for the ultimate shear strength SL and longitudinal shear modulus G12 are 89.9 MPa [69] and 5.29 GPa [53], respectively. The comparison shows a good agreement with the published data with deviation of 1% and 2.5% for both the strength and stiffness, respectively.

100 80 60 40 20 0

0

1

2

3

4

5

Fig. 11. Shear stress-strain behaviour of RVEs at different value of j under longitudinal shear loading. The overall damage is shown to the right whereas, zoom view of a small portion is shown in the box to see the separation between the fibres and matrix and the matrix plasticity. 7

T.A. Sebaey et al.

Composite Structures 248 (2020) 112487 [12] Larranaga-Valsero B, Smith RA, Tayong RB, Fernandez-Lopez A, Guemes A. Wrinkle measurement in glass-carbon hybrid laminates comparing ultrasonic techniques: a case study. Compos A 2018;114:225–40. [13] Naya F, González C, Lopes CS, Van der Veen S, Pons F. Computational micromechanics of the transverse and shear behavior of unidirectional fiber reinforced polymers including environmental effects. Compos A 2017;92:146–57. [14] Catalanotti G, Sebaey TA. An algorithm for the generation of three-dimensional statistically representative volume elements of unidirectional fibre-reinforced plastics: focusing on the fibres waviness. Compos Struct 2019;227(111272). [15] Kathavate KS, Pawar DN, Bagal NS, Adkine AS, Salunkhe VG. Micromechanics based models for effective evaluation of elastic properties of reinforced polymer matrix composites. Mater Today: Proc 2020;21:1298–302. [16] Schmid CF, Switzer LH, Klingenberg DJ. Simulations of fiber flocculation: effects of fiber properties and interfiber friction. J Rheol 2000;44:781–809. [17] Okereke MI, Akpoyomare AI. A virtual framework for prediction of full-field elastic response of unidirectional composites. Comput Mater Sci 2013;70:82–99. [18] Seyedalikhani S, Shokrieh MM, Shamaei-Kashani AR. A novel dynamic constitutive micromechanical model to predict the strain rate dependent mechanical behavior of glass/epoxy laminated composites. Polym Testing 2020;82(106292). [19] Melro AR, Camanho PP, Pinho ST. Generation of random distribution of fibres in long-fibre reinforced composites. Compos Sci Technol 2008;68:2092–102. [20] Yang L, Yan Y, Ran Z, Liu Y. A new method for generating random fibre distributions for fibre reinforced composites. Compos Sci Technol 2013;76:14–20. [21] Vaughan TJ, McCarthy CT. A combined experimental-numerical approach for generating statistically equivalent fibre distributions for high strength laminated composite materials. Compos Sci Technol 2010;70:291–7. [22] Vaughan TJ, McCarthy CT. Micromechanical modelling of the transverse damage behaviour in fibre reinforced composites. Compos Sci Technol 2011;71:388–96. [23] Chowdhury K, Talreja R, Benzerga AA. Effects of manufacturing-induced voids on local failure in polymer-based composites. J Eng Mater Technol 2008;130:021010. [24] Ashouri Vajari D. A micromechanical study of porous composites under longitudinal shear and transverse normal loading. Compos Struct 2015;125:266–76. [25] Mehdikhani M, Petrov NA, Straumit I, Melro AR, Lomov SV, Gorbatikh L. The effect of voids on matrix cracking in composite laminates as revealed by combined computations at the micro- and meso-scales. Compos A 2019;117:180–92. [26] Ghayoor H, Marsden CC, Hoa SV, Melro AR. Numerical analysis of resin-rich areas and their effects on failure initiation of composites. Compos A 2019;117:125–33. [27] Melro A, Camanho P, Pires F, Pinho S. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part I-Constitutive modelling. Int J Solids Struct 2013;50:1897–905. [28] Melro AR, Camanho PP, Andrade Pires FM, Pinho ST. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: Part II- Micromechanical analyses. Int J Solids Struct 2013;50:1906–15. [29] Romanowicz M. A numerical approach for predicting the failure locus of fibre reinforced composites under combined transverse compression and axial tension. Comput Mater Sci 2012;51:7–12. [30] Pulungan D, Lubineau G, Yudhanto A, Yaldiz R, Schijve W. Identifying design parameters controlling damage behaviors of continuous fiber-reinforced thermoplastic composites using micromechanics as a virtual testing tool. Int J Solids Struct 2017;117:177–90. [31] Canal LP, Llorca J, González C, Segurado J. Intraply fracture of fiber-reinforced composites: microscopic mechanisms and modeling. Compos Sci Technol 2012;72:1223–32. [32] Swolfs Y, Morton H, Scott AE, Gorbatikh L, Reed PAS, Sinclair I, Spearing SM, Verpoest I. Synchrotron radiation computed tomography for experimental validation of a tensile strength model for unidirectional fibre-reinforced composites. Compos A 2015;77:106–13. [33] González C, Llorca J. Mechanical behaviour of unidirectional fibre-reinforced polymers under transverse compression: microscopic mechanisms and modelling. Compos Sci Technol 2007;67:2795–806. [34] Trias D, Costa J, Turon A, Hurtado JE. Determination of the critical size of a statistical representative volume element (SRVE) for carbon reinforced polymers. Acta Mater 2006;54:3471–84. [35] Madadi H, Farrokhabadi A. Development a refined numerical model for evaluating the matrix cracking and induced delamination formation in cross-ply composite laminates. Compos Struct 2018;200:12–24. [36] Herráez M, Mora D, Naya F, Lopes CS, González C, Llorca J. Transverse cracking of cross-ply laminates: a computational micromechanics perspective. Compos Sci Technol 2015;110:196–204. [37] Fu C, Wang X. Micro-mechanical analysis of matrix crack-induced delamination in cross-ply laminates in tension. Compos Struct 2020;243:112202. [38] Coenen E, Kouznetsova V, Geers M. Novel boundary conditions for strain localisation analyses in microstructural volume elements. Int J Numer Methods Eng 2012;90:1–21. [39] Barbero EJ, Cosso FA, Campo FA. Benchmark solution for degradation of elastic properties due to transverse matrix cracking in laminated composites. Compos Struct 2013;98:242–52. [40] Sádaba S, Herráez M, Naya F, González C, Llorca J, Lopes C. Special-purpose elements to impose periodic boundary conditions for multiscale computational homogenization of composite materials with the explicit finite element method. Compos Struct 2019;208:434–41. [41] Bahmani A, Li G, Willett TL, Montesano J. Three-dimensional microscopic assessment of randomly distributed representative volume elements for high fiber volume fraction unidirectional composites. Compos Struct 2018;192:153–64.

nor damage initiation and propagation. Under transverse shear, the damage is governed by fibre–matrix debonding whereas, under longitudinal shear, the damage initiates by fibre–matrix debonding and is governed by a combination of the debonding and matrix failure. On the other hand, the predicted mechanical properties (in the three loading cases) are in good agreement (less than 10% induced error) with the data available in the literature for the material under investigation. Data Availability The raw/processed data used herein to justify these findings can be shared upon request. Interested researchers can directly contact the corresponding author. CRediT authorship contribution statement T.A. Sebaey: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Validation, Visualization, Writing ‐ original draft, Writing ‐ review & editing. G. Catalanotti: Conceptualization, Methodology, Software, Writing ‐ review & editing. C.S. Lopes: Software, Supervision, Validation, Writing ‐ review & editing. N. O’Dowd: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Writing ‐ original draft, Writing ‐ review & editing. Acknowledgement This work has received funding from Enterprise Ireland (EI) and the European Union Horizon 2020 Research and Innovative Programme under the Marie Sklodowska‐Curie grant agreement No. 713654. The first author would like to acknowledge the support of Prince Sultan University under the S&M LAB research laboratory grant. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, athttps://doi.org/10.1016/j.compstruct.2020.112487. References [1] Raju B, Hiremath SR, Mahapatra DR. A review of micromechanics based models for effective elastic properties of reinforced polymer matrix composites. Compos Struct 2018;204:607–19. [2] Tan W, Naya F, Yang L, Chang T, Falzon BG, Zhan L, Molina-Aldareguía JM, González C, Llorca J. The role of interfacial properties on the intralaminar and interlaminar damage behaviour of unidirectional composite laminates: Experimental characterization and multiscale modelling. Compos B 2018;138:206–21. [3] Okereke MI, Akpoyomare AI, Bingley MS. Virtual testing of advanced composites, cellular materials and biomaterials: a review. Compos B 2014;60:637–62. [4] Nijssen RPL, Brøndsted P. Advances in Wind Turbine Blade Design and Materials. Cambridge, UK: Woodhead Publishing Limited; 2013. [5] Bonora N, Ruggiero A. Micromechanical modeling of composites with mechanical interface - Part 1: Unit cell model development and manufacturing process effects. Compos Sci Technol 2006;66:314–22. [6] Sebaey TA, Catalanotti G, O’Dowd N. A microscale integrated approach to measure and model fibre misalignment in fibre-reinforced composites. Compos Sci Technol 2019;183(107793). [7] Fedulov BN, Antonov FK, Safonov AA, Ushakov AE, Lomov SV. Influence of fibre misalignment and voids on composite laminate strength. J Compos Mater 2015;49:2887–96. [8] Yuan Y, Yao X, Niu K, Liu B, Wuyun Q. Compressive failure of fiber reinforced polymer composites by imperfection. Compos A 2019;118:106–16. [9] Wilhelmsson D, Gutkin R, Edgren F, Asp LE. An experimental study of fibre waviness and its effects on compressive properties of unidirectional ncf composites. Compos A 2018;107:665–74. [10] Davidson P, Waas AM, Yerramalli CS, Chandraseker K, Faidi W. Effect of fiber waviness on the compressive strength of unidirectional carbon fiber composites. 53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 2012. [11] Mukhopadhy S, Jones MI, Hallett SR. Compressive failure of laminates containing an embedded wrinkle; experimental and numerical study. Compos A 2015;73:132–42.

8

Composite Structures 248 (2020) 112487

T.A. Sebaey et al.

[55] Varna J, Berglund L, Ericson M. Transverse single-fibre test for interfacial debondlng in composites: 2. Modelllng. Compos A 1997;28:317–26. [56] SIMULIA User Assistance. ABAQUS User’s Manual, 2017. Dassault Systemes Simulia Corp.; 2017.. [57] Benzeggagh ML, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass. Compos Sci Technol 1996;56:439–49. [58] Gitman IM, Askes H, Sluys LJ. Representative volume: existence and size determination. Eng Fract Mech 2007;74:2518–34. [59] Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 2003;40:3647–79. [60] Kanit T, N’Guyen S, Forest F, Jeulin D, Reed M, Singleton S. Apparent and effective physical properties of heterogeneous materials: representativity of samples of two materials from food industry. Comput Methods Appl Mech Eng 2006;195:3960–82. [61] Chen G, Ozden UA, Bezold A, Broeckmann C. A statistics based numerical investigation on the prediction of elasto-plastic behavior of WC-Co hard metal. Computat Mater Sci 2013;80:96–103. [62] Budiansky B, Fleck NA. Compressive failure of fibre composites. J Mech Phys Solids 1993;41:183–211. [63] Camanho PP, Arteiro A, Melro AR, Catalanotti G, Vogler G. Three-dimensional invariant-based failure criteria for fibre-reinforced composites. Int J Solids Struct 2015;55:92–107. [64] Zhou L, Zhao L, Liu F, Zhang J. A micromechanical model for longitudinal compressive failure in unidirectional fiber reinforced composite. Results Phys 2018;10:841–8. [65] Hexcel. Hexply 8552 epoxy matrix. Product data sheet, Hexcel Corporation; 2016.. [66] Chamis CC, Abdi F, Garg M, Minnetyan L, Baid H, Huang D, Housner J, Talagani F. Micromechanics-based progressive failure analysis prediction for WWFE-III composite coupon test cases. J Compos Mater 2013;47:1–18. [67] Gan KW, Hallett SR, Wisnom MR. Measurement and modelling of interlaminar shear strength enhancement under moderate through-thickness compression. Compos A 2013;49:18–25. [68] Vogler M. New material modeling approaches for thermoplastics, composites and organic sheets. In: 9th European LS-DYNA Conference. [69] Catalanotti G, Camanho PP, Marques AT. Three-dimensional failure criteria for fiber-reinforced laminates. Compos Struct 2013;95:63–79.

[42] Bai X, Bessa MA, Melro AR, Camanho PP, Guo L, Liu WK. High-fidelity micro-scale modeling of the thermo-visco-plastic behavior of carbon fiber polymer matrix composites. Compos Struct 2015;134:132–41. [43] Yang L, Wu Z, Cao Y, Yan Y. Micromechanical modelling and simulation of unidirectional fibre-reinforced composite under shear loading. J Reinf Plast Compos 2015;34:72–83. [44] Soni G, Singh R, Mitra M, Yan W. Modeling multiple damage mechanisms via a multi-fiber multi-layer representative volume element (M2RVE). Sadhana 2020;45:64. [45] Garnich MR, Karami G. Finite element micromechanics for stiffness and strength of wavy fiber composites. J Compos Mater 2004;38:273–92. [46] Herráez M, Bergan AC, González C, Lopes CS. Modeling fiber kinking at the microscale and mesoscale. Technical Report NASA/TP-2018-220105, NASA, NASA, Langley, Research Center, Hampton; October 2018.. [47] Naya F, Herráez M, Lopes CS, González C, Van der Veen S, Pons F. Computational micromechanics of fiber kinking in unidirectional frp under different environmental conditions. Compos Sci Technol 2017;144:26–35. [48] Paluch B. Analysis of geometric imperfections affecting the fibers in unidirectional composites. J Compos Mater 1996;30:454–85. [49] Bishara M, Rolfes R, Allix O. Revealing complex aspects of compressive failure of polymer composites - Part I: Fiber kinking at microscale. Compos Struct 2017;169:105–15. [50] Catalanotti G. On the generation of RVE-based models of composites reinforced with long fibres or spherical particles. Compos Struct 2016;138:84–95. [51] Varandas LF, Arteiro A, Bessa MA, Melro AR, Catalanotti G. The effect of throughthickness compressive stress on mode II interlaminar crack propagation: a computational micromechanics approach. Compos Struct 2017;182:326–34. [52] Naghdinasab M, Farrokhabadi A, Madadi H. A numerical method to evaluate the material properties degradation in composite rves due to fiber-matrix debonding and induced matrix cracking. Finite Elem Anal Des 2018;146:84–95. [53] Arteiro A, Catalanotti G, Melro AR, Linde P, Camanho PP. Micro-mechanical analysis of the in situ effect in polymer composite laminates. Compos Struct 2014;116:827–40. [54] Arteiro A, Catalanotti G, Melro AR, Linde P, Camanho PP. Micro-mechanical analysis of the effect of ply thickness on the transverse compressive strength of polymer composites. Compos A 2015;79:127–37.

9