A combined use of embedded and cohesive elements to model damage development in fibrous composites

Composite Structures 223 (2019) 110921

Contents lists available at ScienceDirect

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

A combined use of embedded and cohesive elements to model damage development in fibrous composites

T

⁎

Qiang Liu , Larissa Gorbatikh, Stepan V. Lomov Department of Materials Engineering, KU Leuven, Kasteelpark Arenberg 44 bus 2450, 3001 Leuven, Belgium

A R T I C LE I N FO

A B S T R A C T

Keywords: Embedded element technique Cohesive element Interfacial damage Composites

In recent years the embedded element (EE) technique, which is based on the superposition of meshes, has been introduced into finite element (FE) modelling of composite materials. When compared to conventional approaches employing a single continuous mesh this technique has a number of advantages, among which are simplicity in the model creation and improved mesh quality. However, due to the absence of physical interfaces, simulation of interfacial damage with this technique presents a real challenge. In this work we propose an approach in which the EE technique is combined with the cohesive zone model with the purpose to model interfacial damage in fibrous composites. The new hybrid method is validated against conventional FE models in four case studies. Among others the methodology is applied to predict mechanical properties of nanocomposites reinforced by carbon nanotubes showing good agreement with experimental data.

1. Introduction Analysis of the mechanical behaviour of composite materials using finite element (FE) models of their representative volume elements (RVE) is an essential part of material’s design. In conventional FE models that rely on continuous meshes, creation of the RVE model involves in division of the volume into different domains, such as reinforcement and matrix, and this division is performed by partition or subtraction. For brevity we will call these models ‘direct’. The use of the direct models is generally limited by three challenges: (1) the volume partition or subtraction is difficult to perform for materials including a large number of reinforcing elements with high aspect ratios, for example, short fibers, carbon nanotubes (CNTs), graphene flakes and so on; (2) the complex geometry of the RVE after the volume partition or subtraction is detrimental to the quality of FE mesh and to (3) the number of elements/degrees of freedom, which is significantly increased. These drawbacks make direct FE predictions computationally expensive and often lower their accuracy. In recent years the embedded element (EE) technique, employing superposition of meshes, has been introduced to overcome the highlighted drawbacks of direct FE modelling. In an EE-based model the reinforcements and matrix are created and meshed separately, thus allowing for well-structured meshes in all constituents. The reinforcements and matrix are generally defined as embedded and host domains, respectively. The nodal displacements of the embedded domain are

⁎

interpolated from those of the host domain, so that the two domains are coupled together. Up to now the EE-type models have been used for analysis of various composite materials and structures (such as short fiber reinforced [1–5], laminated [6–8], textile based [9–13], CNT reinforced [14–22], concrete [23–26] and biological [27] composites), which cannot be efficiently analysed and in some cases even analysed at all by the direct approach. Good precision of the EE technique for prediction of elastic stress fields and stiffness of composites has been noted in several works and assessed using comparisons against direct modelling results and experimental data [1,4,9,13]. Despite the progress, there is one critical drawback of existing EE models – they cannot model interfacial damage/debonding. This is because the nodal displacements of the reinforcements are directly correlated with those of the matrix and the relative separation and sliding between the two constituents is not allowed. Due to this fact EEbased models are mostly limited to predictions of the elastic stress fields and stiffness of composites, with possible inclusion of non-linear behaviour of the matrix. In several works the EE technique was also invoked to predict failure of fibrous composites, however, only damage in the matrix was considered (via continuum damage mechanics), and interfacial damage not [8,12,19,24]. Very limited attempts have been done so far to introduce interfacial damage into the framework of the EE models. One of the common characteristics in these attempts is that the reinforcements are modelled using one-dimensional (1D) elements. To account for interfacial

Corresponding author. E-mail address: [email protected] (Q. Liu).

https://doi.org/10.1016/j.compstruct.2019.110921 Received 15 January 2019; Received in revised form 16 April 2019; Accepted 24 April 2019 Available online 25 April 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 223 (2019) 110921

Q. Liu, et al.

damage, two types of approaches have been developed. In the first approach, the debonded reinforcements are replaced by new effective reinforcements which are perfectly bonded to the matrix [6,18]. The equivalent reinforcements are modelled using 1D elements, and their behaviour can be calibrated by the direct FE approach [18]. In the second approach [3,7], the interface is discretised as a set of 1D elements, which are embedded into matrix to realize the load transfer from the matrix to the reinforcement. Modelling the interface using 1D elements is efficient from the point of view of computational costs, but this simplification is not always reliable as the mixed-mode damage on the interface cannot be captured by 1D elements. The objective of the current work is to advance the EE technique towards modelling of interfacial damage directly. To achieve this goal we propose to combine the EE technique with the cohesive zone model (CZM) and investigate internal workings of this coupling. Introduction of cohesive zones is a common method to simulate interfacial damage [28,29]. We will call this new hybrid method EE-CZM. To the best of our knowledge, this is the first study to directly model interfacial damage in fibrous composites by combining the EE technique and cohesive elements. The two techniques have been previously combined to model matrix cracking [8] or delamination between plies [12], but the cohesive elements served there as ‘a host domain’ together with the matrix. The most important difference between the current method and those in the Refs. [8,12] is that the current method allows for modelling of interface damage in the framework of superimposed meshes where a physical object such as interface does not even exist. The interface is an elusive entity that lies between the embedded domain (reinforcements) and the host domain (matrix). The accuracy of the new method is evaluated through comprehensive comparisons with direct FE models in several benchmark examples. Its mesh-sensitivity is also investigated. The proposed method is further validated by predicting the mechanical behaviour of realistic composites. Without loss of generality the present study is performed on composites with fibre-like reinforcements, such as short microscopic fibers, CNTs, etc. Therefore the embedded domain will be often called “fibre”, and the host domain “matrix”.

