Energetically-consistent multiscale analysis of fracture in composites materials

Two distinct length scale transition methodologies are developed to establish effective traction-separation relations for fracture in composite materials within a hierarchical multiscale framework. The two methodologies, one kinetics-based and the other kinematics-based, specify effective fracture properties that satisfy a surface-based Hill-Mandel consistency condition. Correspondingly, the total amount of energy dissipated is the same whether a crack is described in detail with micro quantities or in terms of an effective macroscopic crack. Though both methods guarantee consistency in terms of energy rates across length scales, they provide in general distinct effective traction-separation relations. Several representative samples of fiber reinforced composites are analyzed numerically, including the formation and propagation of cracks at mid-ply locations as well as (idealized) ply interfaces. Through post-processing of the microscale results, it is shown that the kinematics-based averaging method provides a macroscopic traction that is prone to rapid fluctuations while the kinetics-based averaging method shows a more smooth response but with openings that can deviate from the surface average of the microscale openings. The two methods are also compared with a previously-proposed scale transition methodology, which is a hybrid method that only satisfies the Hill-Mandel condition approximately. The suitability of the three methods is discussed in light of the results obtained from the simulations.


Introduction
With an ever increasing demand for more efficient lightweight composite materials in the transportation, infrastructure and energy conversion sectors, new types of composite materials are continuously being designed and tested. However, adopting a new type of composites is often hindered by development costs associated to expensive material and structural tests. Designers and certification authorities require a high degree of confidence on the performance of new materials, in particular their actual capacity to safely carry loads, hence a robust fracture theory is a critical aspect of material development. In this context, advanced simulation methods provide a powerful approach to reduce experimental testing costs and shorten design cycle times (i.e., a virtual testing environment (Cox and Yang, 2006;Lopes et al., 2016)). To achieve this goal at the level of coupon testing (material performance), a reliable model is required in order to predict the overall (effective) fracture properties of new and existing composite materials based on the elementary properties of its constituents. Multiscale simulations enable modeling and analysis of fracture processes at the microscale, thus delivering fundamental information about crack processes that purely phenomenological models cannot capture. In particular, the mechanical behavior of a composite depends on the properties of its constituents (e.g., fibers, matrix and fiber-matrix interfaces) as well as its geometrical arrangements (e.g., fiber volume fraction, ply orientation, stacking sequences). By implementing microscale failure modeling at the level of individual constituents, multiscale simulations capture fundamental fracture processes at the microscopic level Melro et al., 2013;Talreja, 2014). This is crucial for composite materials in which fracture entails highly complex nonlinear processes developing at the micro level. Furthermore, thorough understanding of fracture processes is essential to realize the full potential of advanced materials in structural applications. For instance, micro-structural modeling enables researchers and engineers to virtually tailor composite materials, paving the road for microstructural modification and optimization (Okereke et al., 2014). In addition, multiscale simulations allow the incorporation of microscale features that are responsible for statistical variations in terms of failure loads, thus providing a direct link between microscale defects and the quantification of uncertainty in the fracture behavior of composites (Vajari, 2015;Maragoni et al., 2016;Turteltaub and de Jong, 2019). A successful implementation of this framework naturally leads to more accurate simulations and more efficient lightweight designs with the associated benefits in terms of reducing the overall resources required for a given application.
One important issue in the context of multiscale simulations is the notion of scale transition relations, in particular for the so-called hierarchical methods. In a hierarchical approach, the small and large scales are separated, meaning that the kinematical and kinetic fields may be defined separately per scale. Since a classical continuum theory has no intrinsic length scale, the balance principles apply separately in each scale and, consequently, for each set of field quantities. There are some well-established techniques to relate two continuum-based theories, but in general a scale transition relation is meant to preserve some notion of consistency across scales. Several open issues have been identified regarding scale transitions from macroscopic to microscopic in terms of material response (van der Meer, 2016;Chevalier et al., 2019) as well as from microscopic to macroscopic in terms of convergence Gitman et al. (2007); Phu Nguyen et al. (2010); Bosco et al. (2014); Svenning et al. (2016); Goldmann et al. (2018). To guarantee consistency between quantities that are meant to represent the same physical phenomenon across length scales, one has to verify (or otherwise enforce) that the same values are obtained whether one works with effective macroscopic fields or detailed microscopic fields. In particular, the Hill-Mandel condition (in rate form) stipulates that the rate of dissipation of energy should be the same whether one homogenizes the microscale dissipation or computes it directly from the homogenized fields.
The present work focuses on consistency of fracture from the microscale to the macroscale. Satisfaction of a surface-based Hill-Mandel scale transition for fracture in general cannot be guaranteed a priori with a classical multiscale approach in which the effective properties are defined as volume or surface averages of the microscopic quantities. Some of the aforementioned references have dealt with the issue of homogenization of fracture behavior. One possible remedy to this situation was also proposed in Turteltaub et al. (2018), where the effective cohesive traction on a crack was defined as a linear combination of surface and volume-based stress quantities, with an additional parameter to approximately satisfy the scale transition condition. In the present work, two alternative scale transition approaches are considered, namely a kinematics-based approach that relies on an alternative definition of the traction on a cracked surface and a kinetics-based approach that is based on an alternative definition of the effective crack-opening rate. In contrast with the method proposed in Turteltaub et al. (2018), these two alternative scale transition methods exactly satisfy the Hill-Mandel condition for fracture, albeit each method provides a distinct effective macroscopic response (i.e., a distinct effective traction-separation relation). Correspondingly, the objectives of the present work are as follows: (i) to propose and develop two scale transition methods for effective fracture properties, (ii) to compare quantitatively and qualitatively the effective macroscopic responses of these methods as well as a previously-proposed method, (iii) to introduce a methodology to homogenize intersecting cracks, a situation that has hitherto been avoided in previous studies and (iv) to show that this methodology can be used to identify a representative surface element for fracture. With regards to the case of intersecting cracks, it is worth mentioning that this situation is relevant for anisotropic composites under multiaxial loading whereby distinct fracture mechanisms generate intersecting cracks, such as delaminations and transverse ply cracking at ply interfaces in fiber reinforced composites.
The work is organized as follows: the microscale problem is formulated in Sec. 2 and the corresponding requirements for a length scale transition towards a macroscale problem are presented in Sec. 3. The kinetics-based and kinematics-based effective quantities for the macroscale formulation are defined in Sec. 4. A verification of the satisfaction of the scale transition requirements under various loading cases is shown in Sec. 5 together with a comparative analysis of the predictions of each scale transition methodology, including a previously-proposed methodology. This section also contains a demonstration that the methodology can deal with the case of intersecting cracks. Subsequently, the new scale-transition approaches are applied to a multiscale convergence analysis in Sec. 6, where it is shown through examples that the kineticsbased methodology outperforms the kinematics-based approach in terms of establishing effective properties for hierarchical multiscale analysis. Finally, some concluding remarks are given in Sec. 7, including the limits of the proposed methodology.

