A Proper Generalized Decomposition (PGD) approach to crack propagation in brittle materials: with application to random field material properties

Understanding the failure of brittle heterogeneous materials is essential in many applications. Heterogeneities in material properties are frequently modeled through random fields, which typically induces the need to solve finite element problems for a large number of realizations. In this context, we make use of reduced order modeling to solve these problems at an affordable computational cost. This paper proposes a reduced order modeling framework to predict crack propagation in brittle materials with random heterogeneities. The framework is based on a combination of the Proper Generalized Decomposition (PGD) method with Griffith’s global energy criterion. The PGD framework provides an explicit parametric solution for the physical response of the system. We illustrate that a non-intrusive sampling-based technique can be applied as a post-processing operation on the explicit solution provided by PGD. We first validate the framework using a global energy approach on a deterministic two-dimensional linear elastic fracture mechanics benchmark. Subsequently, we apply the reduced order modeling approach to a stochastic fracture propagation problem.

parameters, such as loading conditions, specimen geometry, material properties, etc. Taking into account such uncertainties in an analysis typically implies that the number of times that a solution must be computed increases rapidly with an increase in the number of uncertain parameters. The use of reduced order models is then indispensable as these make it practical to solve the problem for many parameter realizations at an affordable computational effort.
While Reduced Order Modeling (ROM) is a well-established concept in the field of linear elastic solid mechanics [4,6,19], its application to fracture mechanics problems has remained essentially unexplored, with Ref. [25] providing a notable exception. In the present work, a new ROM technique for fracture propagation is presented which allows failure to be studied as a post-processing operation of a parameterized solution that incorporates varying loads, crack lengths and material uncertainties. We propose a parameterization of the crack on the one hand, and a method to take into account the fracture propagation criterion in the reduced order model setting on the other hand. Furthermore, we extend the framework to include random heterogeneities in the material properties.
The reduction method of choice in this work is the Proper Generalized Decomposition (PGD) method, which is a reduced order modeling technique specifically designed to counter the curse of dimensionality induced by the increase in system parameters to be considered in an analysis [10]. The key idea of the PGD technique is to represent the generalized solution in the whole computational vademecum [28,31] (i.e., the high-dimensional parameter space) as a finite sum of terms that involve the product of functions of the system parameters. The computation of this generalized solution is referred to as the offline stage. Once the generalized solution has been obtained, the solution space can be browsed in a computationally efficient way, making it suitable for real time computations [8,22]. This evaluation of the solution space for a particular set of system parameters is referred to as the online stage.
Our work is based on the concept of linear elastic fracture mechanics (LEFM), which is a frequently used model for brittle fracture [20]. We consider Griffith's fracture propagation criterion, which evaluates the stability of a fracture based on an energy balance between the work done by external loads, the elastic energy stored within the system, and the energy dissipated through the fracture surface. Griffith's theory in its basic form is restricted to elastic brittle materials in which there is no plastic deformation near the crack tip. The simulation of fracture evolution in the LEFM framework typically involves a stepwise incrementation of the crack path based on the evaluation of the fracture criterion, which implies that a linear elasticity problem (with a tip singularity) must be solved at each step in the propagation process. This finite element procedure is typically computationally expensive because, on account of accuracy and stability requirements, the crack length increments must generally be small, and because some form of mesh adaptation is required to accommodate changes in fracture geometry. The PGD approach in this work conveniently bypasses these problems, as the fracture length is considered as one of the coordinates of the obtained parametric solution, and differentiation with respect to the fracture length provides a suitable propagation measure in the form of the energy release rate at all configurations in the parametric domain.
This paper is organized as follows. The model problem considered in this work is introduced in Sect. 2. Section 3 demonstrates how a separable form of the problem can be obtained in regard to the fracture length, which is a prerequisite for the application of the PGD method discussed in Sect. 4. We herein adapt the PGD formulation to solve a linear system of equations, which we refer to as the PGD solver [27]. Sect. 5 studies the accuracy of the fracture length parametrization in the setting of a stationary fracture. Section 6 then describes the application of the PGD framework to Griffith's fracture model, along with the consideration of an LEFM benchmark test case [26]. Section 7 then presents an application in the stochastic setting, where we use the Karhunen-Loève expansion [15,23] to discretize random field material properties. A Monte Carlo based stochastic analysis is then performed that demonstrates the efficiency of the PGD framework. Conclusions are presented in Sect. 8.