easy to mesh using structured meshes (e.g. C3D8 in ABAQUS) due to the ideal geometry of the two phases. Matrix domain is created separately, and no explicit partition or subtraction is required. Hence the matrix domain can also be meshed using regular meshes. In order to model load transfer from matrix to fiber by the interface, the outer surface of the interface layer is defined as an embedded domain as shown in Fig. 1(b). This means that only nodes on the outer surface of cohesive elements (e.g. nodes 1 and 2 in Fig. 1(b)) are coupled with the matrix nodes by the EE technique, while nodes on the inner surface (e.g. nodes 3 and 4) are excluded from the embedded domain. The translational displacements of the outer-surface nodes are calculated from those of the adjacent matrix nodes, and the inner-surface nodes have the same displacements as those on the fiber outer surface. There is an important difference between the present approach and other EE-based models and it lies in the definition of the embedded domain. While in the literature this domain is defined as a volume, in the present approach it is a surface, more specifically, the outer-surface nodes of cohesive elements. Combining embedded and cohesive elements this way has not been previously done. If the whole volume of the interface layer is defined as an embedded domain, the interfacial displacements are correlated with those of the matrix and are not independent degrees of freedom. Compared to the two-phase EE models, the present model requires an extra step of introducing the interface layer (i.e. partition of the fiber), but this does not create difficulties in model creation and meshing. Similarly to other types of EE-based models, the EE-CZM model introduces extra matrix material in the region overlapped by the fiber, which would increase the overall composite stiffness compared to the direct model. To eliminate this unphysical effect, the stiffness of the embedded fiber elements should be discounted [2,4], i.e.

E ∗f = Ef − Em,

Ef and Em are Young’s moduli of the embedded elements, where fiber and matrix, respectively. With the use of Eq. (2), the reinforcement effects of real fiber are modelled by the superposition elements: the fiber elements with discounted stiffness and the extra matrix elements overlapped by the discounted fiber. In other words, there are two different ways to model the load transfer from the matrix to the reinforcement: 1) the matrix – cohesive elements – discounted fiber and 2) the matrix to the matrix region overlapped by the discounted fiber. Therefore the cohesive stress in the EE-CZM model only represents the load transferring from the matrix to the discounted fiber, and it is smaller than the cohesive stress in the direct model. To ensure an accurate prediction of the interfacial damage initiation, the strength parameters assigned to cohesive elements should be adjusted according to

2. Methodology 2.1. Combination of the EE technique with CZM For an EE model without interface, fiber is generally considered as an embedded domain, and matrix a host domain. The two domains are coupled through a series of constraint equations, which can be expressed as 8

uiI =

∑ Sk (ξ k, ψk , ηk ) uik , (i = x , y, z ), k=1

(2)

E ∗f ,

(1)

E Si∗ = Si0 ⎛⎜1 − m ⎟⎞, i = n, s, t , Ef ⎠ ⎝

uiI

where represent translational displacements of a node, I, in the embedded domain, uik – the adjacent nodes, k, in the host domain, and Sk (ξ k , ψk , ηk ) are interpolation basis functions in 3D space. With the use of Eq. (1), the stiffness properties of the two constituents are combined. The two-phase (i.e. fiber and matrix) model does not take interface effects into account. To include them, a three-phase model (e.g. Fig. 1(a)), where an interface layer is inserted between fiber and matrix, is generally employed in direct FE modelling. The interface layer is modelled using cohesive elements and has its own material properties. One side of the cohesive elements shares common nodes with the fiber elements, and the other side with the matrix elements. The thickness of the cohesive elements can be set to zero. A novel three-phase model is proposed here, which combines embedded and cohesive elements in a single model and is to be called EECZM model. The basic idea of the new method is demonstrated in Fig. 1(b). An additional layer is added around the fibre (to be called “outer layer”) to represent the interface. Both fiber and interface are

(3)

where Si∗ and Si0 are strength values of the cohesive element and realistic interface, respectively, while n, s, t represents normal, longitudinal shear and transverse shear directions of the cohesive element, respectively. The compensation for the stiffness of fiber elements (Eq. (2)) and interfacial strength of cohesive elements (Eq. (3)) can be neglected when fiber volume fraction is relatively low or when stiffness of the fiber is much larger than that of the matrix. 2.2. Numerical examples for validation Fig. 2 shows four numerical examples used to examine workings and accuracy of the proposed EE-CZM approach. To facilitate fundamental understanding, all numerical examples are based on single fiber composites. The fiber is either modelled as straight or wavy, both of which are common in real composites. Without loss of generality, dimensions 2

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 1. Three-phase FE models for fibrous composites with fiber/matrix interface taken into account: (a) conventional model with a single continuous mesh and (b) hybrid model combining embedded and cohesive elements. Orange colour is used for the fiber, light green for the matrix and semi-opaque for the interface.

Fig. 2. Different numerical examples used for validation of the EE-CZM approach: (a) and (b) for analysis of load transfer mechanism from a single straight and wavy fiber to the matrix, respectively; (c) and (d) for analysis of debonding and pull-out behaviour of a single straight and wavy fiber from the matrix, respectively. Orange colour marks the fiber, light green the matrix and semi-opaque the interface. 3

Composite Structures 223 (2019) 110921

Q. Liu, et al.

and material parameters of the fiber are assumed to be of multi-walled CNTs [15]. Its radius is R = 4.5 nm, the length when it is straight is L = 0.6 μm, and geometrical parameters of the wavy fiber are introduced in Fig. 2(d). The numerical examples in Fig. 2(a) and (b) are used to study the load transfer mechanism from matrix to fiber in the presence of interface but without damage. Numerical examples in Fig. 2(c) and (d) are used to investigate interfacial damage during the pull-out process of a single fiber. For each case study, both the FE models with continuous and embedded meshes are analysed. The mesh configurations in the refined regions are depicted in Fig. 2(a) and (b) as insets. According to the Cox theory, fiber ends have only minor contributions to the load transfer, therefore they are assumed to be detached from the matrix in all direct models. In the direct models, volume partition leads to extremely irregular matrix domain, and hence small-sized tetrahedral (i.e. C3D10) elements are required. In contrast the EE-CZM model can be meshed with high-quality structured elements. Without further specifications, the mesh size is set as 0.3R × 1.0R for the matrix, and 0.2R × 0.5R for the fiber in the four numerical examples.