Microstructural formulation
The smallest length scale of analysis is taken as an aggregate of distinct phases such as fibers embedded in a matrix as shown in Fig. 1. In view of the computational effort required to carry out a large number of multiscale simulations, a plane-strain two-dimensional approach is chosen. This choice provides a computationally-tractable environment to develop and study energetically-consistent averaging methods and can be also applied to the three-dimensional case, albeit at a larger computational cost. For Carbon Fiber Reinforced Polymers (CFRP), a typical two-dimensional Microscopic Volume Element (MVE) is chosen to represent a cross-section with a stacking of uni-directional plies arranged in perpendicular directions. The methodology presented here, however, is not limited to CFRP and may be applied to any composite material.
A typical microstructural domain Ω with boundary ∂Ω is illustrated in Fig. 2. Each individual edge of the domain is denoted by Ω i , with outward normal vector n i ; i 2 ½1; 4� and the global orthonormal basis used is given by e 1 and e 2 . The collection of all cracked surfaces is denoted by Γ and it typically consists of one or more main cracks Γ I , I ¼ 1; 2; …, that may intersect and/or have bifurcated branches. The crack normal is denoted as m with m ¼ m þ pointing to the þ side and m À ¼ À m þ pointing towards the À side. Fig. 2 also displays the periodicity of the cracks, which continue over the edges of the domain due to the application of periodic boundary conditions. While other boundary conditions can be applied Geers et al. (2017);Peri� c et al. (2011), in this derivation periodic boundary conditions (PBCs) are chosen for simplicity.
The fracture process in the MVE is formulated as a quasi-static boundary value problem with equilibrium satisfied at each time t in the (uncracked) bulk material ΩnΓ and traction continuity imposed across the crack surface Γ, i.e., div σðx; tÞ ¼ 0 with σ being the stress tensor, t the traction vector, div the divergence operator and x a point in the reference configuration. As mentioned above, periodic boundary conditions are enforced in the MVE surface ∂Ω for displacements u and anti-periodic conditions for the tractions t, namely uðx þ l 1 e 1 ; tÞ À uðx; tÞ ¼ l 1 εðtÞe 1 tðx þ l 1 e 1 ; tÞ ¼ À tðx; tÞ uðx þ l 2 e 2 ; tÞ À uðx; tÞ ¼ l 2 εðtÞe 2 tðx þ l 2 e 2 ; tÞ ¼ À tðx; tÞ u � ðx � þ l 1 e 1 ; tÞ À u � ðx � ; tÞ ¼ l 1 εðtÞe 1 t � ðx � þ l 1 e 1 ; tÞ ¼ À t � ðx � ; tÞ u � ðx � þ l 2 e 2 ; tÞ À u � ðx � ; tÞ ¼ l 2 εðtÞe 2 t � ðx � þ l 2 e 2 ; tÞ ¼ À t � ðx � ; tÞ with ε being a given macroscopic strain tensor that drives the defor-mation. The quasi-static problem is complemented with initial conditions in a (typically) uncracked material. It is noted that the crack Γ is not known a priori but is in fact an outcome of the simulation. Further, when the crack crosses the MVE boundary, the conditions (3c) and (3d) are enforced. This actually reflects the periodicity of the crack discontinuities. Indeed, denoting the crack opening as EuF with and combining this with (3c) and (3d), yields the following periodic conditions For simplicity, the current formulation is developed for small strains whereby the relevant stretch part (micro-scale strain field) at points away from the crack is given by where r denotes the gradient operator and T the transpose. The constitutive behavior of the composite constituents is assumed to be governed by linear elastic relations up to fracture, i.e., σ ¼ Cε with C being the fourth-order elastic stiffness tensor of the corresponding phase (e.g., fiber or matrix). Moreover, the fracture behavior in the MVE is modeled with a micro-scale traction-separation relation (cohesive relation) expressed as where the traction t on the crack surface Γ typically degrades from the initiation value (fracture strength) to zero for a fully-opened crack. This is typically described by a cohesive relation f coh that depends on the crack surface opening EuF, damage variable(s) κ and the normal vector m ¼ m þ to account for an opening mode and/or anisotropic fracture. Distinct traction-separation relations are used for the phases in the MVE (matrix and fibers) as well as interfaces (e.g., separate relations for sizings in fiber/matrix interfaces representing distinct bonding chemistries). It is worth pointing out that the current modeling approach is limited to brittle fracture since no other inelastic behavior is incorporated in the distinct phases of a composite. The goal of the multiscale analysis is to link the individual tractions separation relations, assumed to be known, to the overall (macroscopic) fracture properties.

Overview of requirements
Since classical continuum mechanics has no intrinsic length scale, all  , denoted as Ω, with main intersecting cracks (Γ 1 and Γ 2 ) and secondary cracks (branches, isolated segments). Each crack surface has a positive and a negative side, indicated by superscripts þ and À , and associated quantities, namely an outward normal m, displacement u and traction t (not shown for clarity).
aspects of the theory should formally be the same at all length scales, except for possible differences in constitutive models. In particular, the same balance principles should be satisfied at all length scales, albeit with the corresponding micro or macro quantities and constitutive models. Thus, in a hierarchical multiscale framework, it is critical to consistently connect two distinct continuum descriptions at two distinct length scales in a way in which there are no contradictions or inconsistencies in the balance principles. Using a cohesive zone approach at both the small and large scales, the length scale transition for power requires that the rate of work dissipated across scales should be equivalent. This requirement is the Hill-Mandel condition applied to the fracture process. Thus, the energy dissipated inside a representative volume element (RVE) and ascribed to the cracking process should coincide with the energy dissipated obtained from a macroscopic continuum point on an equivalent macroscopic crack at the same physical location as the microscale volume element. An important aspect in this procedure is a suitable separation of the bulk material surrounding the crack and the crack itself, with the purpose of identifying a representative surface element (RSE) inside the microscopic volume element. Two alternative approaches to satisfy the scale transition requirement on a representative surface element are presented in the sequel, one that is kinematics-based and the other that is kinetics-based. Each approach is energetically-consistent, but provides a distinct tractionseparation relation.