Model fracture problem
As a model problem we consider a straight fracture in a homogeneous linear elastic two-dimensional (d = 2) continuum, see Fig. 1. The crack propagates in response to an external traction imposed on the system. Inertia, gravity and body forces are neglected. Assuming small deformations and deformation gradients, along with plane strain assumptions, the solid deformation is governed by the momentum balance where the Cauchy stress, σ , follows Hooke's law for isotropic materials σ = 2μ ε + λ tr(ε) I, where u = (u x , u y ) denotes the displacement field, and ε the infinitesimal strain field. The Lamé parameters μ and λ are directly related to the Young's modulus, E, and Poisson's ratio, ν. Exploiting the symmetry of the two-dimensional model, the boundary conditions are given by where n is the outward pointing normal vector and t is the imposed boundary traction. Defining the function space for the vector-valued displacement field as the weak form of the problem reads as follows: The bilinear and linear operators in (2)  where C is the fourth-order elasticity tensor in accordance with Hooke's law (1), i.e., σ = C : ε.
The finite element discretization of the displacement field is given by where {N i (x)} n i=1 denotes the set of n vector-valued finite element basis functions that conform to the space V, and {û i } n i=1 are the corresponding coefficients. Discretization of the weak problem (2) then yields the linear system of equations where the vectorû = (û 1 , . . . ,û n ) contains the solution coefficients, and the coefficients of the stiffness matrix K and load vector f are given by: Evidently, the finite element problem (5) depends on the parameters of the model. In the case that one is interested in a single parameter configuration, this would simply require the assembly of the finite element system for that particular setting, and then to solve that system to find the approximate solution. In the context of this work, however, the central idea is that the system (5) must be assembled and solved for many different parameters. To this end, we introduce the parametric solution to the problem, u(x; μ), where the (scalar) problem parameters μ = (μ 1 , . . . , μ n μ ) are defined over the parameter domains I μ = I μ 1 × · · · × I μ nμ .
The pivotal idea of the PGD method is to attain u(x; μ) as the solution to a problem posed on the higher-dimensional domain × I μ , the spatial semi-discretization of which can be written as: The general PGD strategy to obtaining this solution is to formulate a higher-dimensional weak form problem corresponding to (2), and then to discretize this higherdimensional problem in space and in the parametric dimensions; see, e.g., [9,10] for the fundamentals of PGD. An essential aspect of the PGD framework is that in order to efficiently compute the parametric solution, a separable form of the weak form problem (or its discrete version) must be available. With respect to the spatially discretized system (5) this means that the stiffness matrix and force vector should be of the form, where n k and n f denote the total number of terms needed to represent the parametric stiffness matrix and parametric force vector, respectively. Note that when these affine representations are not available, it is possible to construct affine A non-standard aspect in relation to the fracture problem considered in this work, is that the crack length parameter, l c , enters the problem through the definition of the domain. As a consequence, the separable forms (8), with l c as one of the parameters, will not follow naturally from (5). Obtaining separable forms instead requires recasting of the formulation in a canonical form through a pull back of the problem to a reference configuration. This reformulation of the problem is considered in the next section.