3.1. Load transfer mechanism This subsection explores validity of the EE-CZM approach in predicting load transfer from the matrix to fiber before interfacial damage onset. Hence very small strain is applied to ensure interfacial stress is at a low magnitude. 3.1.1. Straight fiber A tensile strain of 0.05% is applied on the models in Fig. 2(a) along X or Z directions, corresponding to the longitudinal or transverse tension of a single fiber composite. Fig. 3 shows modelling results obtained with the direct and EE-CZM approaches when the single fiber composite is under longitudinal tension. The axial stress along the fiber is plotted in Fig. 3(a), it reaches its maximum at the centre and gradually decreases to zero towards fiber ends. The results obtained in the absence of fiber/matrix interface are also provided for comparison. As the two-phase model is equivalent to a three phase model with an interface of infinite stiffness, the latter predicts a lower stress in the fiber compared to the former. However, the relative error in the maximum stress is only 1.6%, which can be neglected. The interfacial shear stress along the fiber is plotted in Fig. 3(b), it gradually increases from the centre to fiber ends. The stress distribution in the matrix is shown in Fig. 3(c), significant stress concentrations are observed in the matrix around fiber ends, as expected by shear-lag theory. It can be calculated that the average stress in the matrix domain not perturbed by the fiber is 1 MPa; therefore, the stress concentration factor induced by fiber ends for the given mesh size is about 7.0 for the longitudinal tension. As shown in Fig. 3, the EE-CZM approach reproduces the stress profiles and fields predicted by the direct FE model although there are some minor differences. Table 2 lists the elastic strain energy of the fiber, matrix and interface predicted by the direct and EE-CZM approaches. The relative errors between the two methods are less than 2% for the fiber and matrix domains and less than 4.2% for the interface. The computational time of the direct and EE-CZM approaches is about 2.42 and 2.51 min, respectively. Similar comparisons are made for the case of transverse tension (i.e. loading along Y direction). As shown in Fig. 4, both interfacial normal and shear deformations are activated when the single fiber composite is under transverse tension. The profiles of interfacial normal and transverse shear stress (τs) around the fiber are shown in Fig. 4(a) and (b), respectively. Due to the Poisson effect, the material contracts in the longitudinal direction leading to a build-up of the interfacial longitudinal shear deformation. The profiles of interfacial longitudinal shear stress (τt) and axial stress along the fiber are shown in Fig. 4(c) and (d), respectively. For the four types of stress profiles, the EE-CZM predictions show good accordance with those obtained from direct models. The stress fields in the matrix and fiber, corresponding to the middle cross section of the fiber, are shown in Fig. 5 for the single fiber composite subjected to transverse tension. For the stress field obtained from EE-CZM approach, fiber cross section is indicated by a white circle in Fig. 5(b), where a low stress region is observed. It implies that the matrix domain overlapped by the fiber has a negligible contribution to the total strain energy, therefore, this domain can be equivalent as a region without matrix material, i.e. the case in direct model in Fig. 5(a). Some small differences in the stress fields of the matrix (i.e. Fig. 5(a) and (b)) and fiber cross section (i.e. Fig. 5(c) and (d)) are observed, which is partially due to the different meshes in the direct and EE-CZM models (Fig. 2(a)). For example, the stress level in the fiber obtained

2.3. Material parameters For all the computational tests, both fiber and matrix are assumed to be elastic. The material properties of the fiber are those of CNTs [15], and those of the matrix correspond to a typical epoxy resin [30]. The fiber/matrix interfacial behaviour is assumed to follow a bilinear traction-separation law. Before damage the response is governed by the following linear relationship with stiffness, K, (4)

τi = Ki δi, i = n, s, t ,

where τi and δi are interfacial stress and displacement in the three directions, respectively. The damage initiation is governed by a quadratic-stress criterion, which is expressed as: 2

2

2

⎛⎜ 〈τn 〉 ⎟⎞ + ⎜⎛ τs ⎟⎞ + ⎜⎛ τt ⎟⎞ = 1, 0 0 0 ⎝ Ss ⎠ ⎝ St ⎠ ⎝ Sn ⎠

(5)

where 〈 〉 is the Macaulay brackets with 〈x〉 = max(0, x), and Ss0 and St0 are corresponding interfacial strengths. The full failure of fiber/ matrix interface (i.e. complete debonding) is controlled by the respective fracture energies for loading by the normal and shear tractions, GnC , GsC and GtC . The interfacial properties in normal and shear directions are assumed to be the same (namely, Sn0 = St0 , GnC = GtC ), and are set to the experimental data measured by the pull-out of a single CNT from epoxy [31,32]. The interfacial strength and fracture energy are taken as 36 MPa and 0.9 J/m2, respectively. The key material parameters are summarized in Table 1. These parameters are only used for numerical validation of the EE-CZM approach, that is, the results in Section 3. The parameters used in Section 4 are calibrated from the experiments on the benchmarked CNT nanocomposites, which will be introduced later.

Sn0 ,

3. Analysis and discussions The modelling results corresponding to the four numerical examples in Fig. 2 are discussed in this section. Quantitative comparisons between direct and EE-CZM approaches are performed systematically to examine accuracy of the latter. To allow for the comparison of computational time between the two methods, all simulations are carried out using 8 cores (intel Xeon 6140) and 32 GB RAM. Table 1 Key material parameters of fiber (‘f’), matrix (‘m’) and interface. Ef ,1(GPa)