Strain, stress and power relations
Consider a microscopic volume element Ω with boundary ∂Ω and denote as Γ the collection of all cracked surfaces as illustrated in Fig. 2.
The externally applied (macroscopic) strain ε acting on the volume element is defined as where jΩj denotes the volume of the region Ω, the subscript "sym" indicates the symmetric part of the quantity within square brackets and � denotes the tensor product. Observe that the applied strain ε is computed based only on the displacements u on the external boundary ∂Ω. Assuming periodic boundary conditions on the boundary ∂Ω, the macroscopic strain ε applied to the volume element can be decomposed as follows: where the volume averaged strain ε b is computed as and the fracture strain ε f is defined as where m is a unit vector normal to the crack surface (see Fig. 2 for notation).
The externally applied (macroscopic) stress tensor σ acting on the microscopic volume element is defined as with t representing the traction acting on the external boundary of the MVE and x being a position vector in the reference state corresponding to a microscale material point. Similar to the definition of the applied strain, the applied stress tensor only depends on the traction applied on the external surface ∂Ω. Upon application of the balance of linear momentum for a quasi-static process without body forces and taking into account that the traction on the crack surface is continuous (continuous across the crack surface), it can be shown that the externally applied macroscopic stress σ coincides with the volume averaged stress tensor 〈σ〉 Ω , i.e., where the volume averaged stress tensor is defined as The externally applied power density P (per unit volume) done on the MVE is defined as where _ u is the time derivative of the displacement vector (velocity). For a quasi-static process with periodic boundary conditions and without body forces, using the equation of equilibrium and the divergence theorem, it can be shown that the external power per unit volume can be expressed as where the stress power per unit volume P b (also referred to as the bulk stress power density) is given as and the fracture power per unit volume P f is defined as For plane strain or plane stress formulations, the power densities are expressed per unit area and unit depth.

Hill-Mandel condition
At a macroscopic level, the stress power density P M at a continuum point is, by definition, given as where σ M is the macroscopic stress tensor and _ ε M is the macroscopic strain rate. As discussed in Sec. 3.1, the Hill-Mandel consistency condition across length scales indicates that for a microscopic volume element to actually represent a macroscopic continuum point, then the total microscopic rate of work on the volume element must be equal to the local rate of work at the macroscale. In accordance with the theory of multiscale analysis, the macroscopic fields at a continuum point are related to the applied quantities on the boundary ∂Ω of the corresponding microscopic volume element by definition as The consistency condition for the whole MVE ("total" Hill-Mandel condition) is therefore that one has to verify that P M ¼ P : As indicated in Turteltaub et al. (2018), the Hill-Mandel condition for the whole MVE as given in (22), can be satisfied a priori by imposing periodic boundary conditions on the boundary of the domain. Consequently, this condition is automatically satisfied, even in the presence of cracks and regardless of whether the MVE is representative or not. However, in situations where a localized crack appears in the MVE, the traditional approach to multiscale analysis based on successive MVEs of increasing volume generally fails to converge to a RVE. This is often due to the fact that the ratio between the crack surface and the MVE volume is not constant as the MVE volume is increased. This problem may be partially solved if the MVE is varied in volume only by increasing its dimension along the crack in order to keep the surface to volume ratio constant, but this approach requires an a priori knowledge of the crack orientation hence it is typically limited to the analysis of fracture along a pre-determined path such as a weak interface between two distinct materials. In order to analyze a more general case in which the crack orientation is not known a priori, the scale transition for the cracking process requires identifying and separating the deformation mechanisms occurring in the bulk from the process localized on the microscopic cracked surfaces, which is then ascribed to fracture. This approach results in two separate scale transition conditions, namely one for the bulk and one for the crack. This segregation of scale transitions is used to extract information corresponding to the cracking process from the microscale. Since the fracture scale transition is mainly performed averaging along cracked surfaces, it converges to a representative surface element (RSE) and it does not depend on the ratio between the crack surface and the MVE volume.
A key aspect in the scale transition for the fracture process is the format in which the effective traction-separation relation is formulated. In particular, since traction-separation relations are used in the present work at both length scales, the macroscopic (or effective) description of fracture is also expressed in terms of a macroscopic cohesive traction t f acting on an equivalent macrocrack surface and a macroscopic crackopening vector EuF f . The macroscale description of an (equivalent) crack can therefore also be expressed in the same format as the microscale relations (8), i.e., where κ f represents a vector of (internal) damage or history variables and m f is a unit vector normal to the macroscopic crack. The macroscopic crack is associated to an infinitesimal area (or segment per unit depth in plane formulations) that is perpendicular to the normal vector. Denote as � � Γ f � � the length of a straight macroscopic crack infinitesimal segment Γ f for plane formulations that is meant to be the continuum equivalent representation of microcracks in a Representative Surface Element. Multiplying the fracture power per unit volume P f given in (19) by the MVE volume jΩj, and using the terminology introduced above for an effective macroscopic crack, the Hill-Mandel condition for fracture can be expressed as As discussed in Turteltaub et al. (2018); Turteltaub and de Jong (2019), the macroscopic cohesive traction (vector) t f and the macroscopic crack-opening rate vector E _ uF f are not necessarily obtained directly from the volume-averaged stress tensor 〈σ〉 Ω and the fracture strain rate tensor _ ε f . Some alternative definitions for the effective quantities are given in the next section.

Effective quantities
In order to satisfy the surface-based Hill-Mandel scale transition condition for fracture as given in (24), two distinct approaches are proposed in this section, namely a kinematics-based method where the crack opening rate is obtained from an average and the traction is adjusted in accordance with the scale transition and a kinetics-based methods where the traction is obtained from an average and the crack opening rate is adjusted from the scale transition. In addition, a version of a previously-proposed method (as discussed in Turteltaub et al. (2018)) is also included in this section in order to compare the predictions from the two new methods with the existing one.

Effective crack length and crack-based quantities
For subsequent use in the methods presented in this section, the effective crack-length and the crack-averaged traction and projected opening rate are defined in this section. The effective (macroscopic) length of a periodic crack � � Γ f � � is computed using the geometrical interpretation proposed in Turteltaub et al. (2018)), which for completeness is summarized here. The approach is to identify the orientation and number of periodic crossings of a periodic crack with normal unit vector m f in a two-dimensional l 1 � l 2 MVE domain aligned with normal unit vectors n 1 and n 2 , as illustrated in Fig. 2. The orientation of the effective crack is determined from the unit vector normal m f which is defined here as the crack-averaged normal vector, i.e., The effective crack length is defined as where the lengths and the nominal number of periodic crossings, expressed as a real number, is defined as In (26), the quantity r max is a cut-off value for handling near vertical or near horizontal periodic cracks for which the effective length could potentially predict an artificially large value instead of l 1 or l 2 . Similarly, in case of complex crack patterns with a large number of branches, it is convenient to compare the predicted effective crack length with an alternative approach based on the length of a vector connecting directly the entry and exit points of a periodic crack. The magnitude of that vector can be interpreted as the effective length of the main crack without the crack branches and hence equal to a straight segment crossing the computational sample. This value can be used as an alternative definition of the effective crack length; in principle the definition of the effective crack length is somewhat arbitrary since, as shown below, what is relevant is only the product between the effective length and the effective crack opening rate that is in turn consistently defined based on the chosen effective crack length.
The crack-averaged traction t f Γ is computed as while the nominal traction t f Ω associated with the volume-averaged stress tensor is defined as Similarly, define the nominal crack opening rate E _ uF f Γ associated with the fracture strain tensor as where the expression (12) was used to obtain the final expression of the nominal opening rate.
Observe that local crack opening rates E _ uF on crack segments that are perpendicular to the average crack normal vector m f have no net contribution to the value of E _ uF f Γ , which can also be interpreted and referred to as a projected crack-averaged opening rate. With the aforementioned definitions, the distinct scale transition approaches are introduced next.

