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, 20 material properties, etc. Taking into account such uncertain- 21 ties in an analysis typically implies that the number of times 22 that a solution must be computed increases rapidly with an 23 increase in the number of uncertain parameters. The use of 24 reduced order models is then indispensable as these make it 25 practical to solve the problem for many parameter realiza- 26 tions at an affordable computational effort. 27 While Reduced Order Modeling (ROM) is a well-estab- 28 lished concept in the field of linear elastic solid mechanics 29 [4, 6,19], its application to fracture mechanics problems has 30 remained essentially unexplored, with Ref. [25] providing a 31 notable exception. In the present work, a new ROM tech-32 nique for fracture propagation is presented which allows 33 failure to be studied as a post-processing operation of a 34 parameterized solution that incorporates varying loads, crack 35 lengths and material uncertainties. We propose a parame-36 terization of the crack on the one hand, and a method to 37 take into account the fracture propagation criterion in the 38 reduced order model setting on the other hand. Furthermore, 39 we extend the framework to include random heterogeneities 40 in the material properties.
a linear system of equations, which we refer to as the PGD 89 solver [27]. Sect. 5 studies the accuracy of the fracture length an LEFM benchmark test case [26]. Section 7 then presents 94 an application in the stochastic setting, where we use the 95 Karhunen-Loève expansion [15,23] to discretize random 96 field material properties. A Monte Carlo based stochastic 97 analysis is then performed that demonstrates the efficiency 98 of the PGD framework. Conclusions are presented in Sect. 8. 99 2 Model fracture problem 100 As a model problem we consider a straight fracture in a 101 homogeneous linear elastic two-dimensional (d = 2) con-102 tinuum, see Fig. 1. The crack propagates in response to an 103 external traction imposed on the system. Inertia, gravity and 104 body forces are neglected. Assuming small deformations and 105 deformation gradients, along with plane strain assumptions, 106 the solid deformation is governed by the momentum balance 107 108 ∇·σ = 0 in , 109 where the Cauchy stress, σ , follows Hooke's law for isotropic 110 materials 111 σ = 2µ ε + λ tr(ε) I, where n is the outward pointing normal vector and t is the 123 imposed boundary traction. find u ∈ V such that, a(u, v) = ℓ(v) ∀v ∈ V.

133
The finite element discretization of the displacement field 134 is given by  142 where the vectorû = (û 1 ,...,û n ) contains the solution 143 coefficients, and the coefficients of the stiffness matrix K 144 and load vector f are given by: 146 147 Evidently, the finite element problem (5) depends on the 148 parameters of the model. In the case that one is interested in 149 a single parameter configuration, this would simply require 150 the assembly of the finite element system for that particular 151 setting, and then to solve that system to find the approximate 152 solution. In the context of this work, however, the central idea 153 is that the system (5) must be assembled and solved for many 154 different parameters. To this end, we introduce the parametric 155 solution to the problem, u(x; µ), where the (scalar) problem 156 parameters µ = (µ 1 ,...,µ n µ ) are defined over the param-157 eter domains I µ = I µ 1 ×···×I µ nµ .

158
The pivotal idea of the PGD method is to attain u(x; µ) 159 as the solution to a problem posed on the higher-dimensional 160 domain × I µ , the spatial semi-discretization of which can 161 be written as: The general PGD strategy to obtaining this solution is 165 to formulate a higher-dimensional weak form problem 166 corresponding to (2), and then to discretize this higher-167 dimensional problem in space and in the parametric dimen-168 sions; see, e.g., [9,10] for the fundamentals of PGD. An 169 essential aspect of the PGD framework is that in order to 170 efficiently compute the parametric solution, a separable form 171 of the weak form problem (or its discrete version) must be 172 available. With respect to the spatially discretized system (5) 173 this means that the stiffness matrix and force vector should 174 be of the form, where n k and n f denote the total number of terms needed 179 to represent the parametric stiffness matrix and parametric 180 force vector, respectively. Note that when these affine repre-181 sentations are not available, it is possible to construct affine 182   x = 2 l c X for X ≤ 0.5,

214
The Jacobian of this mapping follows as: The inverse of this Jacobian can be obtained analytically and 217 allows for an exact separable representation as the sum of 218 products of matrices that do not depend on the parameter l c 219 and functions that depend only on that parameter: A separable form of the determinant of the Jacobian can 222 similarly be obtained: The matrix and vector components in Eq. (6) can now be 225 transformed via the mapping M(X, l c ) into equivalent inte-226 grals over the reference domain as  The corresponding stiffness matrices are obtained as: Similarly, n f = 2 parametric shape functions are found for 279 the force vector: The corresponding vector components are found as: The system composed of these separable forms for the stiff-288 ness matrix and force vector assumes the canonical form (7). 289