Ef ,2 (GPa)

ν12

Gf ,12 (GPa)

Gf ,23 (GPa)

Em (GPa)

νm

Kn = Ks (GPa/m)

Sn0 = Ss0 (MPa)

GnC = GsC (J/m2)

475

10.3

0.3

27.9

3.8

2.0

0.4

109

36

0.9

4

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 3. Comparisons of the direct and EE-CZM approaches in the analysis of load transfer in a single fiber composite under longitudinal tension: the profiles of (a) axial stress in the fiber and (b) interfacial shear stress and (c) stress distribution in the matrix.

simulation of transverse tension, the direct and EE-CZM approaches take about 2.43 and 2.44 min, respectively.

from the EE-CZM approach is lower than that from the direct model. When a single fiber composite is under transverse tension, a much smaller stress concentration factor in the matrix (about 1.5) is generated compared to the case of the longitudinal tension. This is due to the anisotropic geometry and elastic properties (i.e. the stiffness in the transverse direction is much smaller than that in the longitudinal direction) of the fiber. Table 2 reports the elastic strain energy stored in the fiber, matrix and interface for this case. The EE-CZM approach and the direct model show good agreements. The EE-CZM method predicts a higher strain energy for the fiber than the direct approach, although the former predicts a lower stress in the fiber cross section than the latter as shown in Fig. 5(c) and (d). The reason for this is that the longitudinal contraction due to the Poisson effect has a major contribution to the strain energy of fiber. This conclusion can further be verified by comparing Fig. 4(d) with Fig. 5(c) and (d); the axial compressive stress on the fiber is much larger than the cross-sectional stress. It indicates that the load transfer from the matrix to the fiber is dominated by interfacial longitudinal shearing, even for the transverse tension, while interfacial normal and transverse shear deformations play secondary roles. Therefore, the EE-CZM method slightly overestimates the strain energy in the fiber for both longitudinal and transverse loadings. For the

3.1.2. Wavy fiber When a composite with a single wavy fiber is mechanically loaded, the fiber and the interface are always under mixed-mode deformations. Thus only one type of external load is considered here, namely the tensile displacement along X direction, Fig. 2(b). The axial stress in the fiber, obtained from direct and EE-CZM modelling, is plotted in Fig. 6(a). Fig. 6(b) and (c) show the interfacial normal and longitudinal shear stresses plotted along two representative paths, denoted as paths I and II in Fig. 6(d). The interfacial transverse shear stress along the two paths is of much lower magnitude and hence not given here. As for the three stress profiles reported in Fig. 6(a)–(c), good agreements between direct and EE-CZM modelling results are obtained. Stress distribution in the matrix is shown in Fig. 6(e). Unlike in the case of the straight fiber, the matrix around the wavy fiber is under a complex stress state, and hence the maximum principal stress is shown. The direct and EE-CZM approaches predict nearly identical stress fields, both showing significant stress concentrations around fiber ends. It is found that the stress concentration factor in the matrix predicted by the

Table 2 Comparisons of the elastic strain energy of each constituent in a single (straight) fiber composite obtained by the direct and EE-CZM approaches (unit = 1 × 10−18 J). Under longitudinal deformation

Direct EE-CZM Error (%)

Under transverse deformation

Fiber

Matrix

Interface

Fiber

Matrix

Interface

0.702 0.716 2.08

26.373 26.379 2.31 × 10−2

2.91 × 10−2 3.03 × 10−2 4.20

0.116 0.118 1.52

26.001 26.005 1.42 × 10−2

2.68 × 10−3 2.63 × 10−3 1.90

5

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 4. Comparisons of the direct and EE-CZM approaches on the problem of a single (straight) fiber composite subjected to transverse tension: the profiles of (a) interfacial normal and (b) interfacial transverse shear stress, (c) interfacial longitudinal shear stress and (d) axial stress in the fiber.

Fig. 5. Comparisons of the direct and EE-CZM approaches in analysis of the stress distribution in the matrix ((a), (b)) and fiber cross section ((c), (d)) for the single (straight) fiber composite subjected to transverse tension. 6

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 6. Comparisons of the direct and EE-CZM approaches in analysis of load transfer in the single (wavy) fiber composite under 0.05% tensile deformation in X direction: (a) the axial stress in the fiber, (b) interfacial normal stress and (c) interfacial longitudinal shear stress along the fiber, (d) a schematic illustration of paths I and II locating at fiber surface and (e) stress distribution in the matrix.

The elastic strain energy stored in each constituent in the single (wavy) fiber composite is reported in Table 3. For the strain energies of the fiber and matrix, the EE-CZM approach produces negligible deviations from the direct model. A higher error (i.e. 8.65%) is found for the prediction of interfacial strain energy, but the accuracy of the EE-CZM approach is still acceptable. The computational time taken by the direct and EE-CZM approaches is about 3.30 and 3.60 min, respectively.

Table 3 Comparison of the elastic strain energy of each constituent in the single (wavy) fiber composite obtained by the direct and EE-CZM approaches (unit = 1 × 10−18 J).

Direct EE-CZM Error (%)

Fiber

Matrix

Interface

0.527 0.536 1.65

112.850 112.857 6.20 × 10−3

2.89 × 10−2 3.14 × 10−2 8.65

3.2. Pull-out behaviour of a single fiber direct model is about 31.47. The high value is caused by sharp corners in the matrix, resulting from volume partition during model creation. This problem is commonly present in all direct models of fibrous composites. The concentration is much lower in the EE-CZM approach.

In this subsection, interfacial damage is introduced and the pull-out behaviour of a single straight (Fig. 2(c)) and wavy (Fig. 2(d)) fiber from the matrix is studied. The pull-out displacement applied at the fiber end is 0.5 μm. 7

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 7. The pull-out behaviour of a single straight fiber: (a) force-displacement curve, (b) the profiles of interfacial shear stress corresponding to moments (I)–(IV), stress fields in the matrix at moments of (c) I and (d) IV.