Kinematics-based averaging method
In the kinematics-based averaging method, the effective crack opening rate E _ uF f is chosen as the projected crack-averaged opening rate (30) while the effective cohesive traction t f is computed from the Hill-Mandel condition (24), i.e., and Observe that the scalar β is a function of time during the cracking process such that Hill-Mandel condition (24) is automatically satisfied

Kinetics-based averaging method
In the kinetics-based averaging method, the effective cohesive traction t f is chosen as the crack-averaged traction t f Γ given in (28) while the effective crack opening rate E _ uF f is computed from the Hill-Mandel condition (24), i.e., and As in the previous method, it is noted that the Hill-Mandel condition (24) is automatically satisfied for the pair fE _ uF

Hybrid method: kinematics-based with approximate Hill-Mandel condition method
For comparison purposes, a modified version of a previouslyproposed scale transition method is included in this section. The method is essentially similar to the one presented in Turteltaub et al. (2018); Turteltaub and de Jong (2019) with the exception that the effective crack-opening rate used is the projected crack opening rate E _ uF f Γ (instead of the unprojected rate). This modification facilitates a direct comparison with the aforementioned methods that also use the projected crack opening rate. For convenience, this method will be referred to as the "hybrid method" or as the "approximate Hill-Mandel method" in the sense that it combines features of the kinetics and kinematics-based methods while enforcing the Hill-Mandel condition only approximately. Indeed, similar to the kinematics-based approach presented above, the effective crack opening rate is chosen as However, in this scale transition approach, the effective traction is not directly obtained from the Hill-Mandel condition (24) but, instead, it is written in terms of linear combination of the volume and the surfacebased tractions, namely as where the quantity α can be computed such that the pair fE _ uF f Γ ; t f α g approximately satisfies the Hill-Mandel condition (24). One way to achieve this is to substitute (36) and (35) in (24) and integrate throughout the cracking process from the initial state at t ¼ 0 (typically uncracked) until the fully-cracked state at t ¼ T. This procedure yields the following (constant) value of the scalar α: In some of the simulations carried out in this work, the crack and volume-averaged tractions were similar throughout the cracking process in which case the value of α can become prone to numerical inaccuracies. However, even in that case the proper traction can be recovered from (36), namely t f α � t f Γ � t f Ω , as long as a bounded value of α is obtained, which in practice can be a fixed cut-off value (e.g., α ¼ 0 or α ¼ 1).

Verification of scale transition relations and comparative analysis
In order to verify that the proposed methods and their numerical implementation actually satisfy the scale transition relation, a series of basic simulations are conducted using microscopic volume elements representing typical composite materials. Simultaneously, the predicted macroscopic traction-separation relations from the distinct methods are compared under different loading cases.

Implementation: microstructural samples, material properties and loading cases
The microstructure chosen for verification consists of plies of a unidirectional fiber-reinforced composite in a generic ½0=90� n stacking sequence of n plies with fibers oriented at angles 0 � and 90 � corresponding, respectively, to the e 3 and e 1 directions as indicated in Fig. 2. Due to computational limitations, artificially thin plies are assumed, which only accommodates a few fibers in the thickness direction (� 4 for the simulations considered here). This arrangement is not representative of composites currently used in practice (in which the plies are significantly more thick) but is adopted in the present work in order to keep the overall computational cost within a manageable range based on currently available hardware. Despite this limitation, the benefit of this microscopic volume element is that it is able to capture some of the main physical features of fracture in actual composites (i.e., matrix cracking, fiber cracking, fiber debonding, fiber pull-out and ply delamination).
The results reported in this section pertain to microstructural volume elements representing 75 μm � 75 μm cross-sections with ½0� and ½90� plies. Within each ply, fibers are randomly distributed in each distinct realization, while keeping the fiber volume fraction (approximately) fixed. Samples of other sizes and realizations are analyzed in subsequent sections; in the present section only one typical sample is shown for illustration purposes. The elastic and fracture properties of the composite system used in the simulations are summarized in Tables 1 and 2, respectively. These properties are representative of a commonly used combination in aerospace applications (IM7 carbon fibers and 8552 epoxy matrix), although the method is general and is not limited to this material choice or type of composite. The elastic and some of the fracture properties may be readily obtained from the manufacturer's published data. However, it is important to emphasize that some fracture data are not (publicly) available. Most of the values indicated in Table 2 are estimates and are used in the present study only for computational purposes but should not be used for design purposes. Experimental testing is required to establish accurate fracture data for the constituents (fibers, matrix and fiber-matrix interface), but this falls outside of the scope of the present work. Several representative loading cases are considered to verify and compare the proposed methodologies. The loading cases are defined in terms of the macroscopic strain tensor ε that drives the deformation inside the MVE. In particular, three loading cases are analyzed as shown in Table 3, namely (i) laterally-constrained axial extension in the ½90� fiber direction, (ii) a mixed loading case (equibiaxial extension combined with pure shear) and (iii) equibiaxial extension. The Cartesian components of the macroscopic strain tensor are referred to the tensor bases constructed from the underlying vector basis fe 1 ; e 2 g, with the 1direction aligned with the ½90� fibers in a ½0=90� n stacking sequence as illustrated in Fig. 2 (as indicated above, the ½0� direction is chosen as the global out-of-plane direction e 3 ¼ e 1 � e 2 ). The magnitude of the strain tensor, which depends on the parameter λ > 0, is chosen such that it is large enough to produce complete fracture of the specimen.
The microstructural volume elements are generated and meshed using the open source package Gmsh (Geuzaine and Remacle, 2009). The meshes contain 3-noded, linear, plane strain elements for the bulk deformation and 4-noded two-dimensional cohesive elements, embedded on each edge of the bulk elements, to simulate the fracture process that may initiate at any location where the local fracture criterion of the cohesive element is satisfied. The cohesive stiffness is chosen relatively high to minimize the effect of the artificial compliance that appears from embedding the cohesive elements. The quasi-static problems as described in (1)-(3) are solved using the FEA software Abaqus with implicit time integration (Abaqus Standard version 6.14). Preliminary mesh refinement analyses were conducted to find mesh sizes that provide converged solutions to within a given tolerance, resulting in a characteristic mesh size of 1 μm for most simulations. A numerical viscosity parameter in the cohesive elements is also usually required to find converged solutions. Ideally, this parameter should be as low as possible since it does not represent physical dissipation, but it is typically necessary to regularize the problem. The solutions are monitored in terms of the contribution of the viscous term to the overall energy dissipation. Numerically-converged solutions that contain a large amount of (numerical) viscous dissipation are not considered as converged from the point of view of multiscale analysis. In some cases this only affects parts of the solution (e.g., towards the end of the fracture process when a large number of cohesive elements are simultaneously active). The converged solutions are post-processed with several python-based scripts to detect cracks groups, identify their connectivity and compute the effective traction-separation quantities according to the distinct scale transition approaches.