290
(PGD) method 291 The parametric problem (7)    the Euclidean norms û i and ĝ i j , respectively, such that 318 the modal amplitudes, β i , are given by: 320 We employ the PGD solver algorithm as presented in Ref. resentation. This is not generally the case, as we will 328 discuss, for example, in the stochastic test case con-329 sidered in Sect. 7. In situations where the linear system 330 cannot be separated analytically, it is often replaced by 331 a separable approximation (e.g., [30,31] . An overview of these techniques can be 340 found in, e.g., Ref. [21]. It is noted that in the case of 341 high-dimensional parameter domains, the computation 342 of separable forms can be computationally demanding.

343
-A greedy algorithm [1,8] is used to sequentially compute 344 the terms to the PGD approximationû pgd in Eq. (18). 345 Given the PGD approximation with n pgd − 1 terms, here 346 denoted by an enrichment termû n pgd n µ j=1 G n pgd j is computed as to 349 obtain the PGD approximation with n pgd terms: Each enrichment term is computed one at a time, con-352 structing the summation progressively until the conver-353 gence criterion is met with a user-defined tolerance of ǫ glob . Each step in 356 the greedy algorithm, i.e., computing each of the enrich-357 ment terms, involves the computation of the enrichment 358 modes in space,û i in discrete form, and in the parameter 359 spaces, G i j (µ j ). We herein compute these enrichments 360 iteratively using an alternate direction solver, which is 361 discussed in detail below.

362
-A nalternating direction solution strategy [9]isusedto 363 compute the enrichment termsû n pgd n µ j=1 G n pgd j .Lever-364 aging the separable forms, in this alternating direction 365 strategy the spatial and parametric directions are treated 366 sequentially as to reduce the higher-dimensional para-367 metric problem to a series of low dimensional problems. 368 This iterative process is repeated until a fixed point is 369 reached within a defined tolerance. For the explanation 370 of this alternating direction strategy we will consider 371 n µ = 1 with the fracture length µ 1 = l c as the only 372 parameter.

373
For the alternate direction solution strategy, the paramet-374 ric problem (7) is considered in its weighted residual 375 form: (24)

379
The unknowns in this system are the spatial and paramet-380 ric enrichment modes,û n pgd and G n pgd l c (l c ), respectively.

381
The corresponding test functions are defined as:

384
In the alternate direction strategy, the system (24)i s 385 solved per spatial or parametric dimension:

386
-Given an approximation (or initial guess) for the 387 parametric enrichment mode G n pgd l c , the system (24) 388 reduces to the linear system:

392
Using the separable forms for the stiffness matrix 393 and force vector in equation (9), this system can be 394 rewritten as

408
Equivalently, it holds that for each fracture length l c from which the parametric enrichment mode follows 411 directly as: Substitution of the separable forms for the stiffness 414 matrix and force vector then finally yields: This expression for the parametric enrichment mode 419 can be evaluated quickly by virtue of the fact that 420 the dimensions are separated in the sense that it is 421 not required to reassemble the finite element system 422 for each fracture length. The parametric enrichment 423 mode is represented discretely by projection onto the 424 parametric basis in Eq. (19). Since this discretiza-425 tion pertains to a linear finite element basis, the 426 coefficientsĝ n pgd l c can be computed by evaluation of 427 Eq. (31) in the parametric nodes.

428
The above alternate direction steps are repeated until 429 the relative difference between two successive steps is 430 smaller than a prescribed tolerance, ǫ local , Before considering the application of the PGD framework to 438 fracture problems, in this section we first present a numerical 439 is not exact, as we will consider, for example, in the context dependence of the PGD approximation on the modes is con-477 sidered below (Fig. 4).

478
To study the approximation behavior of the PGD expan-479 sion, we consider the relative energy error with respect to the 480 original finite element solution: In contrast to (34), this error measure provides one scalar 493 error value for the complete parametric solution and has no 494 dependency on l c . Figure 5 displays both error measures for 495 various spatial mesh sizes, h, and a fixed parametric mesh 496 size h l c ≈ 0.015 m. The parameter dependent error (34)dis-497 played in Fig. 5a conveys that for a certain mesh size, the error 498 in the PGD solution is dependent on the crack length. The 499 reason for this is that the uniformity of the mesh in the phys-500 ical domain is affected by the parameter-dependent mapping 501 function (10), which in general causes the error to increase 502 when the crack tip position deviates from l c /H x = 0.5 (i.c., 503 l c = 2) provided that the mesh resolution is of sufficient 504 accuracy. The error e pgd (l c ) is especially significant at the 505 boundaries of the parameter domain, I lc , because at those 506 points the non-uniformity caused by the mapping onto the 507 physical domain (see Fig. 3) is largest.

508
When we compute the mean of the error e pgd (l c ) over 509 the complete parameter domain, i.e., error measure (35), we 510 observe from Fig. 5b that this mean energy error is essentially 511 independent of the mesh size for the finer meshes (h 0.25). 512 This conveys that for these meshes the studied error is dom-513 inated by the PGD approximation, which is expected, as we 514 compare the PGD solution with the FE solution on the same 515 mesh.

516
To study the mesh size contribution to the PGD approxi-517 mation error, in Fig. 6 we display the mean L 2 error between 518 a PGD approximation u pgd (x; l c ) computed with mesh size h 519 and a PGD approximation, u ⋆ pgd (x; l c ), with a high resolution 520 mesh with h ⋆ = 0.03125:    Figure 8 shows that both the parameter-dependent energy 548 error (34) and mean energy error (35) are virtually indepen-549 dent of the parametric mesh size even on parametric meshes 550 as coarse as h l c = 0.125 m (8 elements). This conveys that, 551 in the setting considered here, the accuracy is governed by 552 the number of PGD modes rather than by the resolution of 553 the parametric mesh.     590 We consider Griffith's model [16] for crack propagation in 591 brittle materials. The conceptual idea of this model is that a 592 fracture will propagate if the energy stored in the material is 593 sufficiently large to overcome the fracture energy associated 594 with the creation of new fracture surface. For linear elastic 595 materials an equivalent interpretation of this energy-based 596 model is provided through the concept of stress intensity 597 factors [5]. In the context of the PGD framework we find the 598 energy perspective most suitable, as it provides the possibility 599 to evaluate the propagation criterion directly based on the 600 parametric solution (37).

601
For a fracture in a given configuration, i.e., with a cer-602 tain length l c and a given load scale λ, it can be determined 603 whether or not the fracture will propagate by evaluation of 604 the energy release rate. To derive the PGD form of the energy 605 release rate, we consider the energy of the system: The energy release rate is then defined as :

639
As a means to assess the PGD approximation of the energy 640 release rate, we study the stress intensity factor for a given 641 fracture length l c , and various ratio's of horizontal and ver-642 tical specimen dimensions, H x and H y , respectively. The 643 results presented in this section consider the parameters 644 H x and H y as additional parameters in the PGD expan-645 sion. The separable forms based on these parameters can 646 be obtained without special treatment, and are omitted here 647 for the sake of brevity. The stress intensity factor is defined 648 as  Figure 9 shows the dimensionless stress intensity factors 654 K 1 /K 0 for various parameter configurations, i.e., different 655 l c /H x and H x /H y (see Ref. [26] for a benchmark result). 656 Note that the plotted factors are non-dimensionalized using 657 K 0 = (λt · n) √ πl c , where λt · n gives the magnitude of 658 the applied tensile traction. Figure 9 compares the PGD 659 results based on the settings mentioned in Table 1 for a 660 mesh size h = 0.0625 m. However, note that this plot of 661 non-dimensional stress intensity factors is independent of 662 the input values, i.e., even for different values of geome-663 try and load, similar curves for K 1 /K 0 are obtained. This 664 figure conveys that for the given PGD settings, the stress 665 intensity factor can be computed accurately using the PGD 666 expansion (37). While each point in Fig. 9 would typically 667 represent a finite element simulation in the traditional FEM 668 setting, in the PGD case these are all mere evaluations of the 669 expansion. 670 Fig. 9 Dimensionless stress intensity factors K 1 /K 0 for various crack lengths in specimens of various dimensions loaded in tension. The solid lines represent the results computed through the PGD framework, while the markers indicate the reference values reported in Ref. [26] Now that we have established that the PGD expansion accu-672 rately approximates the stress intensity factor, we will here 3. The unstable propagation region: The energy landscape is plotted in Fig. 10a

699
For a particular load scale, until the critical point is reached 700 the crack is stable (green region in Fig. 10a), and beyond the 701 maximum point the crack is unstable (red region in Fig. 10a).

702
The critical energy states are connected in the form of a curve 703 which gives the critical load value for each fracture length.

704
This curve can be identified in Fig. 10a 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 pro-717 cess is by plotting the load versus the average displacement 718 of the loading boundary, which is depicted in Fig. 10bf o r 719 a initial crack length of l 0 c = 2.495 m. Note that the elastic 720 loading branch (label I. in Fig. 10) corresponds to the region 721 where the crack is stable, i.e, the force varies with ∂E ∂l c < 0. 722 The resultant force at which the crack becomes unstable, i.e., 723 when ∂E ∂l c = 0, is defined as the critical loading force, F c .This 724 corresponds to the maximum force in Fig. 10b. This critical 725 loading force is related to the dimensionless stress intensity 726 factors of Fig. 9 by: The softening branch (label II. in Fig. 10)  which is defined as

772
In the context of the stochastic analysis considered here, 773 we use the PGD framework to compute the parametric solu-774 tion with respect to the fracture length, external load, and 775 with the random variablesz α that parametrize the random 776 Young's modulus field: (47) 778 A prerequisite to apply our framework is to express the stiff-779 ness matrix and force vector also in this separated format. 780 The separable forms of the stiffness matrix and force vec-781 tor required here cannot be obtained in an analytical way 782 like in Sects. 3 and 6. Therefore, in Sect. 7.1 we first discuss 783 how the random heterogeneities, which are parametrized by 784 the random variablesz, can be expressed in a separable form 785 for the stiffness matrix numerically. Furthermore, in Sect. 7.2 786 we outline the computational procedure for a sampling-based 787 stochastic analysis based on the Monte-Carlo method. This 788 stochastic analysis is highly efficient as it leverages the PGD 789 approximation to quickly compute critical force values for 790 realizations of the heterogeneous field of elastic properties. 791 Numerical results for the stochastic test case are presented 792 in Sect. 7.3. where the constant tensor D depends on the Poisson ratio and 802 on the assumed plane strain state. Since the elasticity tensor 803 is evaluated over the reference domain, the KL modes {R • 804 M} n kl α=1 are pulled back to the reference configuration using 805 the geometric mapping function (10). Since this mapping 806 function is dependent on the fracture length parameter l c ,the 807 random elasticity tensor (48) also becomes dependent on the 808 fracture length. elasticity tensor set to C = µ E D, which we denote by  849 where σ (α,β) is the β-th singular value for KL mode α, and 850 whereĥ (α,β) andm (α,β) are the corresponding spatial and 851 parametric modal vectors, respectively. For reasons of effi-852 ciency this singular-value decomposition is truncated to a 853 number of terms, n svd , that is significantly smaller than the 854 total system size. Substitution of this decomposition into Eq. 855 (52) then yields the singular-value decomposition for the KL 856 modal functions, where the modal functions are defined as (56) 872 The components of the matrices K (α,β) (l c ) are given by: Since the spatial modes, h (α,β) (X), are independent of the 877 parameter l c , the matrices K (α,β) can be separated analo-878 gously to the Eqs.   where n sample is the Monte-Carlo sample size.

929
In a typical FE-based Monte Carlo simulation, evaluation 930 of the critical loads is computationally demanding, which 931 practically restricts the sample sizes that can be considered. 932 Therefore, in such cases, a sample size is selected that strikes 933 an adequate balance between the confidence level of the 934 attained statistical moments and the required computational 935 effort. In the PGD setting considered here, the computational 936 effort involved in determining the critical force for a given 937 realization of the random field is negligible compared to the 938 corresponding full finite element simulation. This allows for 939 the consideration of sample sizes that are orders of magni-940 tude larger than those that could be considered using direct 941 FE analysis, which in turn enables the computation of the sta-942 tistical moments with confidence levels that are practically 943 beyond the reach of direct FE analyses. Evidently, the selec-944 tion of the sample size should be based on a trade-off between 945 the error in the PGD approximation and the confidence level 946 of the Monte Carlo method. We consider the same numerical experiment as introduced in 950 Sect. 6.2 (see Table 1), but now with a random field of elastic 951 properties. For the random field (45) we set the mean to µ E = 952 2 GPa and the standard deviation as σ E = 0.2 GPa (a coef-953 ficient of variation of 10%). We consider moderate spatial 954 fluctuations in the random field by selecting the correlation 955 length in Eq. (46)a sl E = 1.5 H x = 6 m. The parameter 956 domain for the load scale is taken as I λ =[6.25, 62.5].

957
We consider a Karhunen-Loève discretization consisting 958 of n kl = 3 modes, which are shown in Fig. 11.I nF i g .12 959 we show two realizations of the KL expansion, as well as a 960 sampling-based reconstruction of the auto-correlation func-961 tion (46). On account of the low spatial frequency of the 962 variations, the KL expansion with only 3 terms is observed 963 to already appropriately reproduce the auto-correlation func-964 tion.

965
Using the tolerances specified in Table 1, the PGD solution 966 (47) is truncated at n pgd = 27 terms. The various compo-967 nents of the PGD solution are displayed in Fig. 13.Fromthe 968 modal amplitudes it can be observed that the PGD approx-969 imation based on 27 terms approximates the finite element 970   Figures 14 and 15 show that the average critical load 989 bearing capacity decreases with an increase in crack length, 990 while a decrease in the standard deviation is observed. The 991 deterministic result is plotted for reference, from which it is 992 observed that the computed mean is slightly smaller than the 993 deterministic value. The observed results from the Monte 994 Carlo simulation are in good agreement with perturbation 995 analysis results (see [17] for an overview) based on the ana-996 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 997 to be expected on account of the considered low spatial fre-998 quency of the random input.

999
The Monte Carlo analysis allows us to inspect which real-1000 izations of the input lead to a certain response in terms of 1001 the fracture load. Figure 16 shows three interesting realiza-  c. The realization with the smallest fracture load of 66.6 MN 1013 corresponds to a Young's modulus field which is very low 1014 throughout the specimen (on average approximately 25% 1015 lower than its mean value), and particularly near the tip. 1016 In the context of the PGD approach employed in this work it 1017 is noted that, in order to inspect these realizations, only the 1018 parameters corresponding to the realization (random vari-1019 able realizations) have to be stored. The input and output 1020 corresponding to these parameters is generated through post-1021 Fig. 13 The seven components of the u pgd (l c ,z 1 ,z 2 ,z 3 ,λ) solution for n pgd = 27. Only a selection of modes is shown for conciseness. Note that all plotted functions are normalized   The suitability and performance of the proposed frame-1036 work is demonstrated using a series of numerical test cases, 1037 starting with a convergence study for the parametric decom-1038 position. This study conveys that the introduced geometric 1039 mapping function for the fracture parameter behaves in 1040 accordance with the well-understood behavior of the PGD 1041 framework. The PGD fracture framework is further demon-1042 strated using two propagating fracture test cases.

1043
In the first test case it is demonstrated how Griffith's prop-1044 agation criterion can be evaluated efficiently using the PGD 1045 approximation. The representation of the fracture length in 1046 the PGD solution enables the straightforward computation 1047 of the energy release rate, which is in contrast with standard 1048 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 1049 numerical techniques for the evaluation of the correspond-1050 ing shape derivative. besides the length, would require the fracture angle as an 1081 additional parameter. In general, however, representing 1082 more complex fracture geometries will rapidly increase 1083 the number of parameters, which is detrimental to the 1084 performance of the PGD framework. This is particularly 1085 the case when one opts to consider a piecewise repre-1086 sentation of fractures, which is natural to finite element 1087 methods.

1088
-For more complex fracture patterns, constructing a suit-1089 able geometric mapping function will be considerably 1090 more challenging than in the prototypical benchmark 1091 considered in this work. Constructing a mapping analyti-1092 cally is very restrictive, but it is very well imaginable that 1093 one can construct discrete mapping operators (mapping 1094 nodal reference coordinates to nodal physical coordi-1095 nates). Such more advanced mappings -the construction 1096 of which evidently warrants further investigation -will, 1097 however, pose several difficulties. For example, the ana-1098 lytical separation of the system of equations as obtained 1099 in this work will not be generally obtainable, which hence 1100 requires the consideration of potentially computation-1101 ally demanding approximations for the separable forms. 1102 Moreover, an open research question remains how to deal 1103 with fractures with changing topology (e.g., branching, 1104 merging), as topological changes can in general not be 1105 captured by the proposed mapping technique.

1106
These complications when extending to more complex frac-1107 tures are evidently very serious. Although future research 1108 developments can ameliorate some of these difficulties, 1109 obtaining PGD approximations that are able to accurately 1110