3.2.1. Straight fiber Fig. 7 shows the pull-out behaviour of a single straight fiber from the matrix obtained by the direct and EE-CZM approaches. The pull-out force as a function of the applied displacement is shown in Fig. 7(a), where a good agreement is observed between the two methods. The relative error in the peak force is 0.06%. The interfacial longitudinal shear stress along the fiber, corresponding to four moments indicated in Fig. 7(a), are shown in Fig. 7(b) to demonstrate the development of interfacial damage. There is no interfacial damage at moment I as interfacial shear stress of all cohesive elements is less than the strength (i.e. 36 MPa). At moment II, damage starts and partially extends along the interface, especially in the region close to the fiber end on which the displacement is applied, but the region between points E and F is still free of damage. At moment III, the damage-free region is greatly reduced to that between points G and H, and the pull-out force reaches its peak value. At moment IV damage is observed along the entire interface, and interfacial shear stress decreases significantly due to the gradual accumulation of damage. Stress distribution in the matrix is shown in Fig. 7(c) and (d) corresponding to the respective moments I and IV, where the fiber positions are indicated. The matrix domain close to the fiber end loaded by the pull-out displacement experiences higher stress than other domains at moment I. With the increase of pull-out displacement, the matrix domain around the whole fiber exhibits high level of stress at moment IV. The direct model, including some sharp geometric features, predicts a higher value for the maximum stress in the matrix than the EE-CZM approach, similarly to the results in Section 3.1. As seen in Fig. 7, the pull-out process of a single straight fiber described by the EE-CZM method consistent with the predictions of the

direct approach in terms of the pull-out force–displacement curves, interfacial stress profiles and stress fields in the matrix. For this studied case, the direct and EE-CZM approaches take 1.46 and 1.38 h, respectively. 3.2.2. Wavy fiber The pull-out behaviour of single wavy fiber predicted by the two methods is shown in Fig. 8. There is a good agreement between the direct and EE-CZM approaches in predicting the force-displacement curves as shown in Fig. 8(a). The relative difference in the peak force is 1.8%. Interfacial shear (τt) and normal (τn) stresses along the two paths depicted in Fig. 6(d) are shown in Fig. 8(b)–(d) for moments I-III, respectively. In Fig. 8(b)–(d), the upper pictures show the stress profiles in the longitudinal shear direction and the lower ones for those in normal direction. Since the compressive stress in the normal direction does not lead to damage, the region with τn smaller than −20 MPa is not shown in Fig. 8(b)–(d) for clarity. The critical locations at which interfacial damage starts to occur are depicted in Fig. 8(b) and (c). During the pull-out process the interface of the wavy fiber undergoes partial damage that transitions into full debonding. In contrast to the case of a straight fiber, interfacial normal stress is herein comparable to shear stress, and hence interfacial damage is governed by a mixed mode. As shown in Fig. 8(b)–(d), the EE-CZM approach succeeds in reproducing the development of interfacial damage predicted by the direct model, although the error between the two methods is slightly higher than that in the case of the straight fiber. Stress fields in the matrix are shown in Fig. 8(e) and (f) corresponding to moments I and III, respectively. The EE-CZM approach 8

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 8. The pull-out behaviour of a single wavy fiber: (a) the force-displacement curve, (b)–(d) the respective profile of interfacial shear stress corresponding to moments (I)–(III), stress fields in the matrix at moments of (e) I and (f) III.

predicts almost the same patterns of the stress fields as the direct model for both moments. A significant stress concentration is observed for the domain where the fiber has a small radius of curvature, but the volume of this high-stress domain is a very small portion compared to the whole matrix. The maximum stresses predicted by the EE-CZM method are 23.6% and 10.2% lower than those from the direct approach at moments I and III, respectively. This difference may be due to deficiencies from both approaches, for example by the unavoidable use of tetrahedral meshes in the direct model. Mesh quality cannot be well guaranteed in the matrix domain including small-sized voids left by volume subtraction, which may lead to numerical singularity and decrease the accuracy of the direct approach. For EE-CZM approach, the stress

concentrations in the matrix induced by stiff inclusions can be modelled, but the ones caused by soft inclusions (i.e. the hole left by a full fiber debonding) cannot be captured. For the pull-out modelling of wavy fiber, the direct and EE-CZM approaches spend 7.60 and 6.63 h, respectively. 3.3. Mesh sensitivity The validity of the EET-CZM approach would depend on the size of fiber (or interface) and matrix meshes. Here the mesh sizes of fiber cross-section and matrix are denoted as Lm and Lf, respectively, and they are expressed as a ratio to the fiber radius R. The mesh size along 9

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 9. Mesh sensitivity in the EE-CZM analysis of the load transfer in a single (wavy) fiber composite: the fiber strain energy for fiber mesh sizes of (a) 0.5R and (b) 1.0R and different matrix mesh sizes, respective interfacial energy in (c) and (d), the illustration of some (e) fiber, interface and (f) matrix meshes.

fiber length is 2 × Lf. Without loss of generality, only the numerical examples of single (wavy) fiber composites, i.e. Fig. 2(b) and (d), are used for investigating the mesh-sensitivity. The mesh size in the numerical example in Fig. 2(b) is changed to study its effects on the load transfer from the matrix to fiber. The applied strain is increased to 8% so that interfacial damage can be activated. The strain energies, consumed by the fiber and the interface, are extracted to quantify the effects of mesh size on the load transfer. The

strain energies in function of external strain are shown in Fig. 9(a)–(d). Three different sizes are considered for fiber meshes, i.e. 0.2R, 0.5R and 1.0R, and five for matrix meshes, i.e. 0.3R, 2.0R, 8.0R, 12R and 20R. The fiber and matrix meshes with some sizes are depicted in Fig. 9(e) and (f), respectively. As the EE-CZM calculations with Lm = 0.3R and Lf = 0.2R match well with those of the direct model, they are considered as references in this mesh-sensitivity study. For the modelling results in Fig. 9(a) and (c), fiber mesh size is fixed as 0.5R 10