Laterally-constrained uniaxial extension in the fiber direction
The first loading case corresponds to extending the specimen in one of the fiber directions while preventing contraction in the perpendicular direction (i.e., laterally-constrained uniaxial extension as indicated in Table 3). The corresponding cracking process and final crack pattern is shown in Fig. 3. Although the details of the crack pattern vary from sample to sample, similar samples show typically a main periodic crack (indicated by a thick line in Fig. 3) that runs through the matrix and the matrix-fiber interfaces on the ½0� layer. Eventually the fibers in the ½90� layer break indicating the complete failure of the sample. Throughout the fracture process, crack bifurcation occurs at multiple sites as well as the formation of isolated crack segments that, after an initial growth, get arrested as the stress decreases due to the formation of a main crack. It is nevertheless important to consider all branches in a multiscale analysis since these isolated crack segments may have a non-negligible contribution to the overall energy dissipation.
In order verify the averaging methods indicated in the previous section, the detailed data of the simulation is post-processed in three different ways, namely using the kinematics-based method, the kineticsbased method and the kinematics-based method with approximate satisfaction of the Hill-Mandel condition (i.e., the hybrid method). To this end, a crack detection algorithm was developed to collect all failed cohesive elements and group them in larger sets (crack segments), which then in turn are collected into connected cracks or isolated segments. Subsequently, the distinct averaging methods are applied to the same set of elements. Distinct power terms for one illustrative simulation are shown in Fig. 4. The results shown in the figure include both surfacebased and volume-based quantities, which are all normalized with respect to the MVE volume jΩj (area per unit depth in plane strain simulations).
The applied nominal strain rate _ ε is constant, which results in a linearly increasing externally-applied power density P as the corresponding externally applied macroscopic stress σ increases linearly with time until the onset of cracking (see Fig. 4). During this initial stage, the energy is stored elastically in the bulk (term P b in the figure). Due to the simulation technique used (embedded cohesive elements), part of the elastic energy is also stored in the cohesive elements with the cohesive stiffness acting as an elastic spring. This numerical elastic strain energy associated to the compliance of the cohesive elements, denoted as P f comp in Fig. 4, is non-negligible because of the large number of cohesive elements. Consequently, the numerical approximation of the elastic power corresponds to the sum of P b and P f comp . It is worth pointing out here that, although the elastic behavior of the computational domain is strongly affected by the large number of embedded cohesive elements (i.e., no elastic convergence for the bulk material), the fracture behavior can converge upon mesh refinement as the fracture process localizes in a convergent crack pattern.   Table 3 Components of the applied macroscopic strain tensors ε expressed in terms of a loading strain parameter λ > 0. Components are referred to the tensor basis derived from the vector basis fe 1 ; e 2 g shown in Fig. 2.

Laterally-constrained uniaxial extension
Mixed deformation: Equibiaxial stretch and pure shear As the cohesive elements reach their cohesive strength, microcracks nucleate in the matrix and in the matrix-fiber interfaces. From that instant, the external power reaches a plateau indicating that the average stress remains approximately constant in time. However, during that time interval there is a redistribution of the loads inside the specimen. The cracked surfaces start to dissipate more energy (at an increased rate P f ), which is indicated by the solid line labeled as P f c , with the sub-index c used to emphasize that it only includes the cohesive elements that are actually cracked (hence excluding the compliant cohesive elements that only deform elastically). During this stage, the net elastic power (bulk plus cohesive compliance) decreases until it reaches zero, indicating that the sample has reached its maximum stored elastic strain energy. Subsequently, the stored elastic strain energy is transferred from the bulk (and from compliant cohesive elements that are not cracked) towards the adjacent cracked surfaces. During the main cracking stage (which occurs while the specimen is still being pulled at a constant rate), the transfer of energy occurs from both the externally-applied power and the adjacent (elastic) material towards the crack. This transfer of energy corresponds to negative values of P b and P f comp , which results in a significant increase in the dissipation rate P f c in the crack. Eventually all power and dissipation rates decrease to zero as the main crack is formed. Further deformation of the sample only results in a translation of the cracked parts which become disconnected. The graphs in Fig. 4 also indicate the points where the viscous regularization becomes noticeable (corresponding to 5% and 10% of the dissipated energy). As can be observed in the figure, the contribution of viscous dissipation is particularly active in the last stage of the cracking process. This indicates that caution has to be exerted interpreting the tail of the response curve as it may contain purely numerical dissipation, albeit a relatively small percentage.
The verification of the scale transition relations is shown in Fig. 4 in terms of the effective rate of dissipation computed using the kinematicsbased method (term P f u indicated with plus symbols) and the kinematicsbased method (term P f t indicated with solid circles). As may be observed, both methods yield the same values as the fracture power P f c , i.e., by construction For clarity, the third method considered (i.e., the hybrid method that approximately satisfies the Hill-Mandel condition) is not shown in the graph but the resulting data points follow approximately the response given by P f c , i.e., This method has been previously verified in Turteltaub et al. (2018).
Through a time-integration scheme of the crack opening rate, effective traction-separation relations can be established for each scale  transition method as shown in Fig. 5 for two distinct values of the viscous regularization, namely 10 À 3 (top figure) and 10 À 5 (bottom figure). The curves correspond to the normal component of the effective cohesive traction as a function of the effective crack opening. The viscous parameter is introduced in the microscale traction-separation relations used in each cohesive element to provide numerical regularization; it acts as a Kelvin-like viscosity that penalizes large values of the crack opening rate, effectively limiting the rate of opening.
Two obvious but important observations may be immediately drawn: (i) each scale transition method provides a distinct traction-separation relation even though all three curves are obtained from the same raw data and (ii) viscous regularization influences the shape of the effective traction-separation relation, particularly for the kinematics-based method. The shape of the effective traction-separation relation obtained from the kinetics-based method is the least affected in terms of the viscous regularization.
In terms of specific features, the kinematics-based method globally preserves geometric information about the onset of matrix cracking at the early stage, which is characterized by a small plateau at a relatively low stress, corresponding to the strength of the matrix material. In contrast, the kinetics-based method eliminates this feature in the effective curve, with the contribution of the matrix cracking being ignored in favor of a stiff initial response dominated by the uncracked fibers. The hybrid method, which approximately enforces the Hill-Mandel condition, provides an average response between the kinematics-based and kinetics-based methods. Indeed, observing the initial response, the hybrid predicts the same fracture strength as the kinetics-based method but reached at the same effective crack opening as the one computed from the kinematics-based method. Correspondingly, in the hybrid method, the contribution of the early matrix cracking is reflected in a less stiff initial response. A qualitative interpretation of the differences between the distinct averaging schemes is provided in the sequel after analyzing two more loading cases.