Fracture length parametrization
In this section we consider the parametrization of the system of equations with respect to the fracture length, l c ∈ I l c = [l min c , l max c ]. For the sake of simplicity, we here consider this fracture length to be the only parameter, such that (8) reduces to: The matrices K i and the vectors f i do not depend on the parameter l c , and the functions φ i (l c ) and ψ i (l c ) depend on the parameter only. In order to determine the parametric forms in (9), a reference domain and a mapping function are introduced as illustrated in Fig. 2. The mapping function, M : ref → , which depends on the parameter l c , transforms the parameterindependent reference domain, ref y), where the length of the crack is equal to l c . Through this mapping, the crack length can be described by applying the corresponding mapping to the reference domain. We here consider the following choice for the mapping x = M(X, l c ): The Jacobian of this mapping follows as: The inverse of this Jacobian can be obtained analytically and allows for an exact separable representation as the sum of products of matrices that do not depend on the parameter l c and functions that depend only on that parameter: A separable form of the determinant of the Jacobian can similarly be obtained: The matrix and vector components in Eq. (6) can now be transformed via the mapping M(X, l c ) into equivalent integrals over the reference domain as where use has been made of the operators defined in (  A fundamental difference between the finite element discretization over the reference grid, Eq. (14), and the system obtained using a direct discretization over the physical domain, equation (6), is that the crack length parameter in (14) appears inside the integrands of the matrix components, and not in the domain boundary (and constraints) definitions. This makes it possible to obtain the separable forms of the stiffness matrix and force vector required for the PGD framework.
Substitution of the definitions of the inverse Jacobian (12), and the determinant of the Jacobian (13) into Eq. (14) yields a system of the form (9). From this substitution it directly follows that the separable form of the stiffness matrix is composed of n k = 4 parametric basis functions: The corresponding stiffness matrices are obtained as: Similarly, n f = 2 parametric shape functions are found for the force vector: The corresponding vector components are found as: The system composed of these separable forms for the stiffness matrix and force vector assumes the canonical form (7).

The Proper Generalized Decomposition (PGD) method
The parametric problem (7) is solved here using the Proper Generalized Decomposition (PGD) method [2,3,8]. The particular use of the PGD method considered here follows the idea presented in [13,27], where the method is applied to a discretized (in both space and parametric dimensions) system of linear equations. This differs from the standard use of PGD, where the method is applied to the weak form of the problem (e.g., [12,24,28,31]). The separated form of the PGD approximation,û pgd (μ), takes a form similar to the separated versions of the stiffness matrix, K, and external force vector, f, in Eq. (8), viz.: where the vectorsû i , for i = 1, . . . , n pgd , are constant vectors of the same size as a standard spatial finite element solution, and the scalar functions G i j (μ j ) are independent of space with μ 1 , μ 2 , . . . , μ n μ as parameters and n μ being the total number of parameters. Note that the parametric functions G i j (μ j ) are represented discretely by a nodal vector associated with a mesh over the parameter domains I μ j in accordance with where is the set of linear finite element basis functions over the parameter domain I μ j , and whereĝ i j = (Ĝ i j,1 , . . . ,Ĝ i j,m j ) is the corresponding vector of coefficients. In Eq. (18) the vectorsū i and functionsḠ i j (μ j ) are the spatial and parametric modes normalized with respect to the Euclidean norms û i and ĝ i j , respectively, such that the modal amplitudes, β i , are given by: We employ the PGD solver algorithm as presented in Ref. [27], the main ingredients of which are: -The PGD algorithm requires the determination of separable forms of the stiffness matrix and force vector as input. As discussed in detail in Sect. 3, the discrete operator K(l c ) for the parametric problem with the crack length l c as a parameter admits an exact separable representation. This is not generally the case, as we will discuss, for example, in the stochastic test case considered in Sect. 7. In situations where the linear system cannot be separated analytically, it is often replaced by a separable approximation (e.g., [30,31]). There exist several methods to compute such separated approximations. For higher-dimensional parameter domains various methods have been proposed in the literature, such as: an approximation based on the PGD concept [14], Singular Value Decomposition (SVD) type approximations [11], approximations based on the CANDECOMP/PARAFAC methods [7,18], and Tucker decomposition type approximations [29]. An overview of these techniques can be found in, e.g., Ref. [21]. It is noted that in the case of high-dimensional parameter domains, the computation of separable forms can be computationally demanding. -A greedy algorithm [1,8] is used to sequentially compute the terms to the PGD approximationû pgd in Eq. (18). Given the PGD approximation with n pgd − 1 terms, here denoted bŷ an enrichment termû n pgd n μ j=1 G n pgd j is computed as to obtain the PGD approximation with n pgd terms: Each enrichment term is computed one at a time, constructing the summation progressively until the convergence criterion is met with a user-defined tolerance of glob . Each step in the greedy algorithm, i.e., computing each of the enrichment terms, involves the computation of the enrichment modes in space,û i in discrete form, and in the parameter spaces, G i j (μ j ). We herein compute these enrichments iteratively using an alternate direction solver, which is discussed in detail below.
-An alternating direction solution strategy [9] is used to compute the enrichment termsû n pgd n μ j=1 G n pgd j . Leveraging the separable forms, in this alternating direction strategy the spatial and parametric directions are treated sequentially as to reduce the higher-dimensional parametric problem to a series of low dimensional problems. This iterative process is repeated until a fixed point is reached within a defined tolerance. For the explanation of this alternating direction strategy we will consider n μ = 1 with the fracture length μ 1 = l c as the only parameter. For the alternate direction solution strategy, the parametric problem (7) is considered in its weighted residual form: The unknowns in this system are the spatial and parametric enrichment modes,û n pgd and G n pgd l c (l c ), respectively. The corresponding test functions are defined as: In the alternate direction strategy, the system (24) is solved per spatial or parametric dimension: -Given an approximation (or initial guess) for the parametric enrichment mode G n pgd l c , the system (24) reduces to the linear system: Using the separable forms for the stiffness matrix and force vector in equation (9), this system can be rewritten as with n k = 4 and n f = 2. An essential idea of the PGD method is that the parametric integrals in this equation can be evaluated efficiently on account of the fact that these are low-dimensional integrals (in this particular case one-dimensional). We herein use a standard trapezoidal integration rule for the evaluation of these integrals. -Given the spatial enrichment modeû n pgd computed through the system (27), the parametric enrichment mode G n pgd l c can be obtained from the system (24).
From (24) it follows that for all δG n pgd l c (l c ): Equivalently, it holds that for each fracture length l c from which the parametric enrichment mode follows directly as: Substitution of the separable forms for the stiffness matrix and force vector then finally yields: This expression for the parametric enrichment mode can be evaluated quickly by virtue of the fact that the dimensions are separated in the sense that it is not required to reassemble the finite element system for each fracture length. The parametric enrichment mode is represented discretely by projection onto the parametric basis in Eq. (19). Since this discretization pertains to a linear finite element basis, the coefficientsĝ n pgd l c can be computed by evaluation of Eq. (31) in the parametric nodes.
The above alternate direction steps are repeated until the relative difference between two successive steps is smaller than a prescribed tolerance, local , with the subscript iter denoting the alternate direction step, and with the norms defined as:

Numerical analysis of the PGD approximation behavior
Before considering the application of the PGD framework to fracture problems, in this section we first present a numerical study on the approximation properties of the PGD expansion introduced above. We specifically study the convergence behavior of the approximation under finite element mesh refinement, and the approximation behavior with respect to the number of PGD terms, n pgd . All results presented in this section are based on the consideration of the fracture length, l c , as the single quantity to be parametrized. Table 1 lists all parameters that are fixed throughout this section.
In the setting considered here, the separable form derived in Sect. 3 is exact up to integration accuracy. Since the integrals are herein evaluated with Gauss schemes of sufficiently high degree, the separable forms are accurate up to floating point precision. In general, however, the separable form (9) is not exact, as we will consider, for example, in the context of the stochastic analysis presented in Sect. 7. An important first step in studying the approximation behavior of the PGD approximation is then to study the accuracy of the separable form (9). This accuracy can be assessed by comparison of the matrix and right hand side obtained through the separable form (9) with their corresponding original finite element counterparts. Evidently, one has to perform this accuracy assessment in such a way that the parameter variations admitted by the PGD expansion are properly taken into account.

Spatial mesh size dependence
We first study the dependence of the PGD approximation (18) on the spatial finite element mesh size parameter, h, defined as the average element size in horizontal direction (h = H x /n elems,x ). For the discretization of the parameter domain, I l c , we consider 136 elements, and we use the PGD solver presented above to obtain an expansion comprising n pgd = 10 terms. In Fig. 13 the various components of this expansion are illustrated, viz. (a) the spatial modesû i , (b) the parameter modes G i l c (l c ), and (c) the amplitudes β i . The amplitudes convey that the influence of the modes decreases significantly for increasing mode numbers, indicating that the displacement of the system is well characterized in the considered setting with 10 modes. A detailed study of the dependence of the PGD approximation on the modes is considered below (Fig. 4).
To study the approximation behavior of the PGD expansion, we consider the relative energy error with respect to the original finite element solution: whereû pgd (l c ) is the parametric solution provided by PGD andû(l c ) is the solution provided by the direct FE analysis (5) when the parameter is fixed to the value l c . Note that while the evaluation ofû pgd (l c ) for a certain crack length l c involves merely the evaluation of the PGD expansion (18), the computation ofû(l c ) involves the assembly and solution of a finite element system. In addition to the parameter-dependent error (34) we consider the mean energy error over the parameter domain: In contrast to (34), this error measure provides one scalar error value for the complete parametric solution and has no dependency on l c . Figure 5 displays both error measures for various spatial mesh sizes, h, and a fixed parametric mesh size h l c ≈ 0.015 m. The parameter dependent error (34) displayed in Fig. 5a conveys that for a certain mesh size, the error in the PGD solution is dependent on the crack length. The reason for this is that the uniformity of the mesh in the physical domain is affected by the parameter-dependent mapping function (10), which in general causes the error to increase when the crack tip position deviates from l c /H x = 0.5 (i.c., l c = 2) provided that the mesh resolution is of sufficient accuracy. The error e pgd (l c ) is especially significant at the boundaries of the parameter domain, I lc , because at those points the non-uniformity caused by the mapping onto the physical domain (see Fig. 3) is largest. When we compute the mean of the error e pgd (l c ) over the complete parameter domain, i.e., error measure (35), we observe from Fig. 5b that this mean energy error is essentially independent of the mesh size for the finer meshes (h 0.25). This conveys that for these meshes the studied error is dominated by the PGD approximation, which is expected, as we compare the PGD solution with the FE solution on the same mesh.
To study the mesh size contribution to the PGD approximation error, in Fig. 6 we display the mean L 2 error between a PGD approximation u pgd (x; l c ) computed with mesh size h and a PGD approximation, u pgd (x; l c ), with a high resolution mesh with h = 0.03125: Both the number of PGD terms and the discretization of the parametric mesh are identical for both of the compared solu-tions, so that this error measure pertains to the mesh size contribution only. For comparison the finite element convergence plots for various settings of the fracture length are displayed in Fig. 6. This comparison conveys that the PGD solution converges with the mesh size with the same rate as

Parametric mesh size dependence
All results presented above were based on a fixed parametric mesh size of h l c ≈ 0.015 and variations in the spatial mesh size. We now consider the influence of variations in the parametric mesh size under a fixed spatial mesh size of h = 0.0625 m. Figure 8 shows that both the parameter-dependent energy error (34) and mean energy error (35) are virtually independent of the parametric mesh size even on parametric meshes as coarse as h l c = 0.125 m (8 elements). This conveys that, in the setting considered here, the accuracy is governed by the number of PGD modes rather than by the resolution of the parametric mesh. In this section we apply the PGD framework outlined above to the simulation of fracture propagation using Griffith's energy criterion [16]. In Sect. 6.1 we commence with the formulation of the propagation criterion based on the PGD solution. Since the evolution of the fracture is driven by the external load, we herein use the PGD framework to compute the parametric solution with respect to both the fracture length (as already considered above) and with respect to the external load, where λ denotes a load scale parameter such that t = λt witht being a load vector defined ast = (0, 1) MPa. For simplicity in notation, from hereon we denoteû pgd forû n pgd pgd . The separable forms of the stiffness matrix and force vector are a straightforward extension of those in Sect. 3 as a consequence of the fact that the external force vector scales linearly with the load scale λ. As a result, we only have to consider a single linear parametric shape function for the load scale parameter for the force vector in Eq. (8b), such that: In Sect. 6.2 we will demonstrate the application of the PGD framework to a fracture propagation benchmark problem, where the advantages of the PGD framework become apparent as it allows for the fast evaluation of the fracture propagation criterion throughout the evolution process of the fracture, without the need for solving additional finite element problems. For all the simulations we assume plane strain conditions with Young's modulus E = 2 GPa and the other input values taken from Table 1. For the parametric domain of the load scale we use I λ = [6.25 , 62.5]. Furthermore, we define the resultant force F = top t · n d as a quantity of interest, where we assume the specimen to be of unit thickness.

The fracture propagation criterion
We consider Griffith's model [16] for crack propagation in brittle materials. The conceptual idea of this model is that a fracture will propagate if the energy stored in the material is sufficiently large to overcome the fracture energy associated with the creation of new fracture surface. For linear elastic materials an equivalent interpretation of this energy-based model is provided through the concept of stress intensity factors [5]. In the context of the PGD framework we find the energy perspective most suitable, as it provides the possibility to evaluate the propagation criterion directly based on the parametric solution (37).
For a fracture in a given configuration, i.e., with a certain length l c and a given load scale λ, it can be determined whether or not the fracture will propagate by evaluation of the energy release rate. To derive the PGD form of the energy release rate, we consider the energy of the system: The energy release rate is then defined as : When the parametric problem K(l c )(l c , λ)û pgd ≈ f(l c , λ) is solved using the PGD solver with sufficient accuracy, i.e., with small enough tolerances, the energy release rate is given by, According to Griffiths energy balance, a crack will propagate when the energy release rate surpasses the critical energy release rate or fracture toughness, G c , i.e.: This implies that for any crack configuration in the parametric space, i.e., (l c , λ) ∈ I l c × I λ , it can be readily evaluated whether or not the crack propagates. The PGD expansion (37) is crucial in this regard as: (i) The expansion allows for the analytical evaluation of the shape derivatives ∂ ∂l c in Eq. (40), this in contrast to the traditional FE setting, in which this derivative is typically evaluated using alternative techniques (e.g., J -integrals [5]). (ii) Evaluation of the fracture criterion at an arbitrary parametric coordinate is merely an evaluation of the expansion, and hence, does not require the solution of an FE model.

Numerical example: a center-crack under tensile loading
The numerical example discussed here demonstrates the PGD-based evaluation of the energy release rate G in two ways: (i) the energy release rate, G, is used to compute the stress intensity factor; (ii) PGD is used to mimic the fracture propagation process while loading the specimen.

Stress intensity factors
As a means to assess the PGD approximation of the energy release rate, we study the stress intensity factor for a given fracture length l c , and various ratio's of horizontal and vertical specimen dimensions, H x and H y , respectively. The results presented in this section consider the parameters H x and H y as additional parameters in the PGD expansion. The separable forms based on these parameters can be obtained without special treatment, and are omitted here for the sake of brevity. The stress intensity factor is defined as and hence is directly related to the energy release rate (40). The material parameter E in Eq. (42) is defined as E = E/(1 − ν 2 ) for the plane strain problems considered herein. Figure 9 shows the dimensionless stress intensity factors K 1 /K 0 for various parameter configurations, i.e., different l c /H x and H x /H y (see Ref. [26] for a benchmark result). Note that the plotted factors are non-dimensionalized using K 0 = (λt · n) √ πl c , where λt · n gives the magnitude of the applied tensile traction. Figure 9 compares the PGD results based on the settings mentioned in Table 1 for a mesh size h = 0.0625 m. However, note that this plot of non-dimensional stress intensity factors is independent of the input values, i.e., even for different values of geometry and load, similar curves for K 1 /K 0 are obtained. This figure conveys that for the given PGD settings, the stress intensity factor can be computed accurately using the PGD expansion (37). While each point in Fig. 9 would typically represent a finite element simulation in the traditional FEM setting, in the PGD case these are all mere evaluations of the expansion.

Fracture propagation
Now that we have established that the PGD expansion accurately approximates the stress intensity factor, we will here use it to predict the evaluation of the loading force under fracture propagation. To this end, we define the energy functional such that we can distinguish between three cases in the energy landscape over the I l c × I λ parameter domain: 1. The region where the crack is stable: 2. The region where the energy balance is critical: 3. The unstable propagation region: The energy landscape is plotted in Fig. 10a along with the values indicating the energy in kJ of the system. Note that plotting this landscape is computationally feasible using the PGD expansion, but would require a large number of FE solves in the case of a non-reduced model. The presented results are based on the assumption of plane strain conditions with material parameter E = 2.01 GPa and the other settings listed in Table 1  For a particular load scale, until the critical point is reached the crack is stable (green region in Fig. 10a), and beyond the maximum point the crack is unstable (red region in Fig. 10a). The critical energy states are connected in the form of a curve which gives the critical load value for each fracture length. This curve can be identified in Fig. 10a as the line separating the green area from the red area. The key insight is to recognize that, for a shorter crack length, which is left of the critical value point, the total energy (43) of the system increases with increasing crack length. Therefore, additional energy must be stored into the material before the crack can propagate, and hence the crack is stable. However, at longer crack lengths, which is right of the maximum value, an increase in crack length leads to a decrease in total energy, which therefore leads to unstable crack propagation. Evidently, the load-bearing capacity of the specimen decreases as the fracture propagates. Fig. 10 Representation of the loading and fracture evolution process in terms of a the energy landscape and b the force-displacement curve. The elastic loading branch is labeled as I., whereas the softening/propagation branch is labeled as II. The observed critical loading force of F c ≈ 36.3 MN is in agreement with equation (44) and the corresponding stress intensity factor reported in Fig. 9 A common way of representing the fracture evolution process is by plotting the load versus the average displacement of the loading boundary, which is depicted in Fig. 10b for a initial crack length of l 0 c = 2.495 m. Note that the elastic loading branch (label I. in Fig. 10) corresponds to the region where the crack is stable, i.e, the force varies with ∂E ∂l c < 0. The resultant force at which the crack becomes unstable, i.e., when ∂E ∂l c = 0, is defined as the critical loading force, F c . This corresponds to the maximum force in Fig. 10b. This critical loading force is related to the dimensionless stress intensity factors of Fig. 9 by: The softening branch (label II. in Fig. 10) reflects the critical values in Fig. 10a for l c ≥ l 0 c . This part of the curve resembles the unstable propagation part of the process. The total area under the force displacement curve represents the energy carried by the system, which, upon complete failure is equal to the total energy dissipated by the fracturing, i.e., G c (H x − l 0 c ). Such force-displacement curves can be plotted for all l 0 c ∈ I l c by virtue of the explicit availability of the energy functional in (43) in the PGD framework.

Application to fracture propagation in random heterogeneous materials
In this section we extend the PGD framework for crack propagation to a stochastic setting. We introduce randomness in the material properties by representation of the Young's modulus by a random fieldẼ(x), where the tilde indicates the randomness. A truncated Karhunen-Loève expansion [15] is used for the parameterization of the Gaussian fieldẼ(x), which is defined as where μ E is the stationary mean of the Young's modulus and where ξ α and r α (x) are the eigenvalues and eigenfunctions corresponding to the spatial covariance function σ 2 E ρ E (x 1 , x 2 ), with σ E the stationary standard deviation. The autocorrelation function is taken as where x 1 and x 2 are two points in the domain and l E is the correlation length. The n kl Karhunen-Loève modes, R α (x) = √ ξ α r α (x), in Eq. (45) are scaled by independent standard normal random variablesz α . On account of (45) the Young's modulus at any fixed location,Ẽ(x), is normally distributed. The variation σ 2 E is selected such that physically impossible negative realizations are avoided. Although not considered herein, the PGD framework can be applied without modification to, e.g., log-normal random fields. It is noted that we herein construct the random field over the computational domain, thereby implicitly assuming that the random material properties adhere to the symmetries of the homogeneous problem. Preservation of the symmetries is in line with the considered parametrization of the fracture problem, as non-symmetries would result in deviations of the fracture path from the x-axis. Although such variations are evidently physical, consideration of these within the PGD framework is beyond the scope of this manuscript.
In the context of the stochastic analysis considered here, we use the PGD framework to compute the parametric solution with respect to the fracture length, external load, and with the random variablesz α that parametrize the random Young's modulus field: A prerequisite to apply our framework is to express the stiffness matrix and force vector also in this separated format. The separable forms of the stiffness matrix and force vector required here cannot be obtained in an analytical way like in Sects. 3 and 6. Therefore, in Sect. 7.1 we first discuss how the random heterogeneities, which are parametrized by the random variablesz, can be expressed in a separable form for the stiffness matrix numerically. Furthermore, in Sect. 7.2 we outline the computational procedure for a sampling-based stochastic analysis based on the Monte-Carlo method. This stochastic analysis is highly efficient as it leverages the PGD approximation to quickly compute critical force values for realizations of the heterogeneous field of elastic properties. Numerical results for the stochastic test case are presented in Sect. 7.3.

Separable representation of the random system of equations
The random field (45) enters the formulation through the elasticity tensor in the bilinear operator (14a), which, in the context of the stochastic setting considered here, is expressed asC where the constant tensor D depends on the Poisson ratio and on the assumed plane strain state. Since the elasticity tensor is evaluated over the reference domain, the KL modes {R • M} n kl α=1 are pulled back to the reference configuration using the geometric mapping function (10). Since this mapping function is dependent on the fracture length parameter l c , the random elasticity tensor (48) also becomes dependent on the fracture length.
Substitution of the random tensor (48) into Eq. (14a) yields a random stiffness matrix of the form with the stiffness matrix contributions defined as where the index 0 corresponds to the mean contribution, and the index α = 1, . . . , n kl to the stiffness contributions of the KL modes. The separable form (8a) of the mean stiffness matrix (50a) is identical to that presented in Eqs. (15) and (16) with the elasticity tensor set to C = μ E D, which we denote by The derivation of an analytical separable form for the KL contributions to the stiffness matrix, Eq. (50b), is obstructed by the appearance of the geometric mapping, M, in the Karhunen-Loève modes, R i . A semi-analytical separable form can, however, be obtained through the singular-value decomposition of the discretized KL modes. For the construction of this decomposition, we first interpolate the KL modes on the spatial mesh and crack length parameter domain mesh used for the PGD approximation as: The coefficients of this interpolation, represented by the matrixR α , are computed using the KL modes constructed on a significantly refined mesh compared to that used for the PGD approximation. Since (bi)linear Lagrangian basis functions are used for both the spatial domain and the parameter domain, the coefficients are determined by evaluation in all nodal coordinates, (X, l c ), in the higher-dimensional parameter domain, where the mapping (10) is used to transfer data between the physical domain and the reference domain. The interpolation (52) on the mesh used for the PGD approximation is convenient from an implementation perspective, but the usage of this specific mesh is not necessary to attain the separable form of the stiffness matrix. A separable form of the discrete KL modes (52) is then obtained through the singular-value decomposition where σ (α,β) is the β-th singular value for KL mode α, and whereĥ (α,β) andm (α,β) are the corresponding spatial and parametric modal vectors, respectively. For reasons of efficiency this singular-value decomposition is truncated to a number of terms, n svd , that is significantly smaller than the total system size. Substitution of this decomposition into Eq.
(52) then yields the singular-value decomposition for the KL modal functions, where the modal functions are defined as The singular value decomposition of the Karhunen-Loève modes (54) involves two approximations, viz.: (i) an approximation related to the interpolation step (52); and (ii) an approximation associated with the truncation of the decomposition (53). Now that we have obtained an approximate separable form for the KL modes in the form of Eq. (54), separation of the stiffness matrix follows from substitution of this decomposition into the KL stiffness matrix contributions (50b): The components of the matrices K (α,β) (l c ) are given by: Since the spatial modes, h (α,β) (X), are independent of the parameter l c , the matrices K (α,β) can be separated analogously to the Eqs. (15) and (16) with the elasticity tensor set to C = Dh (i,β) (X). Similarly to the separable form of the mean stiffness contribution in Eq. (51), we express this separable form as: Substitution of this separable form for the SVD mode β into Eq. (56) then yields with n k = 4 in accordance with Eq. (15). Further substitution into the expansion of the random stiffness matrix (49) gives: Note that this equation is of the same form as the separable form (8a), with the parameter functions given by combinations of the functions in (15), the random variables,z α , and the singular-value modes for the length parameter, m (α,β) . From (60) it is observed that the total number of terms in the separable form is equal to n k (1 + n kl n svd ). Since the stiffness contributions K i 0 and K i (α,β) are independent of the considered parameters, these can be precomputed. Hence, construction of the stiffness matrix in the PGD solver requires evaluation of (60) only, and not the assembly of a finite element system.

Monte Carlo analysis of the critical load
Using the separable form for the stiffness matrix as discussed in Sect. 7.1, the PGD solver discussed in Sect. 4 is used to attain the PGD solution (47). We here use this parametrized solution to perform a Monte Carlo simulation to attain the probability distribution and statistical moments of the critical loading force for specimens with various initial fracture lengths.
To construct the PGD solution (47) it is necessary to consider a finite dimensional domain for the random parameters, z, which parametrize the Karhunen-Loève expansion for the Young's modulus (45). We herein truncate the random domain to Iz i = [−5, 5] for i = 1, . . . , n kl , based on the idea that realizations beyond this range are unlikely and will have a minor effect on the mean and standard deviation of the critical force. We generate realizations of the uncorrelated random variablesz using a random number generator, and we discard realizations outside of the truncated random domain.
Using the realizations of the random variablesz we then employ Griffith's fracture model as discussed in Sect. 6 to compute the corresponding critical forces, F c . The mean and standard deviation for the critical force are then obtained as where n sample is the Monte-Carlo sample size. In a typical FE-based Monte Carlo simulation, evaluation of the critical loads is computationally demanding, which practically restricts the sample sizes that can be considered. Therefore, in such cases, a sample size is selected that strikes an adequate balance between the confidence level of the attained statistical moments and the required computational effort. In the PGD setting considered here, the computational effort involved in determining the critical force for a given realization of the random field is negligible compared to the corresponding full finite element simulation. This allows for the consideration of sample sizes that are orders of magnitude larger than those that could be considered using direct FE analysis, which in turn enables the computation of the statistical moments with confidence levels that are practically beyond the reach of direct FE analyses. Evidently, the selection of the sample size should be based on a trade-off between the error in the PGD approximation and the confidence level of the Monte Carlo method.

Numerical example: a center-crack under tensile loading
We consider the same numerical experiment as introduced in Sect. 6.2 (see Table 1), but now with a random field of elastic properties. For the random field (45) we set the mean to μ E = 2 GPa and the standard deviation as σ E = 0.2 GPa (a coefficient of variation of 10%). We consider moderate spatial fluctuations in the random field by selecting the correlation length in Eq. (46) as l E = 1.5 H x = 6 m. The parameter domain for the load scale is taken as I λ = [6.25, 62.5].
We consider a Karhunen-Loève discretization consisting of n kl = 3 modes, which are shown in Fig. 11. In Fig.12 we show two realizations of the KL expansion, as well as a sampling-based reconstruction of the auto-correlation function (46). On account of the low spatial frequency of the variations, the KL expansion with only 3 terms is observed to already appropriately reproduce the auto-correlation function.
Using the tolerances specified in Table 1, the PGD solution (47) is truncated at n pgd = 27 terms. The various components of the PGD solution are displayed in Fig. 13. From the modal amplitudes it can be observed that the PGD approximation based on 27 terms approximates the finite element problem well, in the sense that the amplitudes of even higherorder modes will be negligible compared to the considered modes. Figure 14 displays the probability distribution of the critical load for various settings of the initial crack length. The displayed results are based on a sample size of 5000. Note that for each of the displayed subplots in Fig. 14 a separate Monte Carlo simulation is required, which would be computationally impractical using a direct FE approach. The efficiency with which realizations can be computed from the PGD approximation (47) allows us to perform Monte Carlo analyses for different settings in the parameter space. This results, for example, in the evaluation of the critical force versus the initial crack length as displayed in Fig. 15a. The confidence level of the mean values displayed in this plot is approximately 98% based on a sample size of 5000 realizations. Such confidence levels are impractical to obtain using direct FE Monte Carlo. Figures 14 and 15 show that the average critical load bearing capacity decreases with an increase in crack length, while a decrease in the standard deviation is observed. The deterministic result is plotted for reference, from which it is observed that the computed mean is slightly smaller than the deterministic value. The observed results from the Monte Carlo simulation are in good agreement with perturbation analysis results (see [17] for an overview) based on the ana- Fig. 12 (a, b) Examples of realizations of the random elasticity field in accordance with (45). c Reconstruction of the auto-correlation kernel (46) lytical fracture loads for homogeneous specimens, which is to be expected on account of the considered low spatial frequency of the random input.
The Monte Carlo analysis allows us to inspect which realizations of the input lead to a certain response in terms of the fracture load. Figure 16 shows three interesting realizations for the case of an initial crack length of l 0 c = 1 m and a coefficient of variation of the Young's modulus of 10%, viz.: a. The realization closest to the mean fracture load of 77.5 MN corresponds to a Young's modulus field which is very close to its mean value everywhere in the specimen. b. The realization with the largest fracture load of 88.5 MN corresponds to a Young's modulus field which is very high throughout the specimen (on average approximately 25% higher than its mean value), and is particularly large near the tip of the initial crack. c. The realization with the smallest fracture load of 66.6 MN corresponds to a Young's modulus field which is very low throughout the specimen (on average approximately 25% lower than its mean value), and particularly near the tip.
In the context of the PGD approach employed in this work it is noted that, in order to inspect these realizations, only the parameters corresponding to the realization (random variable realizations) have to be stored. The input and output corresponding to these parameters is generated through post-

Conclusions
In this work we have proposed a reduced-order modeling technique for a prototypical linear elastic fracture mechanics problem. An essential ingredient in the proposed approach is to introduce the parametrization of the crack through a geometric mapping. For the considered model problem it then follows that a separable form of the stiffness matrix and external force vector can be obtained analytically, which makes it possible to apply the Proper Generalized Decomposition method to obtain a solution to the parametric problem. The suitability and performance of the proposed framework is demonstrated using a series of numerical test cases, starting with a convergence study for the parametric decomposition. This study conveys that the introduced geometric mapping function for the fracture parameter behaves in accordance with the well-understood behavior of the PGD framework. The PGD fracture framework is further demonstrated using two propagating fracture test cases.
In the first test case it is demonstrated how Griffith's propagation criterion can be evaluated efficiently using the PGD approximation. The representation of the fracture length in the PGD solution enables the straightforward computation of the energy release rate, which is in contrast with standard Fig. 15 Dependence of the mean critical force (solid blue line) on the initial crack length with a 98% confidence interval (shaded area) for 10% variation and 5% variation in the Young's modulus finite element methods, which generally require dedicated numerical techniques for the evaluation of the corresponding shape derivative.
In the second test case the PGD approximation is used to efficiently perform a fracture analysis in the presence of random material heterogeneities. Using a singular value decomposition for the interpolation of the random field of elastic properties pulled back to the reference configuration, an approximate separable form of the stiffness matrix is obtained. The random variable coefficients of the Karhunun-Loève field for the modulus of elasticity appear as parameters in this separable form. Since the fracture load can be computed as a post-processing operation on the PGD approximation, Monte-Carlo simulations can be performed with sample sizes (and confidence levels) that are beyond the typical reach of direct sampling-based stochastic finite element analyses.
Although the presented study clearly demonstrates that the PGD framework can be applied efficiently for the simulation of fractures in the considered model problem, the question naturally arises to what extend the proposed technique can be generalized to more complicated fracture problems. In this regard there are two aspects that must be considered in particular: -While the considered fracture is parametrized by a single variable, namely the fracture length, this is evidently not possible in the case of more complex fractures. Of course, the range of applicability of the proposed technique can be extended to a reasonably sized class of fracture problems using a relatively low dimensional parameter space for the fracture geometry. Think for example of slanted fractures in plane strain or plane stress settings, which, besides the length, would require the fracture angle as an additional parameter. In general, however, representing more complex fracture geometries will rapidly increase the number of parameters, which is detrimental to the performance of the PGD framework. This is particularly the case when one opts to consider a piecewise representation of fractures, which is natural to finite element methods. -For more complex fracture patterns, constructing a suitable geometric mapping function will be considerably more challenging than in the prototypical benchmark considered in this work. Constructing a mapping analytically is very restrictive, but it is very well imaginable that one can construct discrete mapping operators (mapping nodal reference coordinates to nodal physical coordinates). Such more advanced mappings -the construction of which evidently warrants further investigation -will, however, pose several difficulties. For example, the analytical separation of the system of equations as obtained in this work will not be generally obtainable, which hence requires the consideration of potentially computationally demanding approximations for the separable forms.
Moreover, an open research question remains how to deal with fractures with changing topology (e.g., branching, merging), as topological changes can in general not be captured by the proposed mapping technique.
These complications when extending to more complex fractures are evidently very serious. Although future research developments can ameliorate some of these difficulties, obtaining PGD approximations that are able to accurately Fig. 16 Realizations of the Young's modulus field corresponding to the mean fracture load, maximum fracture load and minimum fracture load.
All results pertain to an initial fracture length of l 0 c = 1 m parametrize the complete high-dimensional solution space for complex fracture patterns will likely remain impractical. It should, however, be noted that reduced-order models typically do not serve the role of a direct replacement of highfidelity finite element models. Instead, reduced-order models typically play the role of a relatively cheap surrogate to evaluate approximations of the corresponding high-fidelity model. In this regard it is imaginable that the high-dimensional parameter space associated with the fracture geometry in the finite element model can be reduced significantly, without compromising the properties of the reduced-order model to serve as a cheap approximation of the full model or to provide an improved prior.