Composite Structures 223 (2019) 110921

Q. Liu, et al.

are still of good accuracy in the stage prior to the peak force as shown in Fig. 10(b), but show higher deviation for the post-peak response. However, the model with Lf = 1.0R and Lm = 12.0R still has an acceptable accuracy relative to the reference case: relative errors in the peak force, corresponding displacement and pull-out energy are 2.2%, 4.2% and 5.1%, respectively.

and matrix mesh sizes change from 2.0R to 20R except for the reference case. Fiber mesh size is 1.0R for the results in Fig. 9(b) and (d), and other sets are kept the same. Fig. 9(a) and (b) show the fiber strain energy changing with applied strain. The strain energy firstly increases with strain but subsequently decreases due to the onset of interfacial damage, no matter the size of the fiber and matrix meshes. Interfacial energy monotonically increases with strain in the studied range as shown in Fig. 9(c) and (d). This is because interfacial damage gradually develops from fiber ends to its centre: when more cohesive elements engage into the damage process the total interfacial energy increases. As indicated by Fig. 9(a)–(d), coarser fiber and matrix meshes both result in higher strain energy of the fiber and the interface, which is due to the fact that coarser meshes overestimate structural stiffness. When Lf is fixed to 0.5R, the size of matrix mesh in the range from 0.3R to 12R has very limited influences on fiber strain energy (Fig. 9(a)) and interfacial energy (Fig. 9(c)). Compared to the reference case (i.e. Lm = 0.3R and Lf = 0.2R), the model with Lf = 0.5R and Lm = 12.0R produces only 3.9% and 0.6% higher strain energies on fiber and interface, respectively. When Lf equals to 1.0R, in comparison with the reference case, a noticeable difference is observed in fiber strain energy at all matrix mesh sizes (Fig. 9(b)). On the other hand, the maximum strain energy of fiber, predicted by Lf = 1.0R and Lm = 12.0R, is still only 7.4% higher than the reference value, and the error is only 1.7% for interfacial energy. Consequently, EE-CZM with the mesh size of Lf = 1.0R and Lm = 12.0R is still acceptable, especially when the composite properties are dominated by the interface behaviour. Fig. 10 shows the mesh sensitivity of the EE-CZM approach for predicting the pull-out behaviour of a single wavy fiber from matrix. Similarly, fiber mesh size is fixed to 0.5R and 1.0R in Fig. 10(a) and (b), respectively, and matrix mesh size varies from 2.0R to 20.0R except for the reference case. For both Lf = 0.5R and 1.0R, coarser matrix meshes lead to a higher peak force but a lower critical displacement corresponding to it. For Lf = 0.5R, the pull-out force-displacement curves exhibit negligible differences as matrix mesh sizes vary, providing that Lm does not exceed 20.0R. Compared to the reference case, the model with Lf = 0.5R and Lm = 12.0R predicts 1.1% lower value for the peak force, 4.2% lower critical displacement, and 0.5% higher pull-out energy. The model with even coarser meshes, i.e. Lf = 0.5R and Lm = 20.0R, provides a better prediction of the peak force (0.7% error relative to the reference case), a comparable precision for the pull-out energy (0.6% error), but much worse prediction for the critical displacement corresponding to the peak force (12.5% error) as shown in Fig. 10(a). When Lf is fixed to 1.0R, predictions with the coarse matrix meshes

4. Modelling of realistic CNT nanocomposites In this section, the proposed EE-CZM approach is employed to predict the mechanical behaviour of realistic CNT nanocomposites, which can further examine the capability of this novel modelling method. The experimental data for CNT reinforced epoxy in [30] are considered as the benchmark for validation of the EE-CZM approach. A 3D RVE model of a nanocomposite is constructed as shown in Fig. 11. According to experimental data [30], CNTs’ radius and length are 10 nm and 2.5 µm, respectively. The dependence of predicted mechanical properties on the RVE model size was studied in [33,34], with a conclusion that the model size larger than the average length of CNTs is large enough to obtain representative simulation results. Thus, in this study dimensions of the RVE are chosen as 3.0 × 3.5 × 3.0 µm3 in the x, y, z directions. The CNTs are uniformly distributed and randomly orientated as shown in Fig. 11(a). The reader is referred to [15] for a description of how to model wavy CNTs. With the use of EE-CZM, the matrix is meshed by well-structured C3D8R elements in Abaqus as shown in Fig. 11(b). Interface layers are inserted around CNTs and meshed by 8-node cohesive elements (i.e. COH3D8) like Fig. 1(b). The geometric thickness of cohesive elements is one percent of the CNT radius. Based on the mesh-sensitivity study in Section 3.3, the matrix mesh size is selected as 10 times of the CNT radius, and the cross-sectional and longitudinal sizes of CNT meshes are 1 and 4 times of the CNT radius, respectively. The morphologies of CNT and interface elements can be referred to Fig. 9(e) and (f) in Section 3.3. The RVE model is extended in Y direction as shown in Fig. 11 to predict the tensile behaviour of the CNT nanocomposite up to failure. According to the experimental data in [30], the Young’s modulus and Poisson’s ratio of epoxy matrix are 2.0 GPa and 0.4, respectively. The ductile damage behaviour of the matrix is modelled using the material model proposed in [35]. Damage in the matrix is assumed to start when the maximum principal strain exceeds 2.1%, which is determined by calibrating the experimental data of epoxy in [30]. CNTs are assumed to be free of damage and modelled as linear elastic material. They are represented as continuum solids, their multi-wall structure is not modelled explicitly. As the average outer and inner radii of CNTs were reported as 10 and 4 nm in [30], respectively, the number of walls is taken as 9 and each wall has a thickness of 0.34 nm. Assuming the

