- Email: [email protected]
BEM multiscale modelling involving micromechanical damage in fibrous composites
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
BEM multiscale modelling involving micromechanical damage in fibrous composites M.L. Velasco, E. Graciani, L. Távara, E. Correa, F. París∗ Grupo de Elasticidad y Resistencia de Materiales, Escuela Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain
a r t i c l e
i n f o
Keywords: BEM Multiscale model Composites Damage Scale effect
a b s t r a c t Composite laminates are materials where different scales are required for the understanding of the damage and for the calculation of structures made of these materials. The appearance of ultra thin plies of composites has opened the possibility of delaying the onset of damage, which at a first stage appears at a micromechanical level, debondings in between fibre and matrix in the weakest lamina of the laminate. A multiscale BEM model is developed in the present paper with the final purpose of being able to study the effect that the relative size of the laminas of the laminate plays in the appearance of this initial damage. The model presents many difficulties derived from the different scales involved in it and from the non-linear nature of the problem under study. The approach followed involves the solution of the whole problem, with the different scales involved in it, at once. The solution obtained is checked with another already obtained with a much simpler model. The multiscale model developed has been proved to be very efficient, accurate and robust, having been applied to simulate the first stages of damage in the light of the scale effect that is trying to be studied.
1. Introduction Composites have evolved in last decades from being an interesting material used in structures of low responsibility in aeronautical components to be the preferred material for primary structures of modern aircrafts. Although composites have, as one of its main features, an excellent behaviour in terms of damage tolerance, it is obvious that the best possible design is that involving no damage in the material in the expected range of load at which the component has to work. With reference to Carbon Fibre Reinforced Plastics (CFRP), the material object of the present study, the weakest laminas, those oriented 90° with relation to the load direction, present a premature presence of damage, typically viewed as transverse through the thickness of the lamina cracks, not involving the breakage of any fibre. The appearance of this damage was experimentally shown to be dependent on the relative thickness of the 90° lamina with regards to the 0° laminas, see for instance Parvizi et al. [1] and Flaggs and Kural [2]. The availability of having ultra thin plies (UTPs) to design and manufacture composite laminates opens the possibility to delay the appearance of this damage, amplifying the range of work of the material. The understanding of the mechanism of damage of these 90° laminas and its dependence on the relative thicknesses of the 0 and 90° laminas becomes of crucial importance to deal with comfort in using these UTPs.
∗
The necessity of limiting the difficulty associated with structural calculations has traditionally led to simplify the representation and modelling of materials. Thus, metals are typically represented as homogeneous (omitting the representation of grains) and isotropic (representing the case of non-preferential orientation of the grains), whereas a lamina of fibre reinforced composites is represented as homogeneous (omitting in this case the representation of fibres and matrix) though orthotropic (to represent the macromechanical behaviour associated with the presence of the fibres aligned along one direction). However, there is a need for many structural engineering problems (in order to fully understand really them), of descending to the scales where the damage/failure/degradation is taking place, it being necessarily compatible with the modelling, at least of a part of the structure, at the scale at which it is studied. This fact, less common in metals, is very typical in composites, in which the calculations to design the composite structure are made at the level of assumed homogeneity, although the damage appears at the level of fibre and matrix and what is more common at the interfaces between them, a level of detail (micromechanical) not covered in the aforementioned homogeneous representation of a ply of a composite. Associated with this question is the fact that a micromechanical study typically requires a very accurate determination of variables (displacements and stresses) associated with the damage under study (debondings between fibres and matrix, cracks connecting debondings,…), in order to evaluate the magnitudes controlling the damage
Corresponding author. E-mail address: [email protected] (F. París).
https://doi.org/10.1016/j.enganabound.2018.03.012 Received 16 October 2017; Received in revised form 22 February 2018; Accepted 16 March 2018 0955-7997/© 2018 Elsevier Ltd. All rights reserved.
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
(the energy release rate of the crack). This calculation has to be additionally carried out on a very reduced area (in terms of several orders of difference in scale) of the whole problem under study. Bearing all previous reasonings in mind, BEM, París and Cañas [3], has the adequate features to deal with a generic problem like this, as will be specified in Section 2. Not too many papers have been published on the use of multiscale models using BEM. Sfantos and Aliabadi [4] used a multiscale approach for modelling material degradation and fracture of polycrystalline ceramics, although in a different context than that used here. Notice that most of problems requiring a double scale of attention were modelled coupling BEM-FEM, using FEM on the zone of detailed analysis and BEM to model zones with no special interest, e.g. Fernandes et al. [5]. In our case, the zone of interest requires a very accurate solution and therefore BEM has to be used for it. BEM has already been applied for composites at mesomechanical levels to study [0,90] symmetric laminates, see París et al. [6,7], the point of attention being the characterization of the damage at this mesomechanical level, all plies being considered as homogeneous materials. The damage was viewed in the model as transverse cracks appearing in the 90° plies and its potential connections with delamination cracks between 0 and 90° laminas, when the transverse cracks are approaching the interface between the laminas. BEM has also been profusely applied for composites at micromechanical level, using a reduced number of fibres, to understand the mechanisms of damage under tension [8], compression [9,10], or biaxial loading [11–13], as well as to understand problems associated to fatigue [14], curing [15,16] or the interaction in between neighbouring fibres [17,18]. An actual unidirectional lamina can be modelled by a random distribution of fibres included in a large matrix. BEM has proved to be an efficient tool to predict failure loads due to the debonding mechanism in a group of fibres, [19], and also to study the crack path (leading to a macro crack) formed by consecutive debonds [20]. In our case, the nature of the problem under consideration will require a multiscale model involving the presence of isotropic and orthotropic behaviour, multibody formulation and capacity of dealing with contact problems, as damage involving debondings and cracks connecting debondings will be simulated. All these features are contemplated in Graciani et al. [21], whose BEM code will be used in this paper. The final objective of this contribution is to build an efficient and accurate multiscale BEM model able to help in the understanding of the damage appearing in the failure of a 90° lamina of a [0,90n]S laminate: the isolated debonding between fibre and matrix and the abandonment of these debondings of the fibre-matrix interface, kinking and entering into the matrix of the lamina.
Fig. 1. Sequential stages of expected damage at micromechanical level in a 90° ply.
Fig. 2. Multiscale 2D model of a [0,90n]S laminate under tension with a cell in the central part of the 90° lamina.
pure mode II, or, what is more plausible, it will kink and penetrate into the matrix, as represented by stage III in Fig. 1. These will be the damage situations investigated in this paper. Fig. 2 represents the problem that is going to be studied, with the damage presented in Fig. 1, by means of a multiscale 2D model. Due to the symmetry of the laminate and to the symmetry of the modelled damage with respect to axis x in Fig. 2, only one half of the laminate needs to be modelled. Two decisions are involved in the model proposed. The first is to take a single debond which is something widely observed experimentally, París et al. [22], and supported by numerical studies, see for instance Panagiotopoulos and Mantič [23], García et al. [24] or Muñoz-Reja et al. [25]. The second decision is to take the single damage symmetric, which allow the symmetry with respect to axis x to be applied. Obviously, dealing with a non homogeneous material, the damage is very much affected by the presence of nearby fibres, but it can be said that the damage represented here is by far the dominant one observed. The generation of this single damage can also be supported by numerical studies, and although its appearance may depend on the failure criterion employed, it has been observed, Mantič et al. [26], that in any case it tends to adopt, even starting at a nonsymmetrical location, a symmetric disposition. Obviously the size of the cell is magnified in Fig. 2 with respect to its actual size in order to have a better perception of the model. As represented in Fig. 2, an external lamina of thickness t0° with the fibres oriented in the direction of the load (0°) will be modelled as a homogeneous material. The central lamina of a total thickness of 2t90° (t90° in Fig. 2) with the fibres oriented 90° with respect to the load direction (normal to the plane represented in the section appearing in Fig. 2), will be mesomechanically modelled as an isotropic material. Finally, the central part of the model will include a cell covering a fibre, where the damage will be considered. This damage, see Fig. 3, will consist in a first phase on a debonding (denoted in Fig. 3 by the semidebonding 𝜃 d ) between fibre and matrix, located nominally normal to the direction of the load although it is condemned in the model, reproducing the physics of the problem, to be contained along the interface between fibre and matrix. The properties of the material systems used in this paper appear in Tables 1 (carbon system) and 2 (glass system). The discussion about the possible strategies to deal with this problem (notice that the representation of the model associated with Fig. 2 is already giving indication about the final decision taken) requires to notice that the damage will be characterized by the Energy Release Rate
2. The problem: modelling strategies The problem under consideration is finally associated with the possibility of delaying the appearance of damage by using ultra thin plies to design and manufacture laminates. This possibility is based on the observed fact that a 90° ply increases its strength when its thickness tends to a minimum value. This value had been restricted in the past to around 120 μm due to technological reasons, which have been solved giving rise to the possibility of having laminas with a thickness of around 40 μm, having only around 5/6 fibres (of a diameter of 7.5 μm as in the case of carbon fibres) through the thickness of the ply. The onset and progression of damage in an isolated fibre of a 90° ply, under tension, are represented in Fig. 1, París et al. [8]. First, based on the interfacial stress state predicted, a minimum interface crack of about 20° will first appear, stage I in Fig. 1. This interface crack would grow till reaching a total debonding angle in the interval of 130–140°, stage II in Fig. 1, where contact zones at the ends of the interfacial crack appear. Finally, either the interfacial crack grows in a stable way under 2
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 3. Damaged cell in the 90° lamina with a debonding between fibre and matrix. Table 1 Carbon-epoxy system used. Material
Properties
0° ply (orthotropic)
E11 = 135 GPa, E22 = 8.75 GPa, E33 = 8.75 GPa, 𝜈 12 = 0.3, 𝜈 13 = 0.3, 𝜈 23 = 0.4, G12 = 4.75 GPa E(E22 = E33 ) = 8.75 GPa, 𝜈 = 0.4 E = 4.2 GPa, 𝜈 = 0.32 E (E22 ) = 15 GPa, 𝜈 = 0.2; R = 3.75 μm
90° ply (isotropic) Matrix (epoxy, isotropic) Fibre (carbon, isotropic)
Fig. 4. Detail of the model involving a debonding between fibre and matrix and an associated kinked crack.
representative enough to get accurate results. If the area where the debonding takes place is separated from the interface connection, the model for the second step needs to be changed, and a multifibre model for this second step will probably be required. On the contrary, if strategy ST2 is used, the main challenge will be to make compatible the required refinements of the discretization at the two scales involved in the model, avoiding having very similar boundary integral equations that may lead to a bad conditioned final system of equations. Obviously, a computer with more capacity is also required for employing this strategy, which was that finally employed. A similar approach was considered the best option for the analysis of local debondings in analysing composite components, Reinoso et al. [27], where both strategies were compared. Notice that in that case the same scale was required and only the refinement of one part of the component where the damage may appear called for the consideration of the two strategies already mentioned. In the present case, the difficulty is much more extreme as the refinement has to be performed at the microscale level, where the damage appears, combined with a mesoscale model of the rest of the problem. Once the first stage of the expected damage has been identified, it is reasonable to accept that the generation, from these debondings, of a mesocrack through the thickness of a 90° ply, will involve for the debonding crack the abandon of the interface between fibre and matrix, kinking into the matrix [8]. Thus, additionally to the cell already schematized in Fig. 3, a cell as that presented in Fig. 4, showing the presence of a debonding plus a kinked crack entering into the matrix must be studied too. Notice that this cell plays an additional role as that of the previous cell as it covers also the third stage of damage in a 90° lamina of a laminate, represented in Fig. 1. The used self-developed BEM code (which includes all the features necessary to use strategy ST2) is written in Fortran 90 and compiled in a Linux based system. This code is compiled to run in a vectorized manner i.e. able to use multiple processors within a High Performance Computing (HPC) System.
Table 2 Glass-epoxy system used. Material
Properties
Matrix (epoxy, isotropic) Fibre (glass, isotropic)
E = 2.79 GPa, 𝜈 = 0.33 E = 70.9 GPa, 𝜈 = 0.22; R = 7.5 μm
(ERR) of the debonding crack in accordance with interfacial Fracture Mechanics (IFM) [8], which in turns requires to calculate with high accuracy displacements and stresses along a certain circumferential interval (0.5° in this paper) at both sides of the crack tip, i.e. along a distance of 0.098 μm. These dimensions coexist in the model with the thickness of a lamina (oscillating in between 34 and 680 μm) and the length of the laminate (taken as 5000 μm). Notice also the non linear character of the problem, derived from the presence of potential contact zones that may appear between fibre and matrix. There are two strategies to deal with the problem under consideration from a modelling and numerical point of view, and will be designated as ST1 and ST2 in what follows. Strategy ST1: the problem is solved in two steps. In the first step the area of interest (of micromechanical dimensions) where the damage is going to take place is represented with a no very refined discretization, or even the micromechanical representation can be omitted and substituted by a homogeneous material with the same properties that the rest of the lamina. In a second step, the values of stresses or displacements are used as boundary conditions for the isolated cell where fibre, matrix and damage are now represented, Fig. 3. Strategy ST2: the problem is solved directly in one step, modelling in a very refined way the zone where interfacial fracture mechanics calculation will take place in conjunction with areas where refined discretizations are not required (as those as the extreme of laminates where the external action is applied). In the case of using ST1 it is obvious that the problems associated to step 1 and step 2 are reasonably easy to solve (with no too great differences in the size of the elements required on each step), although an iterative procedure, as the solution obtained in step 2 may not be compatible with that obtained in step 1, is required. With reference to this ST1 strategy, if the interface of connection between the first and second problem is not very well discretized, as the damage modelled in the second problem (the debonding between fibre and matrix) is very close to the line of separation between the two models, the results obtained in the second step may be contaminated by the poor discretization of the interface connection carried out in the first step, not being then
3. The BEM model The discretization of the multiscale model is a key factor to get accurate results in the area of interest. In order to evaluate the energy release rate associated to crack propagation, Irwin’s approach will be employed. Therefore, the area of interest is the neighbourhood of the tip of the crack representing the damage under consideration. As has been mentioned above, two different cases will be shown: the first one containing a debonding crack between fibre and matrix and the second one containing a matrix crack initiating from the debond crack and running normal to the applied load. 3
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 5. General view of the discretization carried out in the model.
As will be shown below, a virtual crack closure length of 3.27•10–2 μm (corresponding to a subtended angle in the interface of the fibre with the matrix of 0.5 °) will be employed in the calculation of the ERR which will be calculated, by using a numerical integration scheme in the range 6.54⋅10–6 μm < 𝜌<3.27⋅10–2 , with 𝜌 being the distance to the crack tip. Consequently, smaller elements in the vicinity of the crack tip should be at least one order of magnitude shorter than the lower limit of 𝜌, since linear elements are used and the presence of the singularity spoils the accuracy in the closest vicinity of the crack tip. The discretization of the rest of the model, as the results obtained associated to this zone are not of primary interest, must be focused on reducing as much as possible the number of nodes and on avoiding negative effects on the system of equations, as would be the presence of corners in some of the boundaries (those of the cells for instance) or the presence of boundaries near the damaged area where the calculations must be extremely accurate. Thus, regular rules used in BEM, París and Cañas [3], have guided the discretization away from the crack tip. Fig. 5 represents a scheme of the discretization carried out for a case having a debond crack. Four solids are involved in the model: S1 is the fibre, S2 is the matrix surrounding the fibre in the local cell, S3 is the 90° lamina and S4 is the 0° lamina. Each sector of the boundaries is identified as Sij indicating that is the j-th sector of solid Si . The additional information about using a uniform discretization (UD) or a progressive discretization (PD) appears in between brackets. In this particular case of a progressive discretization the length ratio in consecutive elements has been maintained between 1 and 1.2. Details of the mesh in the local scale will be shown in the next subsections in which the two different local cells will be described separately. The element length (in microns) is also indicated at each corner and in the zones in which a uniform discretization is used. Element length is identical in all elements sharing a corner node and conforming discretizations are employed in all contact zones or bonded interfaces. Notice that the ratio between the size of the model (5⋅103 μm) and the size of the smallest element (∼10–7 μm) is extremely high. Since the accurate evaluation of the boundary integrals when the collocation point is extremely far away from the smallest element is impossible, the use of a local and global scale is mandatory even if the material of solids S2 and S3 was identical. The use of the local model reduces the length ratio between the largest and the smallest elements of each solid allowing an accurate evaluation of the boundary integrals to be carried out. Finally, perfect bonding conditions are employed in the interface between S3 and S4 and between S3 and the local cells described in the following.
The most important information associated with the discretization carried out in the local cell is that regarding to the mesh employed in the neighbourhood of the crack tip. The debonding crack, being an interface crack, can be modelled in accordance with either an open model, Williams [28], which would involve interpenetrations near the crack tip of the faces of the crack, or in accordance with a contact model, Comninou [29], which would involve the appearance of a contact zone near the crack tip. In our case the BEM code employed allows the presence of a contact zone to be developed and consequently BEM solution directly satisfies the conditions of Comninou’s model, given that if interpenetrations tend to appear (using initially open conditions), the iterative procedure employed will adjust the incompatibilities found. Notice that the contact zone to avoid interpenetrations may be extremely small or relatively large depending on the fracture mode mixity, which varies with the length of the crack. However, this fact does not affect the discretization to be carried out, since the evaluation of the ERR has to be carried out at a scale along a distance of growth of physical relevance. Consequently if this distance scale is larger than the contact zone, the presence of extremely small contact zone does not affect the ERR associated to crack propagation at the scale of interest. Since friction between crack faces is neglected, energy release rate G can be evaluated according to Irwin’s Virtual Crack Closure Technique, and decomposed into fracture mode I contribution, GI , and fracture mode II contribution, GII : 𝐺 = 𝐺𝐼 + 𝐺𝐼𝐼 = +
1 Δ𝜃𝑑 ∫0
Δ𝜃𝑑
1 Δ𝜃𝑑 ∫0
Δ𝜃𝑑
[
( ) ( )] 𝜎𝑟𝑟 𝑅, 𝜃𝑑 + 𝜃 ′ Δ𝑢𝑟 𝑅, 𝜃𝑑 + Δ𝜃𝑑 − 𝜃 ′ 𝑑𝜃 ′
[ ( ) ( )] 𝜎𝑟𝜃 𝑅, 𝜃𝑑 + 𝜃 ′ + Δ𝑢𝜃 𝑅, 𝜃𝑑 + Δ𝜃𝑑 − 𝜃 ′ 𝑑𝜃 ′ ,
(1)
where (r,𝜃) is a local co-ordinate system centred with the fibre, 𝜎 rr and 𝜎 r𝜃 are, respectively, the normal and tangential tractions along the interface, while Δur and Δu𝜃 are, respectively, the opening and sliding between crack faces. The previous experience indicates, París et al. [8], that depending on the debonding angle, the interface crack fracture mode evolves from having an extremely small contact zone (for small debondings 𝜃 d < 60°) to showing a significant contact zone (for large debondings 𝜃 d < 60°). When the length of the contact zone is larger than Δa = RΔ𝜃 d (with R being the fibre radius), the scale in which the ERR is being evaluated, pure fracture mode II takes place. Fig. 7 shows a scheme of the discretization carried out (the same for all the cases studied) at both sides of the crack tip and Table 3 gives the values of the angle 𝛼 ij subtended by the element joining nodes i and j, 𝜆ij denoting the length of the element connecting nodes i and j. Since the accuracy of the solution in the closest vicinity of the crack tip is spoiled by the presence of the singularity, the solution of the closest nodes to the crack tip is not employed, and integration is carried out from node 5 up to node 16 (see Table 3), which means that the integrals in (1) are effectively carried out in the range 0.0001°<𝜃′<Δ𝜃 d = 0.5°.
3.1. Interface crack case Fig. 6 represents a scheme of the discretization carried out in the local cell for a case having only a debond crack. Notation employed follows the same rules described in the previous section. 4
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 6. Discretization of the local cell corresponding to a debonding case. Table 3 Subtended central angle and length of the elements in the neighbourhood of the crack tip.
𝛼 ij (0 ) 𝜆ij (μm)
𝛼 ij (0 ) 𝜆ij (μm)
1
2
3
4
5
6
7
8
6.25·10– 6 4.09·10– 7
6.25·10– 6 4.09·10– 7
1.25·10– 5 8.18·10– 7
2.5·10– 5 1.63·10– 6
5·10– 5 3.27·10– 6
2.04·10– 4 1.33·10– 5
4.16·10– 4 2.72·10– 5
8.5·10– 4 5.56·10– 5
9
10
11
12
13
14
15
16
1.73·10– 3 1.13·10– 4
3.53·10– 3 2.31·10– 4
7.21·10– 3 4.72·10– 4
1.47·10– 2 9.62·10– 4
3·10– 2 1.96·10– 3
6.12·10– 2 4·10– 3
1.25·10– 1 8.17·10– 3
2.55·10– 1 1.67·10– 2
open or partially open in all cases (different debondings). For the different debond angles studied, the mesh is identical to that shown in Fig. 6, corresponding to the 𝜃 d = 60° case, except for the fact that the mesh in the neighbourhood of the crack tip is shifted to its corresponding position. When the HPC system already described is used, the computation times are quite reasonable, being of the same order independent from the configuration studied (different length of the debonding). In the cases in which a contact zone is detected, only a few iterations are needed to obtain convergence and only a reduced number of contact equations has to be rebuilt in each iteration. To give an idea of the complexity of the problem, with reference to the case of a debonding of 40°, the total number of elements is 2070, leading to a total system of equations of 4728, with a minimum size of an element of 4,09·10–7 μm and a maximum size of element of 50 μm. When this model is run in the HPC system, which includes 8 cpus, an iteration takes around 40 s. Usually, to get a converged solution of the whole problem, only a reduced number of iterations (9 in the case described here) is needed, the total solution time being around 4–5 min. It is well known how sensitive is BEM, particularly in the presence of singularities, to the discretization carried out, and although the final proof of having a correct discretization will be obtained in the next section comparing with previous results, it has to be mentioned here that this study is not carried out with prediction purposes but with the aim to have a model to compare results associated to different configurations of a laminate. Thus, once the discretization employed is proved to work properly, it will be employed for all cases, then supporting the validity of the comparison for the problem under consideration.
Fig. 7. Discretization performed in the neighbourhood of the crack tip to analyse the first phase of damage (isolated debondings between fibre and matrix).
Note that the problem under consideration is nonlinear, as it implies the contact/separation of surfaces of fibre and matrix. This fact does not interfere with the intrinsic difficulty of solving of the problem, but it is significant during the process of definition of a suitable discretization. Thus, the crack was artificially closed along this process, adjusting the discretization, in order to get a smooth evolution of stresses in the neighbourhood of the crack tip. Once the discretization had been adjusted with the crack closed, it worked properly with faces of the crack
3.2. Kinked matrix crack case The local cell described in the previous subsection corresponds to the modelling of the first stage of damage, that is, the debonding 5
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 9. Predictions on G values obtained with the multiscale model (surrounded by matrix) and with the single fibre model.
Fig. 8. Discretization of the local cell corresponding to the kinked matrix crack case.
of a single fibre surrounded by matrix, both adimensionalized using a reference number G0 , París et al. [8]. Only the values of the total energy release rate are shown for the sake of clarity. It can be clearly observed that they are almost undistinguishable, proving then the validity of the multiscale model developed.
between fibre and matrix. This damage is assumed to continue, once the debonding crack stops, kinking and entering into the matrix following a direction normal to the applied load, what implies for the kinked crack to work under mode I, París et al. [8]. The model described in the previous sections has been adapted, as shown in Fig. 8, for the study of stage III of damage represented in Fig. 1. Although at the global scale, the models representing the two stages of damage have different sizes, the details of the mesh at the global scale of the model corresponding to the third stage identified in Fig. 1 are omitted for the sake of brevity, since the mesh is made following a similar scheme to that described in Fig. 5 for the model of the first and second stages. In the local scale, in order to easily assess stresses ahead of the crack tip, a fictitious interface along the predicted path of the kinked crack has been created, starting at the position of the original crack tip of the debonding and aligned normally to the applied load. Then, for a particular length of the assumed kinked crack, the closest neighbourhood of its crack tip is treated as described in the previous section. Contact conditions are employed between crack faces (from the symmetry plane to the crack tip) and perfect bonding conditions are applied along the fictitious interface where the crack has not extended. Other feasible option to study how the interface crack gets into the matrix (homogeneous material) without including a fictitious interface could be the use of the so-called Dual-BEM [30] or the use of a Symmetric Galerkin Scheme [31,32]. However, due to the fact that in the problem under study the direction of progression of the crack is a priori known, the method selected, in spite of the necessity of using a fictitious interface, was estimated to be the optimum solution. Moreover, the use of those schemes (dual BEM or SGBEM) would imply the use of hypersingular integral equations as well as their corresponding fundamental solutions.
4.2. Multiscale model results Once the validity of the model has been proved, one of the questions that attracts immediate attention is the shape of the hexagonal cell. Hexagonal cells are typically used when representing composites micromechanically as they are reproducible leading to a complete representation of the material at a plane normal to the fibres direction, allowing also in a simple way to reproduce the volume fibre ratio at the level of the cell. However, the fibres are in fact distributed in the material in a complete random way, and consequently any shape of the cell would be acceptable as will represent the actual situation of a certain fibre in the material. As the hexagonal cell has corners, this fact may introduce some alterations in the solution. Notice that in the case of having the corner, as in Section 4.1, located in a unique media, it would not produce any effect when using a FEM approach, but would do when using a BEM approach. However, having the corner of the cell in between dissimilar media (with matrix at one side and the homogeneous model of the 90° ply at the other side), it produces a nominal singular stress state associated to a multimaterial corner, Barroso et al. [33], BEM being much more sensitive than FEM to this situation. Then, an alternative circular cell is now going to be used in order to check if the presence of corners in the hexagonal cell produces significant alterations in the results. Figs. 10 (for a 30° debonding) and 11 (for a 50° debonding), represent the evolutions of the radial stress 𝜎 r ahead of the crack tip and the relative displacements Δur of the faces of the crack near the crack tip, values used to calculate G. It can be clearly observed the similarity of the results found for the two cells, hexagonal (HC) and circular (CC). As a general view of the situation, Fig. 12 represents the values of G predicted with the use of both shapes of the cell for the values of G and their counter parts GI and GII . When comparing the predictions of the results obtained with the two shapes of the cell, it has to be mentioned that they correspond obviously to different discretizations of the boundary of the cell, as uniform discretizations are performed for the circular cell, whereas a progressive discretization (involving more elements) is used in the case of hexagonal cells near the corners of the boundary of the cell. In any case the results are quite similar, having taken in Figs. 10 and 11 two representative situations involving one of them the case of greatest difference in the predictions. All of this supports the use of the BEM model initially represented in Figs. 2 and 3 to study the scale effect in composites by
4. Results 4.1. Validation of the model The key result, as will be used in further studies to discriminate on the role of the scale effect in the damage of composite laminates, is the energy release rate G. As the first question to be addressed is the ability of the multiscale model developed to get accurate results, those obtained with the present model will be compared with the ones obtained with a simpler model, a single fibre surrounded by matrix, París et al. [8]. To make the comparison effective, all the bodies, except for the fibre, appearing in the multiscale model will adopt the properties of the matrix, thus making the comparison advisable. Fig. 9 represents the values of G obtained with the multiscale model and those obtained with the model 6
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 10. Radial stress and relative displacements in the neighbourhood of the crack tip for a debonding of 30° using a hexagonal cell (HC) and a circular cell (CC).
Fig. 13. Evolution of G, and their counterparts GI and GII , for different lengths of the kinked crack.
8, will now be considered. The case corresponding to a debonding of 65° has been selected, as this debonding is in the interval where, as justified in París et al. [8], the abandonment by the debonding crack of the interface fibre-matrix, kinking into the matrix, is more plausible. Different lengths of the kinked crack have been considered, Fig. 13 representing the evolution of the values of G of the kinked crack, and their counterparts GI and GII , as a function of its length for a carbon epoxy-system. First of all, the evolution obtained with the multiscale model is comparable to that obtained by París et al. [8] for a glass-epoxy system, using a single fibre model surrounded by matrix. Several things are noticeable in this evolution. First, coherently with the decision taken to build the model, negligible values of GII are obtained, leading to almost a pure mode I associated to the kinked crack. This is coherent with the reasoning that led in París et al. [8] to select the 60–70° range for the debonding as that more plausible for the kinking of the interface crack to occur. Additionally, it is also noticeable the increasing character of G (and obviously of GI ) with the increase of the length of the kinked crack. This leads to conclude that once the kinking of the debonding crack arises, the growth of the kinked crack will be unstable, based on the increasing character of G and on the unaltered character of the fracture mode of the kinked crack (mode I) along its path. Notice that a maximum kinked crack length of 2R has been considered in Fig. 13, although further extensions of the crack have been explored confirming the same behaviour. This study has required the use of a bigger than in previous cases. The same cell has been used for all cases considered involving a kinked crack of different length, which makes the results comparable among them. The aforementioned lack of stability of the growth of the kinked crack in the matrix detected by the present numerical model is in agreement with the experimental evidence presented in Correa et al. [34]. It is also important to have a look of the deformed shape of the model obtained involving the kinked crack, Fig. 14, where a magnification factor of 1.2 has been employed. In addition to the coherence of the deformed shape obtained (only the fibre and surrounding matrix are represented in order to concentrate in the area of interest), what can be observed is that independently of the length of the kinked crack chosen (the deformed shape represented corresponds to a length of the kinked crack of 2R (7.5 μm)), no interaction between the lips of the two cracks is detected. As soon as a kinked crack is modelled, no matter how short it is, the two faces of the parent debonding crack separate along the whole length, and thus the contact zone that existed prior to kinking for this particular interface crack size (65°) disappears. This result is also in agreement with París et al. [8].
Fig. 11. Radial stress and relative displacements in the neighbourhood of the crack tip for a debonding of 50° using a hexagonal cell (HC) and a circular cell (CC).
Fig. 12. Values of G, GI and GII for the whole range of debondings using hexagonal cells (HC) and circular cells (CC).
means of the energy release rate of the debonding crack, as the results are not affected by the shape of the cell. Once the presence of corners in the cell has been proved to have no influence in the results of interest, it justifies the use of a rectangular cell for the case of including the kinked crack in the model. Once the ability of the multiscale model to accurately predict the values of the parameter used to study the scale effect has been proved, the model generated to study the third phase of damage, Figs. 4 and 7
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
age appears at a micromechanical level (debonding in between fibre and matrix at a first stage) a multiscale model is required. BEM has been found to be the best option to deal with the problem under consideration, using one single model where both scales are modelled at the same time. The problem becomes very complicated, numerically speaking, due mainly to the scales used to represent the laminate, with a huge difference in the size of the elements involved in the model. Additionally, non-linearities at micro simulation level appear, derived from the consideration of the damage at the microscale, i.e. the aforementioned debondings between fibre and matrix and the kinking of these debond cracks entering into the matrix. The model has proved to be extremely accurate, based on the comparison with previous results and supported by experimental evidence. In addition, the ability of the model to deal with the first stages of the damage as a function of the variables that define the scale effect in composites has also been proved to be very robust and agile. Further studies will be carried out, París et al. [35], to check other predictions based on the multiscale numerical model developed in this paper. Additionally, the model is able to cover more complicated cells involving more than one fibre, as well as to study further stages of the damage involved in the scale effect.
Fig. 14. Undeformed and deformed shapes of the model (only the fibre and surrounding matrix are represented).
Acknowledgements This work was supported by the Spanish Ministry of Education, Culture and Sports under Grants MAT2016-80879-P and MAT2013-45069P. References [1] Parvizi A, Garret KV, Bailey JE. Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates. J Mater Sci 1978;13:195–201. [2] Flaggs DL, Kural MH. Experimental determination of the in situ transverse lamina strength in graphite/epoxy laminates. J Compos Mater 1982;16:103–16. [3] París F, Cañas J. Boundary element method. Fundamentals and applications. Oxford: OUP; 1997. [4] Sfantos GK, Aliabadi MH. Multi-scale boundary element modelling of material degradation and fracture. Comput Meth Appl Mech Eng 2007;196:1310–29. [5] Fernandes GR, Pituba JJC, de Souza Neto EA. FEM/BEM formulation for multi-scale analysis of Stretched Plates. Eng Anal Bound Elem 2015;54:47–59. [6] París F, Blázquez A, Mccartney LN, Mantič V. Characterization and evolution of matrix and interface related damage in [0/90]S laminates under tension. Part I: numerical predictions. Compos Sci Technol 2010;70:1168–75. [7] París F, Blázquez A, Mccartney LN, Barroso A. Characterization and evolution of matrix and interface related damage in [0/90]S laminates under tension. Part II: experimental evidence. Compos Sci Technol 2010;70:1176–83. [8] París F, Correa E, Mantič V. Kinking of transversal interface cracks between fiber and matrix. J Appl Mech 2007;74:703–16. [9] Correa E, Mantič V, París F. A micromechanical view of inter-fibre failure of composite materials under compression transverse to the fibres. Compos Sci Technol 2008;68:2010–21. [10] Correa E, París F, Mantič V. Numerical characterisation of the fibre-matrix interface crack growth in composites under transverse compression. Eng Fract Mech 2008;75:4085–103. [11] Correa E, París F, Mantič V. Effect of the presence of a secondary transverse load on the inter-fibre failure under tension. Eng Fract Mech 2013;103:174–89. [12] Correa E, París F, Mantič V. Effect of a secondary transverse load on the inter-fibre failure under compression. Compos Part B Eng 2014;65:57–68. [13] Mantič V, Távara L, Blázquez A, Graciani G, París F. A linear elastic-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads. Int J of Fract 2015;195:15–38. [14] Correa E, Gamstedt EK, París F, Mantič V. Effects of the presence of compression in transverse cyclic loading on fibre–matrix debonding in unidirectional composite plies. Compos Part A 2007;38:2260–9. [15] Correa E, Mantič V, París F. Effect of thermal residual stresses on matrix failure under transverse tension at micromechanical level: a numerical and experimental analysis. Compos Sci Technol 2011;71:622–9. [16] Correa E, París F, Mantič V. Effect of thermal residual stresses on the matrix failure under transverse compression at micromechanical level – a numerical and experimental study. Compos Part A Appl Sci Manuf 2012;43:87–94. [17] Sandino C, Correa E, París F. Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load. Eng Fract Mech 2016;168:58–75. [18] Muñoz-Reja M, Távara L, Mantič V, Cornetti P. Crack onset and propagation at fibre–matrix elastic interfaces under biaxial loading using finite fracture mechanics. Compos Part A – Appl Sci Manuf 2016;82:267–78.
Fig. 15. Values of G, GI and GII for different values of the t0 ° thickness.
4.3. Incidence of the 0° ply thickness in the scale effect Once the ability of the multiscale model developed to work satisfactorily has been proved, a first application of the model to the problem to which the model has been developed for, that is the scale effect in composites, is now going to be developed. The scale effect is sometimes referred to by means of the thickness of the 90° ply, Parvizzi et al. [1], or by a number n representing the proportion of 90° plies with regards to 0° plies in a [0,90n]S laminate, Flaggs and Kural [2]. Three different thicknesses, t0° = 85 μm, 170 μm and 340 μm, which correspond to n = 0.5, 1 and 2, will be considered in this study. The values of G, GI and GII for the three thicknesses mentioned and for the whole interval of debondings considered are represented in Fig. 15, with reference to the multiscale model considering only the first stage of damage. The independence of the values of G, GI and GII observed for the three thicknesses modelled, clearly indicates that the value of n does not play any role in the development of the first phase of damage (stages I and II in Fig. 1), that involving isolated debondings in between fibre and matrix.
5. Conclusions A multiscale model to deal with the onset and propagation of damage in composite laminates, basically those under the denomination [0,90n]S , has been developed. The aim of the development of the model is to have a tool to analyse the scale effect in these laminates, i.e. the onset of damage at different load values when the thickness of the 90° lamina or the value of n changes in the laminate. As the scale effect is controlled by one parameter at the scale of the laminate and the dam8
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
[19] Távara L, Mantič V, Graciani G, París F. Modelling interfacial debonds in unidirectional fibre-reinforced composites under biaxial transverse loads. Compos Struct 2016;136:305–12. [20] Távara L, Mantič V. Crack paths formed by multiple debonds in LFRP composites. Mech Res Commun 2017;84:148–54. [21] Graciani E, Mantič V, París F, Blázquez A. Weak formulation of axi-symmetric frictionless contact problems with boundary elements. Application to interface cracks. Comput Struct 2005;83:836–55. [22] París F, Correa E, Cañas J. Micromechanical view of failure of the matrix in fibrous composite materials. Compos Sci Tech 2003;63:1041–52. [23] Panagiotopoulos CG, Mantič V. Symmetric and nonsymmetric debonds at fiber-matrix interface under transverse loads. An application of energetic approaches using collocation BEM. Anales de Mecánica de la Fractura 2013;30:125–30. [24] García I, Mantič V, Graciani E. Debonding at the fibre–matrix interface under remote transverse tension. One debond or two symmetric debonds? Eur J Mech A/Solids 2015;53:75–88. [25] Muñoz-Reja M, Távara L, Mantič V. Symmetrical or non-symmetrical debonds at fiber–matrix interfaces: a study by bem and finite fracture mechanics on elastic interfaces. J. Multiscale Model 2017;08:1740008. [26] Mantič V, Távara L, Blázquez A, Graciani E, París F. A linear elasti-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads. Int J Fract 2015;195:15–38.
[27] Reinoso J, Blázquez A, Estefani A, París F, Cañas J. A composite runout specimen subjected to tension-compression loading conditions: experimental and global-local finite element analysis. Compos Struct 2013;101:274–89. [28] Williams ML. The stress around a fault or crack in dissimilar media. Bull Seismol Soc Am 1959;49:199 -24. [29] Comninou M. The interface crack. J Appl Mech 1977;44:631–6. [30] Aliabadi MH, Saleh AL. Fracture mechanics analysis of cracking in plain and reinforced concrete using boundary element method. Eng Fract Mech 1998;93:115–44. [31] Sutradhar A, Paulino G, Gray L. Symmetric Galerkin boundary element method. Berlin: Springer; 2008. [32] Távara L, Mantič V, Salvadori A, Gray LJ, París F. Cohesive-zone-model formulation and implementation using the symmetric Galerkin boundary element method for homogeneous solids. Comput Mech 2013;51:535–51. [33] Barroso A, Mantič V, París F. Singularity analysis of anisotropic multimaterial corners. Int J of Fract 2003;119:1–23. [34] Correa E, Valverde MI, Velasco ML. París, F. Microscopical observations of inter-fibre failure under tension. Compos Sci Technol 2018;155:213–20. [35] París F, Velasco ML, Correa E. Micromechanical study of the influence of scale effect in the first stage of damage in composites. Compos Sci Technol 2018;160:1–8.
9
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
BEM multiscale modelling involving micromechanical damage in fibrous composites M.L. Velasco, E. Graciani, L. Távara, E. Correa, F. París∗ Grupo de Elasticidad y Resistencia de Materiales, Escuela Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain
a r t i c l e
i n f o
Keywords: BEM Multiscale model Composites Damage Scale effect
a b s t r a c t Composite laminates are materials where different scales are required for the understanding of the damage and for the calculation of structures made of these materials. The appearance of ultra thin plies of composites has opened the possibility of delaying the onset of damage, which at a first stage appears at a micromechanical level, debondings in between fibre and matrix in the weakest lamina of the laminate. A multiscale BEM model is developed in the present paper with the final purpose of being able to study the effect that the relative size of the laminas of the laminate plays in the appearance of this initial damage. The model presents many difficulties derived from the different scales involved in it and from the non-linear nature of the problem under study. The approach followed involves the solution of the whole problem, with the different scales involved in it, at once. The solution obtained is checked with another already obtained with a much simpler model. The multiscale model developed has been proved to be very efficient, accurate and robust, having been applied to simulate the first stages of damage in the light of the scale effect that is trying to be studied.
1. Introduction Composites have evolved in last decades from being an interesting material used in structures of low responsibility in aeronautical components to be the preferred material for primary structures of modern aircrafts. Although composites have, as one of its main features, an excellent behaviour in terms of damage tolerance, it is obvious that the best possible design is that involving no damage in the material in the expected range of load at which the component has to work. With reference to Carbon Fibre Reinforced Plastics (CFRP), the material object of the present study, the weakest laminas, those oriented 90° with relation to the load direction, present a premature presence of damage, typically viewed as transverse through the thickness of the lamina cracks, not involving the breakage of any fibre. The appearance of this damage was experimentally shown to be dependent on the relative thickness of the 90° lamina with regards to the 0° laminas, see for instance Parvizi et al. [1] and Flaggs and Kural [2]. The availability of having ultra thin plies (UTPs) to design and manufacture composite laminates opens the possibility to delay the appearance of this damage, amplifying the range of work of the material. The understanding of the mechanism of damage of these 90° laminas and its dependence on the relative thicknesses of the 0 and 90° laminas becomes of crucial importance to deal with comfort in using these UTPs.
∗
The necessity of limiting the difficulty associated with structural calculations has traditionally led to simplify the representation and modelling of materials. Thus, metals are typically represented as homogeneous (omitting the representation of grains) and isotropic (representing the case of non-preferential orientation of the grains), whereas a lamina of fibre reinforced composites is represented as homogeneous (omitting in this case the representation of fibres and matrix) though orthotropic (to represent the macromechanical behaviour associated with the presence of the fibres aligned along one direction). However, there is a need for many structural engineering problems (in order to fully understand really them), of descending to the scales where the damage/failure/degradation is taking place, it being necessarily compatible with the modelling, at least of a part of the structure, at the scale at which it is studied. This fact, less common in metals, is very typical in composites, in which the calculations to design the composite structure are made at the level of assumed homogeneity, although the damage appears at the level of fibre and matrix and what is more common at the interfaces between them, a level of detail (micromechanical) not covered in the aforementioned homogeneous representation of a ply of a composite. Associated with this question is the fact that a micromechanical study typically requires a very accurate determination of variables (displacements and stresses) associated with the damage under study (debondings between fibres and matrix, cracks connecting debondings,…), in order to evaluate the magnitudes controlling the damage
Corresponding author. E-mail address: [email protected] (F. París).
https://doi.org/10.1016/j.enganabound.2018.03.012 Received 16 October 2017; Received in revised form 22 February 2018; Accepted 16 March 2018 0955-7997/© 2018 Elsevier Ltd. All rights reserved.
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
(the energy release rate of the crack). This calculation has to be additionally carried out on a very reduced area (in terms of several orders of difference in scale) of the whole problem under study. Bearing all previous reasonings in mind, BEM, París and Cañas [3], has the adequate features to deal with a generic problem like this, as will be specified in Section 2. Not too many papers have been published on the use of multiscale models using BEM. Sfantos and Aliabadi [4] used a multiscale approach for modelling material degradation and fracture of polycrystalline ceramics, although in a different context than that used here. Notice that most of problems requiring a double scale of attention were modelled coupling BEM-FEM, using FEM on the zone of detailed analysis and BEM to model zones with no special interest, e.g. Fernandes et al. [5]. In our case, the zone of interest requires a very accurate solution and therefore BEM has to be used for it. BEM has already been applied for composites at mesomechanical levels to study [0,90] symmetric laminates, see París et al. [6,7], the point of attention being the characterization of the damage at this mesomechanical level, all plies being considered as homogeneous materials. The damage was viewed in the model as transverse cracks appearing in the 90° plies and its potential connections with delamination cracks between 0 and 90° laminas, when the transverse cracks are approaching the interface between the laminas. BEM has also been profusely applied for composites at micromechanical level, using a reduced number of fibres, to understand the mechanisms of damage under tension [8], compression [9,10], or biaxial loading [11–13], as well as to understand problems associated to fatigue [14], curing [15,16] or the interaction in between neighbouring fibres [17,18]. An actual unidirectional lamina can be modelled by a random distribution of fibres included in a large matrix. BEM has proved to be an efficient tool to predict failure loads due to the debonding mechanism in a group of fibres, [19], and also to study the crack path (leading to a macro crack) formed by consecutive debonds [20]. In our case, the nature of the problem under consideration will require a multiscale model involving the presence of isotropic and orthotropic behaviour, multibody formulation and capacity of dealing with contact problems, as damage involving debondings and cracks connecting debondings will be simulated. All these features are contemplated in Graciani et al. [21], whose BEM code will be used in this paper. The final objective of this contribution is to build an efficient and accurate multiscale BEM model able to help in the understanding of the damage appearing in the failure of a 90° lamina of a [0,90n]S laminate: the isolated debonding between fibre and matrix and the abandonment of these debondings of the fibre-matrix interface, kinking and entering into the matrix of the lamina.
Fig. 1. Sequential stages of expected damage at micromechanical level in a 90° ply.
Fig. 2. Multiscale 2D model of a [0,90n]S laminate under tension with a cell in the central part of the 90° lamina.
pure mode II, or, what is more plausible, it will kink and penetrate into the matrix, as represented by stage III in Fig. 1. These will be the damage situations investigated in this paper. Fig. 2 represents the problem that is going to be studied, with the damage presented in Fig. 1, by means of a multiscale 2D model. Due to the symmetry of the laminate and to the symmetry of the modelled damage with respect to axis x in Fig. 2, only one half of the laminate needs to be modelled. Two decisions are involved in the model proposed. The first is to take a single debond which is something widely observed experimentally, París et al. [22], and supported by numerical studies, see for instance Panagiotopoulos and Mantič [23], García et al. [24] or Muñoz-Reja et al. [25]. The second decision is to take the single damage symmetric, which allow the symmetry with respect to axis x to be applied. Obviously, dealing with a non homogeneous material, the damage is very much affected by the presence of nearby fibres, but it can be said that the damage represented here is by far the dominant one observed. The generation of this single damage can also be supported by numerical studies, and although its appearance may depend on the failure criterion employed, it has been observed, Mantič et al. [26], that in any case it tends to adopt, even starting at a nonsymmetrical location, a symmetric disposition. Obviously the size of the cell is magnified in Fig. 2 with respect to its actual size in order to have a better perception of the model. As represented in Fig. 2, an external lamina of thickness t0° with the fibres oriented in the direction of the load (0°) will be modelled as a homogeneous material. The central lamina of a total thickness of 2t90° (t90° in Fig. 2) with the fibres oriented 90° with respect to the load direction (normal to the plane represented in the section appearing in Fig. 2), will be mesomechanically modelled as an isotropic material. Finally, the central part of the model will include a cell covering a fibre, where the damage will be considered. This damage, see Fig. 3, will consist in a first phase on a debonding (denoted in Fig. 3 by the semidebonding 𝜃 d ) between fibre and matrix, located nominally normal to the direction of the load although it is condemned in the model, reproducing the physics of the problem, to be contained along the interface between fibre and matrix. The properties of the material systems used in this paper appear in Tables 1 (carbon system) and 2 (glass system). The discussion about the possible strategies to deal with this problem (notice that the representation of the model associated with Fig. 2 is already giving indication about the final decision taken) requires to notice that the damage will be characterized by the Energy Release Rate
2. The problem: modelling strategies The problem under consideration is finally associated with the possibility of delaying the appearance of damage by using ultra thin plies to design and manufacture laminates. This possibility is based on the observed fact that a 90° ply increases its strength when its thickness tends to a minimum value. This value had been restricted in the past to around 120 μm due to technological reasons, which have been solved giving rise to the possibility of having laminas with a thickness of around 40 μm, having only around 5/6 fibres (of a diameter of 7.5 μm as in the case of carbon fibres) through the thickness of the ply. The onset and progression of damage in an isolated fibre of a 90° ply, under tension, are represented in Fig. 1, París et al. [8]. First, based on the interfacial stress state predicted, a minimum interface crack of about 20° will first appear, stage I in Fig. 1. This interface crack would grow till reaching a total debonding angle in the interval of 130–140°, stage II in Fig. 1, where contact zones at the ends of the interfacial crack appear. Finally, either the interfacial crack grows in a stable way under 2
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 3. Damaged cell in the 90° lamina with a debonding between fibre and matrix. Table 1 Carbon-epoxy system used. Material
Properties
0° ply (orthotropic)
E11 = 135 GPa, E22 = 8.75 GPa, E33 = 8.75 GPa, 𝜈 12 = 0.3, 𝜈 13 = 0.3, 𝜈 23 = 0.4, G12 = 4.75 GPa E(E22 = E33 ) = 8.75 GPa, 𝜈 = 0.4 E = 4.2 GPa, 𝜈 = 0.32 E (E22 ) = 15 GPa, 𝜈 = 0.2; R = 3.75 μm
90° ply (isotropic) Matrix (epoxy, isotropic) Fibre (carbon, isotropic)
Fig. 4. Detail of the model involving a debonding between fibre and matrix and an associated kinked crack.
representative enough to get accurate results. If the area where the debonding takes place is separated from the interface connection, the model for the second step needs to be changed, and a multifibre model for this second step will probably be required. On the contrary, if strategy ST2 is used, the main challenge will be to make compatible the required refinements of the discretization at the two scales involved in the model, avoiding having very similar boundary integral equations that may lead to a bad conditioned final system of equations. Obviously, a computer with more capacity is also required for employing this strategy, which was that finally employed. A similar approach was considered the best option for the analysis of local debondings in analysing composite components, Reinoso et al. [27], where both strategies were compared. Notice that in that case the same scale was required and only the refinement of one part of the component where the damage may appear called for the consideration of the two strategies already mentioned. In the present case, the difficulty is much more extreme as the refinement has to be performed at the microscale level, where the damage appears, combined with a mesoscale model of the rest of the problem. Once the first stage of the expected damage has been identified, it is reasonable to accept that the generation, from these debondings, of a mesocrack through the thickness of a 90° ply, will involve for the debonding crack the abandon of the interface between fibre and matrix, kinking into the matrix [8]. Thus, additionally to the cell already schematized in Fig. 3, a cell as that presented in Fig. 4, showing the presence of a debonding plus a kinked crack entering into the matrix must be studied too. Notice that this cell plays an additional role as that of the previous cell as it covers also the third stage of damage in a 90° lamina of a laminate, represented in Fig. 1. The used self-developed BEM code (which includes all the features necessary to use strategy ST2) is written in Fortran 90 and compiled in a Linux based system. This code is compiled to run in a vectorized manner i.e. able to use multiple processors within a High Performance Computing (HPC) System.
Table 2 Glass-epoxy system used. Material
Properties
Matrix (epoxy, isotropic) Fibre (glass, isotropic)
E = 2.79 GPa, 𝜈 = 0.33 E = 70.9 GPa, 𝜈 = 0.22; R = 7.5 μm
(ERR) of the debonding crack in accordance with interfacial Fracture Mechanics (IFM) [8], which in turns requires to calculate with high accuracy displacements and stresses along a certain circumferential interval (0.5° in this paper) at both sides of the crack tip, i.e. along a distance of 0.098 μm. These dimensions coexist in the model with the thickness of a lamina (oscillating in between 34 and 680 μm) and the length of the laminate (taken as 5000 μm). Notice also the non linear character of the problem, derived from the presence of potential contact zones that may appear between fibre and matrix. There are two strategies to deal with the problem under consideration from a modelling and numerical point of view, and will be designated as ST1 and ST2 in what follows. Strategy ST1: the problem is solved in two steps. In the first step the area of interest (of micromechanical dimensions) where the damage is going to take place is represented with a no very refined discretization, or even the micromechanical representation can be omitted and substituted by a homogeneous material with the same properties that the rest of the lamina. In a second step, the values of stresses or displacements are used as boundary conditions for the isolated cell where fibre, matrix and damage are now represented, Fig. 3. Strategy ST2: the problem is solved directly in one step, modelling in a very refined way the zone where interfacial fracture mechanics calculation will take place in conjunction with areas where refined discretizations are not required (as those as the extreme of laminates where the external action is applied). In the case of using ST1 it is obvious that the problems associated to step 1 and step 2 are reasonably easy to solve (with no too great differences in the size of the elements required on each step), although an iterative procedure, as the solution obtained in step 2 may not be compatible with that obtained in step 1, is required. With reference to this ST1 strategy, if the interface of connection between the first and second problem is not very well discretized, as the damage modelled in the second problem (the debonding between fibre and matrix) is very close to the line of separation between the two models, the results obtained in the second step may be contaminated by the poor discretization of the interface connection carried out in the first step, not being then
3. The BEM model The discretization of the multiscale model is a key factor to get accurate results in the area of interest. In order to evaluate the energy release rate associated to crack propagation, Irwin’s approach will be employed. Therefore, the area of interest is the neighbourhood of the tip of the crack representing the damage under consideration. As has been mentioned above, two different cases will be shown: the first one containing a debonding crack between fibre and matrix and the second one containing a matrix crack initiating from the debond crack and running normal to the applied load. 3
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 5. General view of the discretization carried out in the model.
As will be shown below, a virtual crack closure length of 3.27•10–2 μm (corresponding to a subtended angle in the interface of the fibre with the matrix of 0.5 °) will be employed in the calculation of the ERR which will be calculated, by using a numerical integration scheme in the range 6.54⋅10–6 μm < 𝜌<3.27⋅10–2 , with 𝜌 being the distance to the crack tip. Consequently, smaller elements in the vicinity of the crack tip should be at least one order of magnitude shorter than the lower limit of 𝜌, since linear elements are used and the presence of the singularity spoils the accuracy in the closest vicinity of the crack tip. The discretization of the rest of the model, as the results obtained associated to this zone are not of primary interest, must be focused on reducing as much as possible the number of nodes and on avoiding negative effects on the system of equations, as would be the presence of corners in some of the boundaries (those of the cells for instance) or the presence of boundaries near the damaged area where the calculations must be extremely accurate. Thus, regular rules used in BEM, París and Cañas [3], have guided the discretization away from the crack tip. Fig. 5 represents a scheme of the discretization carried out for a case having a debond crack. Four solids are involved in the model: S1 is the fibre, S2 is the matrix surrounding the fibre in the local cell, S3 is the 90° lamina and S4 is the 0° lamina. Each sector of the boundaries is identified as Sij indicating that is the j-th sector of solid Si . The additional information about using a uniform discretization (UD) or a progressive discretization (PD) appears in between brackets. In this particular case of a progressive discretization the length ratio in consecutive elements has been maintained between 1 and 1.2. Details of the mesh in the local scale will be shown in the next subsections in which the two different local cells will be described separately. The element length (in microns) is also indicated at each corner and in the zones in which a uniform discretization is used. Element length is identical in all elements sharing a corner node and conforming discretizations are employed in all contact zones or bonded interfaces. Notice that the ratio between the size of the model (5⋅103 μm) and the size of the smallest element (∼10–7 μm) is extremely high. Since the accurate evaluation of the boundary integrals when the collocation point is extremely far away from the smallest element is impossible, the use of a local and global scale is mandatory even if the material of solids S2 and S3 was identical. The use of the local model reduces the length ratio between the largest and the smallest elements of each solid allowing an accurate evaluation of the boundary integrals to be carried out. Finally, perfect bonding conditions are employed in the interface between S3 and S4 and between S3 and the local cells described in the following.
The most important information associated with the discretization carried out in the local cell is that regarding to the mesh employed in the neighbourhood of the crack tip. The debonding crack, being an interface crack, can be modelled in accordance with either an open model, Williams [28], which would involve interpenetrations near the crack tip of the faces of the crack, or in accordance with a contact model, Comninou [29], which would involve the appearance of a contact zone near the crack tip. In our case the BEM code employed allows the presence of a contact zone to be developed and consequently BEM solution directly satisfies the conditions of Comninou’s model, given that if interpenetrations tend to appear (using initially open conditions), the iterative procedure employed will adjust the incompatibilities found. Notice that the contact zone to avoid interpenetrations may be extremely small or relatively large depending on the fracture mode mixity, which varies with the length of the crack. However, this fact does not affect the discretization to be carried out, since the evaluation of the ERR has to be carried out at a scale along a distance of growth of physical relevance. Consequently if this distance scale is larger than the contact zone, the presence of extremely small contact zone does not affect the ERR associated to crack propagation at the scale of interest. Since friction between crack faces is neglected, energy release rate G can be evaluated according to Irwin’s Virtual Crack Closure Technique, and decomposed into fracture mode I contribution, GI , and fracture mode II contribution, GII : 𝐺 = 𝐺𝐼 + 𝐺𝐼𝐼 = +
1 Δ𝜃𝑑 ∫0
Δ𝜃𝑑
1 Δ𝜃𝑑 ∫0
Δ𝜃𝑑
[
( ) ( )] 𝜎𝑟𝑟 𝑅, 𝜃𝑑 + 𝜃 ′ Δ𝑢𝑟 𝑅, 𝜃𝑑 + Δ𝜃𝑑 − 𝜃 ′ 𝑑𝜃 ′
[ ( ) ( )] 𝜎𝑟𝜃 𝑅, 𝜃𝑑 + 𝜃 ′ + Δ𝑢𝜃 𝑅, 𝜃𝑑 + Δ𝜃𝑑 − 𝜃 ′ 𝑑𝜃 ′ ,
(1)
where (r,𝜃) is a local co-ordinate system centred with the fibre, 𝜎 rr and 𝜎 r𝜃 are, respectively, the normal and tangential tractions along the interface, while Δur and Δu𝜃 are, respectively, the opening and sliding between crack faces. The previous experience indicates, París et al. [8], that depending on the debonding angle, the interface crack fracture mode evolves from having an extremely small contact zone (for small debondings 𝜃 d < 60°) to showing a significant contact zone (for large debondings 𝜃 d < 60°). When the length of the contact zone is larger than Δa = RΔ𝜃 d (with R being the fibre radius), the scale in which the ERR is being evaluated, pure fracture mode II takes place. Fig. 7 shows a scheme of the discretization carried out (the same for all the cases studied) at both sides of the crack tip and Table 3 gives the values of the angle 𝛼 ij subtended by the element joining nodes i and j, 𝜆ij denoting the length of the element connecting nodes i and j. Since the accuracy of the solution in the closest vicinity of the crack tip is spoiled by the presence of the singularity, the solution of the closest nodes to the crack tip is not employed, and integration is carried out from node 5 up to node 16 (see Table 3), which means that the integrals in (1) are effectively carried out in the range 0.0001°<𝜃′<Δ𝜃 d = 0.5°.
3.1. Interface crack case Fig. 6 represents a scheme of the discretization carried out in the local cell for a case having only a debond crack. Notation employed follows the same rules described in the previous section. 4
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 6. Discretization of the local cell corresponding to a debonding case. Table 3 Subtended central angle and length of the elements in the neighbourhood of the crack tip.
𝛼 ij (0 ) 𝜆ij (μm)
𝛼 ij (0 ) 𝜆ij (μm)
1
2
3
4
5
6
7
8
6.25·10– 6 4.09·10– 7
6.25·10– 6 4.09·10– 7
1.25·10– 5 8.18·10– 7
2.5·10– 5 1.63·10– 6
5·10– 5 3.27·10– 6
2.04·10– 4 1.33·10– 5
4.16·10– 4 2.72·10– 5
8.5·10– 4 5.56·10– 5
9
10
11
12
13
14
15
16
1.73·10– 3 1.13·10– 4
3.53·10– 3 2.31·10– 4
7.21·10– 3 4.72·10– 4
1.47·10– 2 9.62·10– 4
3·10– 2 1.96·10– 3
6.12·10– 2 4·10– 3
1.25·10– 1 8.17·10– 3
2.55·10– 1 1.67·10– 2
open or partially open in all cases (different debondings). For the different debond angles studied, the mesh is identical to that shown in Fig. 6, corresponding to the 𝜃 d = 60° case, except for the fact that the mesh in the neighbourhood of the crack tip is shifted to its corresponding position. When the HPC system already described is used, the computation times are quite reasonable, being of the same order independent from the configuration studied (different length of the debonding). In the cases in which a contact zone is detected, only a few iterations are needed to obtain convergence and only a reduced number of contact equations has to be rebuilt in each iteration. To give an idea of the complexity of the problem, with reference to the case of a debonding of 40°, the total number of elements is 2070, leading to a total system of equations of 4728, with a minimum size of an element of 4,09·10–7 μm and a maximum size of element of 50 μm. When this model is run in the HPC system, which includes 8 cpus, an iteration takes around 40 s. Usually, to get a converged solution of the whole problem, only a reduced number of iterations (9 in the case described here) is needed, the total solution time being around 4–5 min. It is well known how sensitive is BEM, particularly in the presence of singularities, to the discretization carried out, and although the final proof of having a correct discretization will be obtained in the next section comparing with previous results, it has to be mentioned here that this study is not carried out with prediction purposes but with the aim to have a model to compare results associated to different configurations of a laminate. Thus, once the discretization employed is proved to work properly, it will be employed for all cases, then supporting the validity of the comparison for the problem under consideration.
Fig. 7. Discretization performed in the neighbourhood of the crack tip to analyse the first phase of damage (isolated debondings between fibre and matrix).
Note that the problem under consideration is nonlinear, as it implies the contact/separation of surfaces of fibre and matrix. This fact does not interfere with the intrinsic difficulty of solving of the problem, but it is significant during the process of definition of a suitable discretization. Thus, the crack was artificially closed along this process, adjusting the discretization, in order to get a smooth evolution of stresses in the neighbourhood of the crack tip. Once the discretization had been adjusted with the crack closed, it worked properly with faces of the crack
3.2. Kinked matrix crack case The local cell described in the previous subsection corresponds to the modelling of the first stage of damage, that is, the debonding 5
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 9. Predictions on G values obtained with the multiscale model (surrounded by matrix) and with the single fibre model.
Fig. 8. Discretization of the local cell corresponding to the kinked matrix crack case.
of a single fibre surrounded by matrix, both adimensionalized using a reference number G0 , París et al. [8]. Only the values of the total energy release rate are shown for the sake of clarity. It can be clearly observed that they are almost undistinguishable, proving then the validity of the multiscale model developed.
between fibre and matrix. This damage is assumed to continue, once the debonding crack stops, kinking and entering into the matrix following a direction normal to the applied load, what implies for the kinked crack to work under mode I, París et al. [8]. The model described in the previous sections has been adapted, as shown in Fig. 8, for the study of stage III of damage represented in Fig. 1. Although at the global scale, the models representing the two stages of damage have different sizes, the details of the mesh at the global scale of the model corresponding to the third stage identified in Fig. 1 are omitted for the sake of brevity, since the mesh is made following a similar scheme to that described in Fig. 5 for the model of the first and second stages. In the local scale, in order to easily assess stresses ahead of the crack tip, a fictitious interface along the predicted path of the kinked crack has been created, starting at the position of the original crack tip of the debonding and aligned normally to the applied load. Then, for a particular length of the assumed kinked crack, the closest neighbourhood of its crack tip is treated as described in the previous section. Contact conditions are employed between crack faces (from the symmetry plane to the crack tip) and perfect bonding conditions are applied along the fictitious interface where the crack has not extended. Other feasible option to study how the interface crack gets into the matrix (homogeneous material) without including a fictitious interface could be the use of the so-called Dual-BEM [30] or the use of a Symmetric Galerkin Scheme [31,32]. However, due to the fact that in the problem under study the direction of progression of the crack is a priori known, the method selected, in spite of the necessity of using a fictitious interface, was estimated to be the optimum solution. Moreover, the use of those schemes (dual BEM or SGBEM) would imply the use of hypersingular integral equations as well as their corresponding fundamental solutions.
4.2. Multiscale model results Once the validity of the model has been proved, one of the questions that attracts immediate attention is the shape of the hexagonal cell. Hexagonal cells are typically used when representing composites micromechanically as they are reproducible leading to a complete representation of the material at a plane normal to the fibres direction, allowing also in a simple way to reproduce the volume fibre ratio at the level of the cell. However, the fibres are in fact distributed in the material in a complete random way, and consequently any shape of the cell would be acceptable as will represent the actual situation of a certain fibre in the material. As the hexagonal cell has corners, this fact may introduce some alterations in the solution. Notice that in the case of having the corner, as in Section 4.1, located in a unique media, it would not produce any effect when using a FEM approach, but would do when using a BEM approach. However, having the corner of the cell in between dissimilar media (with matrix at one side and the homogeneous model of the 90° ply at the other side), it produces a nominal singular stress state associated to a multimaterial corner, Barroso et al. [33], BEM being much more sensitive than FEM to this situation. Then, an alternative circular cell is now going to be used in order to check if the presence of corners in the hexagonal cell produces significant alterations in the results. Figs. 10 (for a 30° debonding) and 11 (for a 50° debonding), represent the evolutions of the radial stress 𝜎 r ahead of the crack tip and the relative displacements Δur of the faces of the crack near the crack tip, values used to calculate G. It can be clearly observed the similarity of the results found for the two cells, hexagonal (HC) and circular (CC). As a general view of the situation, Fig. 12 represents the values of G predicted with the use of both shapes of the cell for the values of G and their counter parts GI and GII . When comparing the predictions of the results obtained with the two shapes of the cell, it has to be mentioned that they correspond obviously to different discretizations of the boundary of the cell, as uniform discretizations are performed for the circular cell, whereas a progressive discretization (involving more elements) is used in the case of hexagonal cells near the corners of the boundary of the cell. In any case the results are quite similar, having taken in Figs. 10 and 11 two representative situations involving one of them the case of greatest difference in the predictions. All of this supports the use of the BEM model initially represented in Figs. 2 and 3 to study the scale effect in composites by
4. Results 4.1. Validation of the model The key result, as will be used in further studies to discriminate on the role of the scale effect in the damage of composite laminates, is the energy release rate G. As the first question to be addressed is the ability of the multiscale model developed to get accurate results, those obtained with the present model will be compared with the ones obtained with a simpler model, a single fibre surrounded by matrix, París et al. [8]. To make the comparison effective, all the bodies, except for the fibre, appearing in the multiscale model will adopt the properties of the matrix, thus making the comparison advisable. Fig. 9 represents the values of G obtained with the multiscale model and those obtained with the model 6
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
Fig. 10. Radial stress and relative displacements in the neighbourhood of the crack tip for a debonding of 30° using a hexagonal cell (HC) and a circular cell (CC).
Fig. 13. Evolution of G, and their counterparts GI and GII , for different lengths of the kinked crack.
8, will now be considered. The case corresponding to a debonding of 65° has been selected, as this debonding is in the interval where, as justified in París et al. [8], the abandonment by the debonding crack of the interface fibre-matrix, kinking into the matrix, is more plausible. Different lengths of the kinked crack have been considered, Fig. 13 representing the evolution of the values of G of the kinked crack, and their counterparts GI and GII , as a function of its length for a carbon epoxy-system. First of all, the evolution obtained with the multiscale model is comparable to that obtained by París et al. [8] for a glass-epoxy system, using a single fibre model surrounded by matrix. Several things are noticeable in this evolution. First, coherently with the decision taken to build the model, negligible values of GII are obtained, leading to almost a pure mode I associated to the kinked crack. This is coherent with the reasoning that led in París et al. [8] to select the 60–70° range for the debonding as that more plausible for the kinking of the interface crack to occur. Additionally, it is also noticeable the increasing character of G (and obviously of GI ) with the increase of the length of the kinked crack. This leads to conclude that once the kinking of the debonding crack arises, the growth of the kinked crack will be unstable, based on the increasing character of G and on the unaltered character of the fracture mode of the kinked crack (mode I) along its path. Notice that a maximum kinked crack length of 2R has been considered in Fig. 13, although further extensions of the crack have been explored confirming the same behaviour. This study has required the use of a bigger than in previous cases. The same cell has been used for all cases considered involving a kinked crack of different length, which makes the results comparable among them. The aforementioned lack of stability of the growth of the kinked crack in the matrix detected by the present numerical model is in agreement with the experimental evidence presented in Correa et al. [34]. It is also important to have a look of the deformed shape of the model obtained involving the kinked crack, Fig. 14, where a magnification factor of 1.2 has been employed. In addition to the coherence of the deformed shape obtained (only the fibre and surrounding matrix are represented in order to concentrate in the area of interest), what can be observed is that independently of the length of the kinked crack chosen (the deformed shape represented corresponds to a length of the kinked crack of 2R (7.5 μm)), no interaction between the lips of the two cracks is detected. As soon as a kinked crack is modelled, no matter how short it is, the two faces of the parent debonding crack separate along the whole length, and thus the contact zone that existed prior to kinking for this particular interface crack size (65°) disappears. This result is also in agreement with París et al. [8].
Fig. 11. Radial stress and relative displacements in the neighbourhood of the crack tip for a debonding of 50° using a hexagonal cell (HC) and a circular cell (CC).
Fig. 12. Values of G, GI and GII for the whole range of debondings using hexagonal cells (HC) and circular cells (CC).
means of the energy release rate of the debonding crack, as the results are not affected by the shape of the cell. Once the presence of corners in the cell has been proved to have no influence in the results of interest, it justifies the use of a rectangular cell for the case of including the kinked crack in the model. Once the ability of the multiscale model to accurately predict the values of the parameter used to study the scale effect has been proved, the model generated to study the third phase of damage, Figs. 4 and 7
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
age appears at a micromechanical level (debonding in between fibre and matrix at a first stage) a multiscale model is required. BEM has been found to be the best option to deal with the problem under consideration, using one single model where both scales are modelled at the same time. The problem becomes very complicated, numerically speaking, due mainly to the scales used to represent the laminate, with a huge difference in the size of the elements involved in the model. Additionally, non-linearities at micro simulation level appear, derived from the consideration of the damage at the microscale, i.e. the aforementioned debondings between fibre and matrix and the kinking of these debond cracks entering into the matrix. The model has proved to be extremely accurate, based on the comparison with previous results and supported by experimental evidence. In addition, the ability of the model to deal with the first stages of the damage as a function of the variables that define the scale effect in composites has also been proved to be very robust and agile. Further studies will be carried out, París et al. [35], to check other predictions based on the multiscale numerical model developed in this paper. Additionally, the model is able to cover more complicated cells involving more than one fibre, as well as to study further stages of the damage involved in the scale effect.
Fig. 14. Undeformed and deformed shapes of the model (only the fibre and surrounding matrix are represented).
Acknowledgements This work was supported by the Spanish Ministry of Education, Culture and Sports under Grants MAT2016-80879-P and MAT2013-45069P. References [1] Parvizi A, Garret KV, Bailey JE. Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates. J Mater Sci 1978;13:195–201. [2] Flaggs DL, Kural MH. Experimental determination of the in situ transverse lamina strength in graphite/epoxy laminates. J Compos Mater 1982;16:103–16. [3] París F, Cañas J. Boundary element method. Fundamentals and applications. Oxford: OUP; 1997. [4] Sfantos GK, Aliabadi MH. Multi-scale boundary element modelling of material degradation and fracture. Comput Meth Appl Mech Eng 2007;196:1310–29. [5] Fernandes GR, Pituba JJC, de Souza Neto EA. FEM/BEM formulation for multi-scale analysis of Stretched Plates. Eng Anal Bound Elem 2015;54:47–59. [6] París F, Blázquez A, Mccartney LN, Mantič V. Characterization and evolution of matrix and interface related damage in [0/90]S laminates under tension. Part I: numerical predictions. Compos Sci Technol 2010;70:1168–75. [7] París F, Blázquez A, Mccartney LN, Barroso A. Characterization and evolution of matrix and interface related damage in [0/90]S laminates under tension. Part II: experimental evidence. Compos Sci Technol 2010;70:1176–83. [8] París F, Correa E, Mantič V. Kinking of transversal interface cracks between fiber and matrix. J Appl Mech 2007;74:703–16. [9] Correa E, Mantič V, París F. A micromechanical view of inter-fibre failure of composite materials under compression transverse to the fibres. Compos Sci Technol 2008;68:2010–21. [10] Correa E, París F, Mantič V. Numerical characterisation of the fibre-matrix interface crack growth in composites under transverse compression. Eng Fract Mech 2008;75:4085–103. [11] Correa E, París F, Mantič V. Effect of the presence of a secondary transverse load on the inter-fibre failure under tension. Eng Fract Mech 2013;103:174–89. [12] Correa E, París F, Mantič V. Effect of a secondary transverse load on the inter-fibre failure under compression. Compos Part B Eng 2014;65:57–68. [13] Mantič V, Távara L, Blázquez A, Graciani G, París F. A linear elastic-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads. Int J of Fract 2015;195:15–38. [14] Correa E, Gamstedt EK, París F, Mantič V. Effects of the presence of compression in transverse cyclic loading on fibre–matrix debonding in unidirectional composite plies. Compos Part A 2007;38:2260–9. [15] Correa E, Mantič V, París F. Effect of thermal residual stresses on matrix failure under transverse tension at micromechanical level: a numerical and experimental analysis. Compos Sci Technol 2011;71:622–9. [16] Correa E, París F, Mantič V. Effect of thermal residual stresses on the matrix failure under transverse compression at micromechanical level – a numerical and experimental study. Compos Part A Appl Sci Manuf 2012;43:87–94. [17] Sandino C, Correa E, París F. Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load. Eng Fract Mech 2016;168:58–75. [18] Muñoz-Reja M, Távara L, Mantič V, Cornetti P. Crack onset and propagation at fibre–matrix elastic interfaces under biaxial loading using finite fracture mechanics. Compos Part A – Appl Sci Manuf 2016;82:267–78.
Fig. 15. Values of G, GI and GII for different values of the t0 ° thickness.
4.3. Incidence of the 0° ply thickness in the scale effect Once the ability of the multiscale model developed to work satisfactorily has been proved, a first application of the model to the problem to which the model has been developed for, that is the scale effect in composites, is now going to be developed. The scale effect is sometimes referred to by means of the thickness of the 90° ply, Parvizzi et al. [1], or by a number n representing the proportion of 90° plies with regards to 0° plies in a [0,90n]S laminate, Flaggs and Kural [2]. Three different thicknesses, t0° = 85 μm, 170 μm and 340 μm, which correspond to n = 0.5, 1 and 2, will be considered in this study. The values of G, GI and GII for the three thicknesses mentioned and for the whole interval of debondings considered are represented in Fig. 15, with reference to the multiscale model considering only the first stage of damage. The independence of the values of G, GI and GII observed for the three thicknesses modelled, clearly indicates that the value of n does not play any role in the development of the first phase of damage (stages I and II in Fig. 1), that involving isolated debondings in between fibre and matrix.
5. Conclusions A multiscale model to deal with the onset and propagation of damage in composite laminates, basically those under the denomination [0,90n]S , has been developed. The aim of the development of the model is to have a tool to analyse the scale effect in these laminates, i.e. the onset of damage at different load values when the thickness of the 90° lamina or the value of n changes in the laminate. As the scale effect is controlled by one parameter at the scale of the laminate and the dam8
M.L. Velasco et al.
Engineering Analysis with Boundary Elements 93 (2018) 1–9
[19] Távara L, Mantič V, Graciani G, París F. Modelling interfacial debonds in unidirectional fibre-reinforced composites under biaxial transverse loads. Compos Struct 2016;136:305–12. [20] Távara L, Mantič V. Crack paths formed by multiple debonds in LFRP composites. Mech Res Commun 2017;84:148–54. [21] Graciani E, Mantič V, París F, Blázquez A. Weak formulation of axi-symmetric frictionless contact problems with boundary elements. Application to interface cracks. Comput Struct 2005;83:836–55. [22] París F, Correa E, Cañas J. Micromechanical view of failure of the matrix in fibrous composite materials. Compos Sci Tech 2003;63:1041–52. [23] Panagiotopoulos CG, Mantič V. Symmetric and nonsymmetric debonds at fiber-matrix interface under transverse loads. An application of energetic approaches using collocation BEM. Anales de Mecánica de la Fractura 2013;30:125–30. [24] García I, Mantič V, Graciani E. Debonding at the fibre–matrix interface under remote transverse tension. One debond or two symmetric debonds? Eur J Mech A/Solids 2015;53:75–88. [25] Muñoz-Reja M, Távara L, Mantič V. Symmetrical or non-symmetrical debonds at fiber–matrix interfaces: a study by bem and finite fracture mechanics on elastic interfaces. J. Multiscale Model 2017;08:1740008. [26] Mantič V, Távara L, Blázquez A, Graciani E, París F. A linear elasti-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads. Int J Fract 2015;195:15–38.
[27] Reinoso J, Blázquez A, Estefani A, París F, Cañas J. A composite runout specimen subjected to tension-compression loading conditions: experimental and global-local finite element analysis. Compos Struct 2013;101:274–89. [28] Williams ML. The stress around a fault or crack in dissimilar media. Bull Seismol Soc Am 1959;49:199 -24. [29] Comninou M. The interface crack. J Appl Mech 1977;44:631–6. [30] Aliabadi MH, Saleh AL. Fracture mechanics analysis of cracking in plain and reinforced concrete using boundary element method. Eng Fract Mech 1998;93:115–44. [31] Sutradhar A, Paulino G, Gray L. Symmetric Galerkin boundary element method. Berlin: Springer; 2008. [32] Távara L, Mantič V, Salvadori A, Gray LJ, París F. Cohesive-zone-model formulation and implementation using the symmetric Galerkin boundary element method for homogeneous solids. Comput Mech 2013;51:535–51. [33] Barroso A, Mantič V, París F. Singularity analysis of anisotropic multimaterial corners. Int J of Fract 2003;119:1–23. [34] Correa E, Valverde MI, Velasco ML. París, F. Microscopical observations of inter-fibre failure under tension. Compos Sci Technol 2018;155:213–20. [35] París F, Velasco ML, Correa E. Micromechanical study of the influence of scale effect in the first stage of damage in composites. Compos Sci Technol 2018;160:1–8.
9