Mixed equibiaxial and pure shear deformation
The second loading case involves the simultaneous application of equibiaxial extension together with a pure shear deformation aligned with the fiber directions as indicated in Table 3 (i.e., with the principal shear strain directions oriented �45 � with respect to the ½90�-fiber direction. This strain is equivalent to a laterally-constrained extension of magnitude 2λ in the direction þ45 � (measured clockwise from the ½90�-fiber direction). A typical crack pattern is shown in Fig. 6, which shows a periodic crack that involves several types of failure, namely matrix cracking, fiber cracking and fiber pull-out. The main crack branch is indicated in bold whereas the thin lines represent crack branches or isolated crack segments. The local orientation is indicated with the local normal vector. A typical feature, also encountered in similarly-oriented and loaded samples, is that the fiber cracking occurred in the fiber direction and not in the (average) principal strain direction. Matrix cracking occurred mostly in a plane perpendicular to the principal strain direction while a significant fiber pull-out (or fiber separation) is observed for the [90]-oriented fibers. After postprocessing it was found that the effective normal vector m f is oriented in the extension direction, as may be verified through visual inspection in Fig. 6, where the fiber cracking and fiber pullout portions compensate each other to obtain the same effective crack orientation as the crack portion through the matrix. Some samples showed a less significant amount of fiber pull-out, but this was accompanied by fiber cracks with local normal oriented at 0 � (i.e., [90]-direction) and matrix cracking between fibers at þ45 � such that the effective normal also remained at þ  As in the previous example, one can verify that the scale transition methodologies consistently preserve the dissipation of energy across length scales. Indeed, as shown in Fig. 7, the total energy P f c dissipated in the crack coincides with the effective energy dissipated according to the kinematics-based method (P f u indicated with plus symbols) and the kinematics-based method (term P f t indicated with solid circles). For this loading case, the matrix cracks undergo a large opening before the fibers fail. The bulk (and the compliant cohesive elements) initially store elastic strain energy, which is then transferred and dissipated in the main crack as shown in Fig. 7. This process is similar to the one observe in the previous loading case (uniaxial extension). However, the differences in the effective traction-separation relations are more significant, as can be observed in Fig. 8.
Although the traction-separation relations effectively dissipate the same amount of energy (although only approximately for the hybrid method), the main features of the effective relation are significantly different. This is directly related to the definitions of the effective properties since all three relations are obtained from the same raw data (microscale data). As in the previous example, the kinematics-based method preserves the kinematics of the onset and propagation of the initial matrix cracking, with a relatively large effective crack opening and, in order to compensate for the amount of initial amount of dissipation, it predicts a relatively low effective cohesive traction. In contrast, the kinetics-based method predicts a stiff response initially, thus emphasizing the large force transmitted through the (yet uncracked) fibers and consequently ignoring the large matrix crack opening. Subsequently, the onset and evolution of the fiber cracking is also predicted in two rather distinct ways. The kinematics-based method has on average a small increment in its effective opening rate during fiber cracking, hence it compensates by predicting a relatively large cohesive traction. The kinematics-based model slowly decreases the cohesive traction by increasing the effective crack opening, eventually predicting a relatively large final crack opening. The hybrid method, which preserves the kinematics-based approach but only approximately enforces the Hill-Mandel condition, provides an intermediate response, with a predicted fracture strength comparable to the one obtained from the kinetics-based method and a final crack opening similar to the one predicted by the kinematics based method. Similarly, the initial response is a combination of the two other methods in terms of cohesive   crack running horizontally is dominated by the separation of the matrix between plies (delamination) whereas the vertical crack passes through the matrix and matrix-fiber interfaces but is mostly characterized by the fact that it breaks the fibers (fiber cracking). In this case a segment that appears to be common to both cracks was assigned by the identification algorithm to the vertical crack based on the average orientation of that segment. stiffness as shown in Fig. 8.

Equibiaxial deformation: intersecting cracks
The third and last loading case is equibiaxial extension as indicated in Table 3. This case has a new feature compared to the previous loading examples, namely that two intersecting periodic cracks appear in the microscopic volume element as shown in Fig. 9.
A python-base script was developed to identify multiple cracks and segregate them according to their average orientation. The result of the post-processing operation is the identification of a crack running horizontally, which is dominated by the separation of the matrix between the [0] and [90] plies (delamination) and a vertical crack that passes through the matrix and matrix-fiber interfaces but is mostly characterized by the fact that it breaks the fibers (fiber cracking). A small segment that appears to be common to both cracks is assigned by the identification algorithm to the vertical crack based on the average orientation of that segment, although for some other simulations it was assigned to the horizontal crack.
Based on this detection and segregation algorithm, the scale transition methods can be applied separately to each crack, resulting in two distinct traction-separation relations as shown in Fig. 10. The tractionseparation relation on the left figure corresponds to the delamination crack, which is characterized by a relatively small effective strength and a relatively small dissipation of energy. This can be traced back to the fracture properties of the matrix material that are lower than of the fibers. In contrast, the traction-separation shown on the right figure, corresponding to fiber cracking as the dominant failure mechanism, has a relatively large effective fracture strength and fracture energy.
The individual features of the traction-separation relations are similar to the ones observed in the previous two loading cases, namely an overprediction of the fracture strength using the kinematics-based model and an overprediction of the ultimate crack opening from the kinetics-based model. It is also worth pointing out that the hybrid method that approximates the Hill-Mandel condition under-predicts the energy dissipated for fiber cracking as may be observed in terms of the areas under the curves in Fig. 10 (right), which are a visual measure of the energy dissipated during cracking since the tangential components (not shown) have only a small contribution. This example illustrates that it is feasible to extract information about intersecting cracks due to the capacity of the identification algorithm to identify separate cracks. Further, it also illustrates a typical feature of anisotropic fracture mechanics for composite materials, namely that the cohesive relations are dependent upon material orientation, hence a given cohesive relation must be specified only on the corresponding plane of fracture (e.g., delamination or fiber cracking).

Interpretation of the distinct scale transition methodologies
As illustrated in the previous examples, the choice of the effective quantities in a traction-separation relation leads to energeticallyequivalent relations (same fracture energy) but otherwise have distinct features (such as the maximum cohesive traction or the critical crack opening, which may be used as measures for the effective fracture strength). Consequently, it is relevant to provide a simple interpretation for the differences between the kinematics-based method and the kinetics-based method. The hybrid method (approximate Hill-Mandel) is not treated explicitly in this section since it is a combination of the kinematics the kinetics-based methods.
As mentioned before, the surface-based Hill-Mandel condition for hierarchical multiscale analysis of localized mechanical response cannot be satisfied a priori using the methodologies commonly-used for a volume-based approach (e.g., uniform stress or strain on the MVE boundary or periodic boundary conditions). However, both conjugated pairs f½½ _ u�� f Γ ; t f;HM Γ g (kinematics-based) and f½½ _ u�� f;HM Γ ; t f Γ g (kinetics-based) satisfy the Hill-Mandel condition by construction. The kinematics-based method preserves the kinematics of the microstructural element (i.e., the effective crack opening rate is directly a surface average of the corresponding microscopic quantity), but in order to enforce consistency in terms of power, the corresponding effective cohesive traction needs to adapt and hence deviate from the surface average. This framework may be seen as similar to the Voigt-Taylor constant deformation assumption used in micromechanics (in time rate form), but keeping in mind two important distinctions: in the present framework the balance of linear momentum is satisfied (as opposed to the Voigt-Taylor assumption in which it is not) and in the present method the Hill-Mandel scale transition is enforced a posteriori (as opposed to the Voigt-Taylor method in which it is satisfied a priori). Similarly, the kinetics-based method preserves the cohesive tractions (the effective traction is directly a surface average) and enforces the Hill-Mandel condition at the expense of the effective crack opening. In line with the previous comment, it can be mentioned that this framework is somewhat analogous to the Reuss-Sachs constant stress assumption but with the distinction that the kinetics-based method satisfies the kinematic conditions for cracking (which the Reuss-Sachs method does not) and the Hill-Mandel scale transition is satisfied only a posteriori (which is satisfied a priori in the Reuss-Sachs method).
A schematic illustration of the differences between the two averaging methods is shown in Fig. 11 for a composite material in which the matrix material fails first while the fibers fail later, a situation analogous to the examples presented above.
As indicated in the figure, the matrix cracking stage is characterized by an increase in the crack opening in the matrix while the fibers increase their share of the load bearing distribution. Consequently, the kinematics-based method determines a higher effective crack opening and computes a lower cohesive traction compared to the kinetics-based method. Subsequently, as the load increases and reaches its peak, the second stage in the fracture process involves fiber breaking. During this second stage, the situation is reversed, i.e., the kinematics-based method determines a lower effective crack opening and computes a higher cohesive traction compared to the kinetics-based method. This is due to the fact that both methods need to compensate with one adjustable quantity (either traction of opening rate) the significant increase in dissipation rate as the fibers break. Indeed, the kinematics-based method, which fully respects the geometric opening, adapts the cohesive traction by increasing it significantly during fiber breaking, above the value predicted by equilibrium. Conversely, for the kinetics-based method, in which the effective traction decreases during fiber breaking, the adjusted crack opening rate becomes higher than the geometrically accurate one to compensate for the increased dissipation. In closing this section, it is worth pointing out that the hybrid method (approximate Hill-Mandel) typically preserves most of the relevant microscale information from the crack opening and cohesive tractions at the expense of the dissipated energy.
Although each method has some desirable properties, it is shown in the next section that the kinematics-based method should be avoided in general for multiscale analysis.

Multiscale convergence analysis and comparison between scale transition methods
In a hierarchical multiscale analysis, an important step is to verify the convergence towards a representative volume element (RVE) as the microscale volume elements are increased in size. For fracture, this procedure actually requires identifying a representative surface element (RSE) where the phenomenon is localized, which is embedded in a sufficiently large microscopic volume element. As shown in Fig. 12, two types of volume elements are considered in this section, namely ½0� and ½0=90� n layouts. For the convergence analysis, square domains of increasing size are considered (25, 50, 75 and 100 μm). The material properties used are the same as in the previous section as given in Tables 1 and 2. For each size, several realizations were tested with random distributions of the fibers in the ½0� direction while the fibers in the ½90� were kept fixed. Mesh refinement was conducted for each size, with energy convergence being achieved typically with a mesh of about 1 μm or 2 μm without the need for large viscous regularization. However, due to the modeling technique adopted in the present analysis (i.e., embedded cohesive elements) further mesh refinement generated solutions that had a large amount of artificial viscous dissipation in order to obtain numerically converged solutions. This was monitored in terms of the percentage of viscous dissipation compared to the total dissipation.
Consequently, there is a lower limit in terms of mesh size below which divergence of actual fracture dissipation occurs. In that case, the viscous dissipation would need to be excluded from the total dissipation in order to preserve the physically-based fracture quantities. This, however, falls beyond the scope of the present work, hence only results that did not contain a large amount of viscous dissipation were used in the postprocessing.

Fracture on ½0� samples
The first type of samples studied are transversely isotropic arrangements subjected to laterally-constrained uniaxial loading. For each characteristic size of the microscopic volume element, a mesh refinement analysis (not shown here for conciseness) indicated that a mesh with characteristic element size of 1 μm provided a sufficiently converged solution, both in terms of the energy as well as the corresponding traction-separation relation obtained from both methods. In addition, for each characteristic size of the microscopic volume element, three distinct realizations were analyzed. This simulation setup, although limited due to the overall computational cost, allows to study variations both across realizations and across volume element sizes. The results are shown in Fig. 13 in terms of the energy dissipated due to cracking (top figure), the effective traction-separation relation obtained from the kinematics-based method (middle figure) and the traction-separation relation obtained from the kinetics-based method (bottom figure). As may be observed, there is a relatively fast multiscale convergence in terms of the energy dissipated, with the 50� 50 μm samples already providing sufficiently converged results. The accumulated dissipation at the end of the fracture process, which corresponds to the composites' fracture energy for transverse ply cracks, is about 210 J/ m 2 , which is slightly above the fracture energy of the matrix and the matrix-fiber interface (see Table 2 recalling that 1 J/m 2 ¼ 10 À 3 MPa mm). The fact that the effective value (per unit macroscopic length and depth) is slightly above the microscale value (per unit microscopic length and depth) also reflects the fact that the path of the crack is typically not a straight line at the microscopic level.
Similarly, the effective traction-separation relations obtained from both methods converge relatively fast, with the average response of three 50 � 50 μm samples being within the (discrete) standard deviation of the three 75 � 75 μm samples. The predicted effective fracture strength, as represented by the maximum value reached in the tractionseparation relation, shows some fluctuations for the kinematics-based model (middle figure), but does converge relatively fast across sample sizes for the kinetics-based method. As shown in the bottom figure, the converged fracture strength is close to 85 MPa, which indicates that the crack is dominated by the separation between fibers and matrix (see Table 2). Both methods predict similar values for the ultimate effective crack-opening of about 5 μm.
As may be observed in Fig. 13 (middle), the traction-separation relation from the kinematics-based method is prone to numerical inaccuracies with sudden jumps in the cohesive traction that can be traced back to inaccurate instantaneous values of the crack-opening rates. While it is of course possible to apply a smoothing technique, the results are shown in their unfiltered form to highlight the limitations of applying the kinematics-based approach directly. In contrast, the traction-separation response predicted from the kinetics-based approach is relatively smooth, which reflects that instantaneous inaccuracies in the calculated value of ½½ _ u�� have a less noticeable effect on ½½u�� since the total crack opening is obtained by integration in time. Although both methods appear to be viable in terms of predicting effective properties, it is shown in the next example that the kinematics-based may not be a reliable approach, at least not in its current form.

Fracture on ½0=90� n samples
The second type of samples analyzed are cross-sections of ½0=90� n layouts subjected to laterally-constrained axial extension in the direction of the [90] fibers. The results of the simulations for samples of various sizes are shown in Fig. 14 in terms of energy dissipated (top figure), effective traction-separation relations using the kinematicsbased method (middle figure) and effective traction-separation relations using the kinetics-based method (bottom figure).
In this case, converge of the microscopic volume element within a given tolerance is achieved for the samples of 75 � 75 μm in the sense of the effective fracture energy. Indeed, the predicted energy dissipated for the 75 � 75 μm and the 100 � 100 μm were both close to about 2400 J/ m 2 at the end of the fracture process (i.e., composite fracture energy). However, since energy density is a global measure for the volume element, its convergence across MVEs of distinct sizes does not imply convergence of a traction-separation relation, which is a more stringent requirement as it contains more detailed information. From Fig. 14, it can be seen that post-processing the data in distinct ways has an important effect: whereas the kinematics-based method shows no convergence, the same data shows a convergent traction-separation relation using the kinetics-based method. In fact, the kinetics-based method had no convergent traction-separation relation even at the level of mesh refinement (not shown here for conciseness). However, using the same data, the kinetics-based method revealed a convergent traction-separation relation upon mesh refinement up to a lower limit when viscous regularization prevented global convergence.  As indicated above, the data deteriorated upon further mesh refinement due to viscous effects, providing non-convergent results also when using energy as a measure. That problem, however, is related to the modeling technique (i.e., embedded cohesive elements) and not to the scale transition methods. The conclusions drawn from the simulations pertain to cases without a significant viscous regularization. In this case, from the kinetics-based method, the simulations predict a composite fracture strength of about 1330 MPa in the [90]-direction, which shows the significant influence of the fibers properties (see also Table 2). In contrast, the kinematics-based method does not yield any conclusive value for the fracture strength.

Concluding remarks
In the present work, two distinct length scale transition methods are developed and implemented for the post-processing of microscale data into macroscale quantities, in particular in the form of an effective traction-separation relation for a composite material. The two methodologies, termed the kinematics-based and kinetics-based approach, satisfy the necessary conditions to account for energy dissipation from the micro to the macroscale. Results from the two methods are compared to a previously-proposed scale transition method, termed the hybrid method. Furthermore, the two methodologies are used to analyze the detailed microscale data across volume elements of increasing size in a study to establish numerically the existence of a representative surface element. The following conclusions may be drawn based on the numerical experiments with samples of fiber-reinforced composites: � The preferred approaches to establish effective quantities rely on averaging the tractions across the cracked surface (kinetics-based method) or a combination of this with a volume-averaged stress in the neighborhood of the crack surface (hybrid method). � The kinematics-based method that establish the effective crack opening as an average of the crack opening along the cracked surface should in general be avoided as it fails to provide conclusive (convergent) results. However, the average of the crack opening along the cracked surface may be used as effective crack opening in the hybrid method (i.e., by relaxing the condition of the scale transition requirement and only enforcing it approximately). � The choice of the effective quantities in a multiscale fracture mechanics problem is critical since the same microscale data was shown to converge for the kinematics-based method and diverge for the kinematics-based method, despite that both methods satisfy the same Hill-Mandel scale transition condition. � For intersecting cracks, it is possible to extract separate effective traction-separation relations that may be associated with a different fracture mechanism in a composite material that is anisotropic from the fracture point of view.
In addition to the aforementioned conclusions, it is worth indicating that the use of viscous regularization needs to be carefully monitored to distinguish between physically-relevant quantities and purely numerical artifacts in the simulations. In principle, regularization of the original formulation must be kept to the minimum to avoid non-physical results, but large enough to avoid singularities as non-convergent results cannot in general be relied upon. Another possible source of inaccuracies is related to the boundary conditions used (in the present case periodic boundary conditions). This topic falls outside of the scope of the present work, but it is worth mentioning that, since in general the predictions are influenced by the type of boundary conditions chosen, an exhaustive analysis is still required across different types of conditions. It is also relevant to mention that the present analysis does not include failure by compression, which typically involves mechanisms such as fiber buckling, fiber kinking and/or shear sliding. In those cases, an additional contact algorithm is required to avoid overlapping of (failed) cohesive elements, which is however outside of the scope of the current simulations. Despite the aforementioned list of open challenges, the present work demonstrates that choosing the right scale transition method is critical for the proper interpretation of the microscale data.
In terms of implementation in a multiscale computational environment of the scale transition averaging methods developed here, there are two general approaches that may be used: (1) an interactive hierarchical multiscale framework with simultaneous macro and micro-level calculations (often referred to as an "FE-squared" framework) and (2) a single-scale macroscopic framework together with a pre-computed response (sometimes referred to as an "offline" approach), in which a macroscopic model has been calibrated through suitably-chosen microscale simulations. In both approaches the averaging methods developed here may be used. In the first case it is typically used in an incremental or rate-form, possibly embedded in an iterative macroscopic algorithm for implicit time-integration schemes. In the second case, it may be applied as indicated in the present work, with the goal of calibrating a macroscopic model (see, e.g., the approach van Hoorn (2016) for the case of a macroscopic cohesive relation which also includes characteristics of the microstructure).
It is also worth mentioning that, depending on the multiscale framework used, care has to exerted to guarantee consistency between the macroscopic quantities used as input and the effective quantities obtained through postprocessing of the microscopic simulations. In particular, an iterative method would be required in an FE-squared framework to guarantee that an increment in the macroscopic crack opening displacement used as input for a microscale calculation actually matches the effective crack opening obtained from the kinetics-based method. This issue, however, falls outside the scope of the current work.