Fig. 10. Mesh sensitivity in the EE-CZM analysis of the single (wavy) fiber pull-out behaviour: the fiber mesh size of (a) 0.5R and (b) 1.0R and different sizes of the matrix meshes. 11

Composite Structures 223 (2019) 110921

Q. Liu, et al.

Fig. 11. RVE model of a CNT nanocomposite corresponding to CNT content of 0.25 wt%: (a) CNT distribution and (b) matrix meshes.

Fig. 12. Simulation results for CNT nanocomposites: (a) stress-strain curves, damage contours in (b) the matrix and (c) CNT/matrix interface, and (d) maximum principal stress (SP1) inside CNTs corresponding to the failure of the nanocomposite with CNT content of 0.25 wt%. Experimental data is cited from Ref. [30]. 12

Composite Structures 223 (2019) 110921

Q. Liu, et al.

stiffness of each wall is 1000 GPa [30] the Young’s modulus of an equivalent CNT is calculated to be 849 GPa. The effect of CNT’s anisotropy has been examined in [36] showing limited influence. For the purpose of computational efficiency, the anisotropy is not included here. The interfacial strength and fracture energy are assumed to be 75 MPa and 5.0 J/m2 referring to the experimental data from single CNT pull-out experiments in [32]. Fig. 12 shows simulation results for the RVE model in Fig. 11. As shown in Fig. 12(a), the EE-CZM approach leads to good predictions for the stress-strain response of the matrix and the nanocomposites with CNT contents of 0.15 wt% and 0.25 wt%. A ductile damage behaviour is observed for the matrix and the two nanocomposites. The simulation results predict embrittlement of the nanocomposites, which is consistent with the experimental data. With the use of EE-CZM, both microand nano-scale damage mechanisms inside the CNT nanocomposite can be captured. Fig. 12(b) and (c) show damage contours corresponding to the moment of the sharp stress drop on the stress-strain curve, i.e. the failure of nanocomposite. Only damaged domains are shown for clarity and only for the case of 0.25 wt% of CNTs. The nanocomposite fails when a global crack at the micro-scale crosses through the matrix (Fig. 12(b)). The interfacial damage at the nano-scale is also clearly observed in Fig. 12(c), which is more significant around the matrix crack. Interestingly, as indicated by the zoom-in picture in Fig. 12(c), the interfacial damage is easier activated at the CNT ends and locations at which CNTs’ local orientation coincides with the direction of external loading. As shown in Fig. 12(d), the maximum value of the stress inside CNTs is on the order of 20 GPa, which is relatively high but is still below the often reported value of the CNT tensile strength, ∼50 GPa [37]. This supports our assumption of no CNT breakage in the present model.

[2] Lu Z, Yuan Z, Liu Q. 3D numerical simulation for the elastic properties of random fiber composites with a wide range of fiber aspect ratios. Comput Mater Sci 2014;90:123–9. [3] Lu Z, Yuan Z, Liu Q, Hu Z, Xie F, Zhu M. Multi-scale simulation of the tensile properties of fiber-reinforced silica aerogel composites. Mater Sci Eng A 2015;625:278–87. [4] Liu H, Zeng D, Li Y, Jiang L. Development of RVE-embedded solid elements model for predicting effective elastic constants of discontinuous fiber reinforced composites. Mech Mater 2016;93:109–23. [5] Harper L, Qian C, Luchoo R, Warrior N. 3D geometric modelling of discontinuous fibre composites using a force-directed algorithm. J Compos Mater 2017;51:2389–406. [6] Blacklock M, Joosten MW, Pingkarawat K, Mouritz AP. Prediction of mode I delamination resistance of z-pinned laminates using the embedded finite element technique. Compos A Appl Sci Manuf 2016;91:283–91. [7] Pappas G, Canal L, Botsis J. Characterization of intralaminar mode I fracture of AS4/PPS composite using inverse identification and micromechanics. Compos A Appl Sci Manuf 2016;91:117–26. [8] Joosten MW, Dingle M, Mouritz A, Khatibi AA, Agius S, Wang CH. A hybrid embedded cohesive element method for predicting matrix cracking in composites. Compos Struct 2016;136:554–65. [9] Tabatabaei S, Lomov SV, Verpoest I. Assessment of embedded element technique in meso-FE modelling of fibre reinforced composites. Compos Struct 2014;107:436–46. [10] Tabatabaei S, Lomov SV. Eliminating the volume redundancy of embedded elements and yarn interpenetrations in meso-finite element modelling of textile composites. Comput Struct 2015;152:142–54. [11] Tabatabaei SA, Swolfs Y, Wu H, Lomov SV. Full-field strain measurements and meso-FE modelling of hybrid carbon/self-reinforced polypropylene. Compos Struct 2015;132:864–73. [12] Muñoz R, Martínez-Hergueta F, Gálvez F, González C, LLorca J. Ballistic performance of hybrid 3D woven composites: experiments and simulations. Compos Struct 2015;127:141–51. [13] Vorobiov O, Tabatabaei S, Lomov SV. Mesh superposition applied to meso-FE modelling of fibre-reinforced composites: cross-comparison of implementations. Int J Numer Meth Eng 2017;111:1003–24. [14] Romanov V, Lomov SV, Verpoest I, Gorbatikh L. Inter-fiber stresses in composites with carbon nanotube grafted and coated fibers. Compos Sci Technol 2015;114:79–86. [15] Romanov VS, Lomov SV, Verpoest I, Gorbatikh L. Modelling evidence of stress concentration mitigation at the micro-scale in polymer composites by the addition of carbon nanotubes. Carbon 2015;82:184–94. [16] Romanov VS, Lomov SV, Verpoest I, Gorbatikh L. Stress magnification due to carbon nanotube agglomeration in composites. Compos Struct 2015;133:246–56. [17] Savvas D, Stefanou G, Papadopoulos V, Papadrakakis M. Effect of waviness and orientation of carbon nanotubes on random apparent material properties and RVE size of CNT reinforced composites. Compos Struct 2016;152:870–82. [18] Savvas D, Papadopoulos V, Papadrakakis M. The effect of interfacial shear strength on damping behavior of carbon nanotube reinforced composites. Int J Solids Struct 2012;49:3823–37. [19] Rai A, Subramanian N, Chattopadhyay A. Investigation of damage mechanisms in CNT nanocomposites using multiscale analysis. Int J Solids Struct 2017;120:115–24. [20] Romanov VS, Lomov SV, Verpoest I, Gorbatikh L. Can carbon nanotubes grown on fibers fundamentally change stress distribution in a composite? Compos A Appl Sci Manuf 2014;63:32–4. [21] Mehdikhani M, Matveeva A, Aravand MA, Wardle BL, Lomov SV, Gorbatikh L. Strain mapping at the micro-scale in hierarchical polymer composites with aligned carbon nanotube grafted fibers. Compos Sci Technol 2016;137:24–34. [22] Liu Q, Lomov SV, Gorbatikh L. Spatial distribution and orientation of nanotubes for suppression of stress concentrations optimized using genetic algorithm and finite element analysis. Mater Des 2018;158:136–46. [23] Oudjene M, Meghlat E, Ait-Aider H, Lardeur P, Khelifa M, Batoz J. Finite element modelling of the nonlinear load-slip behaviour of full-scale timber-to-concrete composite T-shaped beams. Compos Struct 2018;196:117–26. [24] Papadopoulos V, Impraimakis M. Multiscale modeling of carbon nanotube reinforced concrete. Compos Struct 2017;182:251–60. [25] Nzabonimpa J, Hong W-K, Kim J. Nonlinear finite element model for the novel mechanical beam-column joints of precast concrete-based frames. Comput Struct 2017;189:31–48. [26] Kang J. Composite and non-composite behaviors of foam-insulated concrete sandwich panels. Compos B Eng 2015;68:153–61. [27] Yousefsani SA, Shamloo A, Farahmand F. Micromechanics of brain white matter tissue: a fiber-reinforced hyperelastic model using embedded element technique. J Mech Behav Biomed Mater 2018;80:194–202. [28] Yang L, Yan Y, Liu Y, Ran Z. Microscopic failure mechanisms of fiber-reinforced polymer composites under transverse tension and compression. Compos Sci Technol 2012;72:1818–25. [29] Li D, Yang Q-S, Liu X, He X-Q. Experimental and cohesive finite element investigation of interfacial behavior of CNT fiber-reinforced composites. Compos A Appl Sci Manuf 2017;101:318–25. [30] Zarasvand KA, Golestanian H. Determination of nonlinear behavior of multi-walled carbon nanotube reinforced polymer: experimental, numerical, and micromechanical. Mater Des 2016;109:314–23. [31] Cooper CA, Cohen SR, Barber AH, Wagner HD. Detachment of nanotubes from a polymer matrix. Appl Phys Lett 2002;81:3873–5.

5. Summary and conclusion A novel approach combining the embedded element technique and cohesive zone model is proposed in this work to model interface damage in fibrous composites. The main advantage of the EE-CZM approach is to allow for high-quality structured meshes on both reinforcements and matrix, which is generally difficult to achieve in conventional direct FE models with continuous meshes. This new approach has been systematically compared with the conventional method in several numerical case studies. It is found that the EE-CZM accurately reproduces results of the direct FE modelling, including stress distribution and strain energy in the constituents, as well as interfacial damage development in fibrous composites. The proposed EECZM method is validated against experimental data available for CNT nanocomposites. The latter are challenging to model with the direct FE approach due to the high aspect ratios of the reinforcements. The proposed method not only succeeds in predicting the overall mechanical behaviour of the nanocomposite but also captures the concurrent micro- and nano-scale damage mechanisms inside the nanocomposite. Sensitivity of the EE-CZM approach to the mesh size is thoroughly evaluated, and it is found that the modelling results are sufficiently accurate when the fiber and matrix mesh sizes do not exceed 1 and 12 times of the fiber radius, respectively. The proposed EE-CZM approach opens a new way to model damage development in fibrous composites. Acknowledgements The research leading to these results has received funding from KU Leuven in the framework of the C2 PERMEA project. S.V. Lomov holds Toray Chair for Composite Materials in KU Leuven. References [1] Harper L, Qian C, Turner T, Li S, Warrior N. Representative volume elements for discontinuous carbon fibre composites–Part 1: Boundary conditions. Compos Sci Technol 2012;72:225–34.

13

Composite Structures 223 (2019) 110921

Q. Liu, et al.

[35] Melro A, Camanho P, Pires FA, Pinho S. Micromechanical analysis of polymer composites reinforced by unidirectional fibres: part I-constitutive modelling. Int J Solids Struct 2013;50:1897–905. [36] Romanov V. Modeling tools for micro-scale stress analysis of nano-engineered fiberreinforced composites. University of Leuven; 2015. [37] Gorbatikh L, Wardle BL, Lomov SV. Hierarchical lightweight composite materials for structural applications. MRS Bull 2016;41:672–7.

[32] Barber AH, Cohen SR, Eitan A, Schadler LS, Wagner HD. Fracture transitions at a carbon-nanotube/polymer interface. Adv Mater 2006;18:83–7. [33] Liu Q, Lu Z, Hu Z, Li J. Finite element analysis on tensile behaviour of 3D random fibrous materials: model description and meso-level approach. Mater Sci Eng A 2013;587:36–45. [34] Yuan Z, Lu Z. Numerical analysis of elastic–plastic properties of polymer composite reinforced by wavy and random CNTs. Comput Mater Sci 2014;95:610–9.

14