A Manifold Learning Approach to Accelerate Phase Field Fracture Simulations in the Representative Volume Element
SSN Applied Sciences manuscript No. (will be inserted by the editor)
A Manifold Learning Approach to Accelerate Phase Field FractureSimulations in the Representative Volume Element
Yangyuanchen Liu · Kexin Weng · Yongxing Shen
Received: date / Accepted: date
Abstract
The multiscale simulation of heterogeneous ma-terials is a popular and important subject in solid mechanicsand materials science due to the wide application of compos-ite materials. However, the classical FE (finite element )scheme can be costly, especially when the microproblem isnonlinear. In this paper, we consider the case when the mi-croproblem is the phase field formulation for fracture. Weadopt the locally linear embedding (LLE) manifold learningapproach, a method for non-linear dimension reduction, toextract the manifold that contains a collection of phase-field-represented initial microcrack patterns in the representativevolume element (RVE). Then the output data correspondingto any other microcrack pattern, e.g., the evolved phase fieldat a fixed load, can be accurately reconstructed using thelearned manifold with minimum computation. The methodhas two features: a minimum number of parameters for thescheme, and an input-specific error bar. The latter featureenables an adaptive strategy for any new input on whetherto use the proposed, less expensive reconstruction, or to usean accurate but costly high-fidelity computation instead. Keywords
Multiscale simulation · Manifold learning · Locally linear embedding · Phase field for fracture
This work is supported by the Guangdong Province Key Area R&DProgram, grant No. 2019B010940001, by the Shanghai Natural Sci-ence Foundation, grant No. 19ZR1424200, and by the National NaturalScience Foundation of China, grant No. 11972227.Y. Liu · K. Weng · Y. ShenUniversity of Michigan – Shanghai Jiao Tong University Joint InstituteShanghai Jiao Tong University800 Dongchuan Road, Shanghai, 200240, ChinaE-mail: [email protected]. WengCollege of EngineeringUniversity of MichiganAnn Arbor, MI 48109-2102, USA
Heterogeneous materials such as composites have been widelyapplied in various industries such as aircraft and automobilemanufacturing. The multiscale simulation of heterogeneousmaterials is therefore a crucial task in computational me-chanics.Such simulation is usually facilitated by the classicalFE scheme [1], as illustrated in Fig. 1. In a typical FE scheme, the finite element method is applied at the microscaleand the macroscale concurrently, and hence the name. Moreprecisely, at the macroscale, the entire composite part is dis-cretized into continuum finite elements, each of which hasseveral Gauss quadrature points for numerical integration.For each Gauss quadrature point, the effective constitutivebehavior for the macroscale is obtained through a homog-enization process via a finite element analysis at the mi-croscale. The computational domain at the microscale is calleda representatixve volume element (RVE). Take the mechan-ical simulation for a fiber-reinforced composite as an exam-ple, a typical RVE consists of a fiber and the surroundingmatrix [2], possibly with defects such as cracks. Normallythe desired effective responses include the stress tensor andthe elasticity tensor, and the simplest way of homogeniza-tion is by volume averaging.Among available numerical methods for the analysis atthe RVE with crack propagation, the phase field approach tofracture [3], also known as the regularized variational the-ory for fracture, shows clear advantages. This approach isbuilt on Griffiths theory for brittle fracture [4]. The key ideais to use a scalar field, called phase field, to represent thecrack path, instead of incorporating the explicit geometryof the crack path in the computational domain. The advan-tages include obviating the need for explicitly tracking thecrack path geometry, and the ability to predict crack nucle-ation and bifurcation without extra criterion. The method a r X i v : . [ c s . C E ] J u l Liu, Weng, and Shen
Microscale stress 𝝈𝝈 Boundary condition EquilibriumCompatibilityConstitutive equationLocalization Periodic medium
Macroscopic strain �𝜺𝜺 ⋅ Homogenization
Macroscopicstress �𝝈𝝈 = 𝝈𝝈 Fig. 1
Flowchart illustrating the FE scheme. has since been applied to fracture modeling in Euler-Bernoullibeams [5], thin shells [6], composite materials [7,8], cement-based materials [9], layered structures [10], and CO fractur-ing [11].However, solving the equations arising from the phasefield method for fracture can be costly. Since the strain en-ergy functional to minimize in this approach is not con-vex, the required number of iterations for convergence isnot known a priori . The RVE analysis is, of course, no ex-ception. Many efforts have been devoted to accelerating thephase field fracture solution procedure. Heister et al. [12]and Li et al. [13] constructed mesh adaptivity approachesfor the problem. Ziaei-Rad and Shen [14] developed a mas-sively parallel algorithm for the phase field approach withtime adaptivity. Gerasimov and De Lorenzis [15] proposeda line search procedure for the monolithic scheme to over-come the iterative convergence issues of non-convex mini-mization. Wick [16,17] developed modified Newton-Raphsonschemes for fully monolithic quais-static brittle phase fieldfracture propagation. Farrell and Maurini [18] reformulatedthe staggered algorithm of the phase field analysis as a non-linear Gauss-Seidel iteration and employed over-relaxationto accelerate convergence. Wu et al. [19] developed a quasi-Newton monolithic methodwith the Brodyen-Fletcher-Goldfarb-Shanno (BFGS) algorithm. Kopaniˇc´akov´a and Krause [20]proposed a trust region method with application to mono-lithic phase-field fracture models.We aim to accelerate the multiscale simulation from an-other perspective. In fact, in many cases, the RVEs are sim-ilar within the same multiscale analysis. This similarity canbe exploited to accelerate computation, for example, via man-ifold learning.In the machine learning context, manifold learning isemployed to extract the manifold that represents high-dimensionaldata points and to perform data reconstruction with a min-imum amount of computation. Manifold learning has beenwidely applied to multiscale analysis [21,22,23,24], see alsothe review by Matouˇs et al. [25]. An instance of manifold learning techniques is locally linear embedding (LLE). Pro-posed by Roweis and Saul [26], LLE is an unsupervisedlearning algorithm that computes low-dimensional, topology-preserving embeddings of high-dimensional data points. Asan instance of kernel principal component analysis (kernelPCA), LLE has many attractive properties. For example, thelocal geometry of high-dimensional data is preserved in thelow-dimensional manifold. LLE is particularly suitable forproblems with a large amount of similar high-dimensionaldata.However, LLE assumes that the data all reside on a sin-gle continuous manifold [27], which poses certain restric-tions on the application. For example, in image-based simu-lations [28], each RVE is represented as a vector containing,e.g., pixel values. In this case, if the dimension of this vec-tor varies between RVEs, the nonuniform data structure willmake LLE training and interpolation impossible. This is be-cause the neighborhood finding and interpolation operationsof the LLE algorithm requires that the linear combination ofdata points to be well defined.Despite such restrictions, the advantages of LLE make itideal for random RVE computation and computational ho-mogenization [28,29,30] for multiscale analysis of hetero-geneous materials.Inspired by [28] for heat conduction problems, for theproblem of multiscale fracture simulation at hand, we aim tolearn a manifold that contains a collection of similar crackedRVEs, and to efficiently compute any desired output de-pendent on such microstructure using LLE reconstructions.Concretely speaking, the input is chosen as the phase fieldpattern at the beginning of a certain time step (termed “ini-tial phase field” for short), and the output can be the phasefield at the end of the time step – so as to make a closed loopfor the analysis of the next step – or any other derived quan-tity from such phase field solution such as the homogenizedstress. In the discrete picture, we construct a finite elementmesh to describe the RVE, interpolate the phase field for thecrack pattern using the finite element basis functions, andvectorize the description of the initial crack pattern of eachRVE using the nodal values of the phase field. The desiredoutput is the phase field solution corresponding to a certainboundary condition.Compared with recent contributions on applying machinelearning techniques, neural networks in particular, for con-stitutive modeling [31,32,33,34,35,36] and similar compu-tations for RVEs [37,38,39,40,41], the adopted method pos-sesses the following features.First, the number of hyperparameters is minimal: onlythe size of the neighborhood and the number of reduceddimensions need to be input by the user. The selection ofsuch hyperparameters is determined by a systematic cross-validation approach. Manifold Learning Approach to Accelerate Phase Field Fracture Simulations in the Representative Volume Element 3
Second, there is no limit on the dimension of the desiredoutput, as long as it is a continuous functional of the mi-crostructure, while a typical neural network would have oneset of thresholds and weights per scalar output.Finally, for any new input, the uncertainty (“error bar”)for the reconstructed output can be obtained, as a strong cor-relation is observed between the reconstruction error and aparameter solely dependent on the input information. In thiscase, the parameter is the distance from the new input tothe learned data manifold. This last feature enables a cri-terion to be developed to assess the reconstruction error apriori ; in other words, a criterion to decide whether to usethe reconstruction which is less expensive, or resort to thehigh-fidelity computation which is more accurate. This alsoserves as an indicator of whether the collection of inputsshould be augmented with the new input in question, in agreedy sampling fashion, should some kind of adaptivity isto be implemented.However, it is still worth noting that, just like many othermachine learning techniques, the LLE approach requires enoughdata points to guarantee the accuracy of predictions. Hencethe training set should be dense and large enough. More-over, as inherited from the general LLE technique, the pro-posed approach requires the data structure to be homoge-neous, making the distance function and linear combinationbetween data points well-defined. Finally, the output shouldcontinuously depend on the input data, which is also a nec-essary condition for a well-posed problem anyway.The content of this paper is structured as follows. InSection 2, the FE scheme and phase field method are in-troduced. In Section 3, the manifold learning and LLE tech-niques are explained in detail. In Section 4, numerical imple-mentations and results are illustrated with error assessments.Finally, in Section 5, a summary of the proposed computa-tional strategy is presented. Scheme Applied to Composite Fracture
In this section, we introduce the FE scheme in the mul-tiscale fracture simulation of a fiber-reinforced composite.The FE is a two-scale modeling scheme which applies FEdiscretizations at both macro and micro scales, the formertaking input from the latter through the analysis of the RVE.In our case, as shown in Fig. 2, the RVE is composed ofa strong fiber in the center with a weaker matrix. We aim toperform the fracture simulation of the cracked RVE at themicroscale. Once the local behavior is determined, the over-all macroscopic response of the RVE can be obtained usingany well-established homogenization theory and be used forthe macroscopic simulation.For simplicity, we only consider the microcrack evolu-tion in the matrix and ignore all other defects, such as crackson the interface (debonding) and in the fiber, see Fig. 3. Macroscopic Microscopic RVEMicrocrack
Fig. 2
Modeling a macroscopic composite as a collection of RVEs.
𝑆𝑆ℬ 𝑠𝑠 𝑑𝑑 ( 𝒙𝒙 , 𝑡𝑡 ) Fig. 3
The simplified RVE to be analyzed in this work. In this RVEthere is a strong fiber inside a weaker matrix. The only allowed formof failure is matrix cracking.
Phase Field Approach for RVE Cracking.
Among many cracksimulation methods, we adopt the phase field method to sim-ulate the microcrack evolution in RVE. The phase field mod-eling of brittle fracture has shown its advantages on simulat-ing complex fracture process, such as obviation of remesh-ing, see [3,42,43]. The phase field approach of fracture isbased on the variational energy formulation proposed by[44], which can be considered as a generalization of Grif-fiths theory [4].As shown in Fig. 4(b), the phase field method uses a dif-fuse field d to represent the cracked microstructures where d = d = Liu, Weng, and Shen vantageous in terms of data structure for the manifold learn-ing approach, as each cracked microstructure can be uni-formly represented as a vector consisting of the nodes’ phasefield values. This feature is favorable in the manifold learn-ing process introduced in Section 3, as we can adopt a datastructure for the inputs (and outputs) as vectors of the samelength. (a) 𝑙𝑙𝑙𝑙𝑑𝑑 ( 𝒙𝒙 , 𝑡𝑡 ) (b) (c) Fig. 4
Representations of a unit cracked microstructures: (a) discretecrack model; (b) phase field corresponding to (a); (c) pixel represen-tation of the phase field model with a structured quadrilateral mesh ofwith h / l = . Fig. 4(c) and Fig. 5(a) show the pixel representationsof microcracks with the phase field approach. At the firstsight, there are at least two possible alternatives to trans-late a cracked microstructure into a numerical representa-tion: (1) using the characteristic function of the cracks, i.e.,1 for the crack and 0 otherwise, as shown in Fig. 5(b), (2) us-ing the distance function to the cracks, as shown in Fig. 5(c).Considering that we will need to quantify the “distance” ofsuch microstructures, both alternatives present severe draw-backs: method (1) would not be able to tell the distance ofnon-overlapping cracks, while method (2) would weight toomuch on the difference of crack pattern pairs in areas faraway from the cracks.The adopted variant of the phase field formulation isas follows. In a plane strain setting, let B = ( − L , L ) bethe area initially occupied by the RVE. Within the RVE, let S ⊂⊂ B be the fiber, and B s = B \ S be the matrix, seeFig. 3. In the absence of body force and traction boundarycondition, the phase field formulation for the RVE is [45] Π l [ u , d ] = (cid:90) B s Ψ [ ε , d ] d B + (cid:90) S Ψ ( ε ) d B + g c (cid:90) B s (cid:18) d l + l | ∇ d | (cid:19) d B , (1)where the arguments u ∈ H ( B , R ) and d ∈ H ( B s ) arethe displacement field and the phase field, respectively, andthe strain tensor is defined as ε = ( ∇ u + ∇ u T ) /
2. Here weset the convention for the phase field d as d = d = ( λ , µ ) and ( λ , µ ) be the Lam´e constants of the matrix and the fiber, respec-tively, then the strain energy density for the fiber is given RVE (1)
RVE (2) RVE (3) (a) Phase field (chosen) RVE (1) RVE (3)RVE (2) (b) Characteristic function (not recommended)
RVE (2)
RVE (1) RVE (3) (c) Distance function (not recommended)
Fig. 5
Numerical representations of cracked microstructures with 5 × far away fromthe cracks. by Ψ ( ε ) = λ ( tr ε ) + µ ε : ε , while that for the matrix also depends on d , for which weadopt the formulation proposed by Amor et al. [42]. Thismodel splits the strain energy density Ψ into volumetric anddeviatoric parts: Ψ ( ε , d ) = g ( d ) Ψ + ( ε ) + Ψ − ( ε ) , Manifold Learning Approach to Accelerate Phase Field Fracture Simulations in the Representative Volume Element 5 where Ψ + ( ε ) = K (cid:104) tr ε (cid:105) + + µ (cid:107) dev ε (cid:107) , (2a) Ψ − ( ε ) = K (cid:104) tr ε (cid:105) − , (2b) σ ( ε , d ) = g ( d ) (cid:0) K (cid:104) tr ε (cid:105) + + µ dev ε (cid:1) + K (cid:104) tr ε (cid:105) − , (2c)where K = λ + µ / ε : = ε − ( / )( tr ε ) , (cid:104) a (cid:105) ± : = ( a ± | a | ) / g ( d ) = ( − d ) + k , where k is a small positive number.The positive numbers g c and l are the energy release rateof crack propagation and the regularization length scale, re-spectively.The strong form of the governing equations, except thedisplacement boundary condition at ∂ B , readdiv σ = , in B s ∪ S , (3a) σ = ∂Ψ∂ ε , in B s , (3b) σ = ∂Ψ ∂ ε , in S , (3c) ∂Ψ∂ d g c (cid:18) dl − l ∆ d (cid:19) = , in B s , (3d) σ · n (cid:12)(cid:12) B s = σ · n (cid:12)(cid:12) S on ∂ S (3e) u (cid:12)(cid:12) B s = u (cid:12)(cid:12) S on ∂ S (3f) ∇ d · n = ∂ B s , (3g)where n denotes the outward unit normal vector of ∂ S or ∂ B .The general quasi-static calculation for each load stepof the microcrack evolution is shown in Fig. 6: the inputsare the crack configuration (represented by a phase field) attime t and the boundary conditions for u and d at the nexttime step t + ∆ t , and the output is the updated phase field at t + ∆ t . Here t represents a time-like variable to indicate theprocess of load increment, and likewise t + ∆ t . 𝑆𝑆ℬ 𝑠𝑠 𝑑𝑑 ( 𝒙𝒙 , 𝑡𝑡 ) 𝑆𝑆𝑑𝑑 ( 𝒙𝒙 , 𝑡𝑡 + Δ𝑑𝑑 ) 𝒖𝒖 𝒙𝒙 , 𝑡𝑡 + Δ𝑡𝑡 � 𝜕𝜕ℬ (a) (b) (c)
ℬ ℬ 𝑠𝑠 Fig. 6 (a) RVE with micro cracks; (b) the boundary conditions of RVE;(c) RVE with evolved micro cracks
For simplicity, we fix the following boundary conditionson ∂ B and focus on the effect of the crack path at t on itsupdated counterpart at t + ∆ t . Let ε ∈ R × be the imposed macroscopic strain tensor, then the boundary conditions areset to be u = ε · x , on ∂ B . (3h) The FE scheme introduced in Section 2 requires an un-predictable number of iterations for convergence due to thenon-convexity of the functional Π l . In order to reduce com-putational cost, we adopt the so-called manifold learningmethod. The manifold learning scheme uses techniques tra-ditionally designed for machine learning purposes to extractthe manifold that represents high-dimensional data pointsand perform reconstruction with minimum amount of com-putation [28,30]. The main idea is to generate enough inputsand pre-compute their outputs offline, in this case the phasefields at t and t + ∆ t , respectively, then provides the desiredoutput for any input by reconstruction.In this section, we will elaborate on the manifold learn-ing approach and the LLE technique [26], specialized to theproblem stated in Section 2. In particular, as we fix the loadshown in Fig. 6(b), the only input to consider is the initialcrack path (i.e. the initial phase field) [Fig. 6(a)], and theoutput is the evolved phase field [Fig. 6(c)] upon equilib-rium.3.1 Locally Linear EmbeddingLocally linear embedding (LLE), proposed by Roweis andSaul [26], is an unsupervised learning algorithm that com-putes low-dimensional, topology-preserving embeddings ofhigh-dimensional data points. LLE is an instance of kernelprincipal component analysis (kernel PCA), which handlesnonlinear dimensionality reduction [46]. As illustrated inFig. 7, LLE maps high-dimensional data into a single globalcoordinate system of lower dimensionality.In this paper, we use LLE to accelerate the computationof the phase field. The main idea is that from the offline cal-culation of enough cracked microstructures, we will be ableto reconstruct crack evolution due to various initial crackpatterns with minimal computation online.The specific process of LLE is as follows. Suppose thatthere are N input data points X i ∈ R D where i = , ..., N ,each X i containing the phase field values representing a spe-cific cracked microstructure. According to [26], under theassumption that all inputs are on the same manifold, we canlinearly reconstruct each data point X i by its k ( (cid:28) N ) near-est neighbors, say X i = ∑ j ∈ S i W i j X j , (4) Liu, Weng, and Shen (a) (b) (c)
Fig. 7
The illustration of locally linear embedding. (a) A two-dimensional manifold; (b) the three-dimensional data points sampledfrom (a), colored according to the z -coordinates; (c) the data pointsafter dimensionality reduction by LLE. where W i j are the weights to be determined and S i representsthe set of the k nearest neighbors of X i in the l -norm.To compute these weights W i j , we minimize the costfunction which measures the reconstruction errors: F ( W ) = N ∑ i = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X i − ∑ j ∈ S i W i j X j (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) . (5)The minimization of F ( W ) is subjected to two constraints:(i) each data point X i is reconstructed only from its neigh-bors: W i j = X j / ∈ S i . (ii) the rows of the weight matrix sum to 1: ∑ j ∈ S i W i j = i = , ..., N . An important featureis, for any data point, the weights are invariant to rotation,rescaling and translation of that data point with respect to itsneighbors [26].Now we suppose that all data points are mapped into alower dimensional embedding space (manifold) of dimen-sion L , L (cid:28) D . The reconstruction weights W i j remainunchanged in such transformation. Therefore, each high di-mensional data point X i is mapped to a low dimensionalvector Y i representing coordinates on the manifold. We com-pute Y : = { Y i } by minimizing the embedding cost function G ( Y ) = N ∑ i = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) Y i − ∑ j ∈ S i W mi Y j (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , Y : = { Y i } . (6)During this minimization, the weights W i j are fixed. To fullydetermine { Y i } , certain constraints have to be imposed sothat the solution is unique [26]. The resulting constrainedminimization problem can be solved via an N × N eigen-value problem.3.2 Training and Output ReconstructionAs previously discussed, the offline procedure of this mani-fold learning scheme consists of two stages: (1) dataset gen-eration with the phase field analysis for the RVE, (2) datamanifold construction with LLE. Then for any given phasefield under the same load, the online reconstruction proce-dure readily delivers the phase field evolution.To generate the training data, we subject a series of RVEswith an initial crack at various locations to the unilateral ten-sion test. The configuration and mesh with an initial phasefield are shown in Fig. 8. The mesh shown in Fig. 8(b) con-tains D nodes, so every input data point X i as well as thecorresponding output data point Z i is a column vector with D phase field values. 𝒖𝒖 𝐷𝐷 = �𝜺𝜺 ⋅ 𝒙𝒙 𝑆𝑆 ℬ (a) (b) Fig. 8 (a) Setup of the boundary value problem for the RVE; (b) meshand a typical initial phase field. Manifold Learning Approach to Accelerate Phase Field Fracture Simulations in the Representative Volume Element 7
Here we made some simplifications for the micro cracksimulation so that we can better illustrate the main idea: (1)as mentioned in Section 2, the load is a unilateral tensionwith given displacement as shown in (3h), where the macro-scopic strain is ε = ε e ⊗ e ; (2) we only consider cracksin the matrix and ignore those on the interface and in thefiber; (3) the initial crack consists of two edges and threeconnected nodes, but nodes belonging to the same elementare forbidden.With the phase field values d = u ≡ , we minimize (1) to get an “equilibrated”phase field as a typical input X i . The totality of such inputsis termed the training set . The process of construction of thedata manifold with the training set is illustrated in Fig. 9. FEM SolutionTraining Set 𝑿𝑿 … 𝑿𝑿 𝑖𝑖 … 𝑿𝑿 𝑁𝑁 LLE interpolation LLEPrediction 𝒁𝒁 … 𝒁𝒁 𝑖𝑖∗ … 𝒁𝒁 𝑁𝑁∗
Error FEMOutput 𝒁𝒁 … 𝒁𝒁 𝑖𝑖 … 𝒁𝒁 𝑁𝑁 Fig. 9
The process of manifold learning using LLE.
For each input X i in the training set, we generate thehigh-fidelity solution of the evolved phase field through afinite element program, and the result is denoted Z i . Noticethat only input data are used during the LLE construction,while the output data { Z i } are only used for reconstruction.The output data are not limited to be the phase field solutionat the given load, nor need it have the same dimension as theinput datapoints.Once we obtain the data manifold, we reconstruct theoutput, marked by Z ∗ i , for every new input X ∗ i not in thetraining set through the following process:1. We find k ( (cid:28) N , which can be the same as k , see Sec-tion 4 for more details) nearest neighbors of X ∗ i in X and the corresponding weights in the high dimensionalspace R D , then we map X ∗ i to the low dimensional man-ifold Y ∗ i ∈ R L .2. We find the k nearest neighbors of Y ∗ i in Y , called S ∗ i ,and their weights W i j in the low dimensional manifold.Note that these neighbors may not correspond to thosein the previous step. 3. Locally linear reconstruct the output with weights andits k nearest neighbors in high dimensional data space: Z ∗ i = ∑ j ∈ S ∗ i W i j Z j . In this section, we detail the numerical implementation andresults of the proposed manifold learning method. In addi-tion, we provide a validation check for the computationalstrategy.4.1 Data GenerationIn our high-fidelity finite element analysis, the material con-stants are chosen as according to Table 1. The RVE size L = ε = ε e ⊗ e where ε = . × − . The regularized length scale parameter l is chosen such that h ≤ l /
2, where h is the mesh size. Werandomly generated 496 initial phase fields as detailed inSection 3.2, which correspond to 496 data points for train-ing (manifold learning). Table 1
Material parameters used in the high-fidelity finite elementsimulations λ (GPa) µ (GPa) λ (GPa) µ (GPa) g c ( mJ / mm ) l (mm)121.15 80.77 105.58 172.27 2.7 40 k and L , while the reconstruction process is defined by onehyperparameter k . Hence, the complete manifold model forthe problem requires three hyperparameters ( k , k , L ).The adopted parameter selection method is called crossvalidation (CV). Through CV we will select the best combi-nation of hyperparameters which leads to a balance of costand accuracy. The CV process is proceeded as follows. Wesplit the whole dataset ( N =
490 datapoints) to be n = X ( ) ,. . . , X ( n ) ,then we choose n − j th subset X ( j ) . Let Z ( j ) = { Z ( j ) i } denote the corre-sponding output phase field data for the validation set, and Liu, Weng, and Shen Z ∗ ( j ) = { Z ∗ ( j ) i } the LLE reconstruction. Then the final CVerror R reads R = n n ∑ j = ∑ i (cid:107) Z ∗ ( j ) i − Z ( j ) i (cid:107) l (cid:107) Z ( j ) i (cid:107) l . This procedure is illustrated in Fig. 10.
Validation Set Validation Set Validation Set Validation Set … Round 2Round 1 Round nRound 3
Validation Error 0.6 0.2 0.5 0.3Final CV Error = Average(Validation Error)
Fig. 10
The process of cross validation.
The procedure to select hyperparameters consists of twostages: (1) the dimension reduction process involving k ,and (2) the reconstruction process involving k . Iteratingthrough combinations of ( L , k , k ) with a fixed k value,an error matrix is deduced with columns denoting values of k / k , and rows denoting values of L as shown in Table 2.We find that k = k will yield a low CV error, which is rea-sonable, as the case k > k will lead to information loss inthe reconstruction process, and k < k will add noise to thereconstruction process. Table 2
CV error with different combinations of k / k and L , with k = k / k \ L
20 40 60 80 100 1201/4 0.4798 0.4042 0.3769 0.3599 0.3529 0.34421/2 0.4021 0.3645 0.3467 0.3338 0.3261 0.32221 0.3754 0.3472 0.3260 0.3114 0.3059 0.30362 0.3768 0.3472 0.3254 0.3103 0.3024 0.30074 0.3774 0.3457 0.3241 0.3096 0.3026 0.3001
Then we fix k = k and perform more CV to obtain Ta-ble 3, from which we determine that k = k =
20 gives arelatively low CV error for each L .Then, we plot the CV error as a function of L in Fig. 11.This figure indicates that an increase in L will reduce theaverage error, as expected. However, using a larger L in-creases the training time. Therefore, we follow the standardway to make the trade-off, i.e., to get the critical turningpoint at approximately the elbow, where L =
80. When L Table 3
CV error with different combinations of k (= k ) and L . k \ L
20 40 60 80 100 1205 0.6293 0.4905 0.4452 0.4229 0.4015 0.385010 0.4109 0.3665 0.3461 0.3343 0.3275 0.328115 0.3715 0.3432 0.3244 0.3147 0.3082 0.306120 0.3754 0.3472 0.3260 0.3114 0.3059 0.303625 0.3914 0.3582 0.3323 0.3148 0.3070 0.303630 0.4073 0.3694 0.3383 0.3175 0.3081 0.3046 is beyond this value, the error decreases at a very slow rate,while the training efficiency continually decreases. C V e rr o r Fig. 11
CV error vs. L . In conclusion, the chosen hyperparameters are ( k , k , L ) =( , , ) .4.3 Results and DiscussionTo remove the data bias, we generate a new set of 496 datapoints, which shares no data points with the set used in theparameter selection process. With the selected hyperparame-ters ( k , k , L ) = ( , , ) , we build the model using 464data points for training, and use the remaining 32 data pointsfor testing. To visualize the manifold built by the trainingdata, and together showing the test data, we perform an LLEreduction again for the 80-dimensional manifold to 2 dimen-sions, as in Fig. 12. It can be observed that the test datapoints are not far from the manifold trained from the trainingdata.Next we extract and visualize the nearest neighbors of acertain data point, as shown in Fig. 13(a) and (b).In Fig. 13(a), we observe that the nearest neighbors inthe training set are close to the chosen test data point (PointNo. 11). In Fig. 13(b), however, the nearest neighbors of Manifold Learning Approach to Accelerate Phase Field Fracture Simulations in the Representative Volume Element 9
Fig. 12
2D visualization of the 80D manifold(a)(b)
Fig. 13 (a) Nearest 20 neighbor points of test point No. 11; (b) Nearest20 neighbor points of test point No. 13. the chosen test data point (Point No. 13) appear scatteringaround. This phenomenon is still acceptable since the dis-tances between points in the remaining 78 dimensions arenot seen in the figures.We next visualize the cracked microstructures in Fig. 14,where we can observe a pattern that similar microstructurewill cluster in a continuous mode, showing the dimensionreduction is reasonable.4.4 Reconstruction Error Analysis for the Phase FieldIn this subsection, the output is the evolved phase field Z ∗ i = (cid:8) d j (cid:9) i , where j = , , ..., D . Therefore, the output Z ∗ i andinput X i have the same dimension. A histogram showingthe reconstruction errors is given in Fig. 15, where we usethe normalized l -norm to represent the error magnitude inthe output phase field, i.e., (cid:107) Z ∗ i − Z i (cid:107) l (cid:107) Z i (cid:107) l . (7)From this figure it can be seen that the LLE reconstructionerror for the phase field is acceptable.To examine the deciding factor of such error, we plot thenormalized error in l -norm of the 32 test points versus theirdistance to the manifold in Fig. 16. Here the distance of X ∗ i to the manifold is given by (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) X ∗ i − ∑ j ∈ S ∗ i W i j X j (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) l . A positive relationship between the reconstruction error andthis distance is observed, without outliers. Thus we can safelysay that if a test data point is close enough to the manifold,the reconstruction error of its microcrack propagation resultwill be small, guaranteeing the validity of this LLE manifoldlearning method.4.5 Reconstruction Error Analysis for the HomogenizedStressIn this subsection, the output is the homogenized stress Z ∗ i = σ i , where σ i = (cid:8) σ x , σ y , σ z , σ xy (cid:9) i , where for the plane straincase, σ z = ν ( σ x + σ y ) for the matrix and likewise for thefiber. As Fig. 1 shows, the homogenized stress is obtainedfrom the RVE through the volume average, σ = | B | (cid:90) ∂ B σ d B . Then the normalized reconstruction error in l -norm (7) be-comes (cid:107) σ ∗ i − σ i (cid:107) l (cid:107) σ i (cid:107) l . (8) Fig. 14
Crack microstructures mapped into the manifold described by the first two coordinates of LLE. Representative microstructures are shownnext to the square points. The solid line shows a continuous mode change of microstructures. l -norm of 32 test data points0246810 C o un t Fig. 15
Normalized l reconstruction error of the evolved phase field,i.e., (cid:107) Z ∗ i − Z i (cid:107) l / (cid:107) Z i (cid:107) l , of the 32 test data points. The normalized reconstruction error of the homogenized stressis shown in Fig. 17. It shows that the normalized reconstruc-tion error is smaller than 0 .
05, which is very small. Fig. 18shows that the reconstruction error is bounded by a factortimes the distance to the manifold, indicating a similar con-clusion, i.e., an a priori error estimate can be obtained. N o r m a li z e d e rr o r i n l n o r m regressionLLE prediction Fig. 16
Normalized l error of the evolved phase field vs. the distanceto the manifold. The l errors have a positive correlation with the dis-tance to the manifold. Remark.
Through the correlation of the reconstruction er-ror and the distance to the manifold, we can pre-determinewhether a new input data point X ∗ i is suitable for the man-ifold learning approach: if X ∗ i is close enough to the man-ifold, the reconstruction of the phase field at the given loadwill be accurate; otherwise, if it is far away from the man-ifold, we should either not use the manifold reconstruction Manifold Learning Approach to Accelerate Phase Field Fracture Simulations in the Representative Volume Element 11 l -norm of 32 test data points012345678 C o un t Fig. 17
Normalized l reconstruction error of the homogenized stress,i.e., (cid:107) σ ∗ i − σ i (cid:107) l / (cid:107) σ i (cid:107) l , of the 32 test data points. N o r m a li z e d e rr o r i n l n o r m Fig. 18
Normalized l error of the homogenized stress vs. the distanceto the manifold. The l errors are bounded by a factor times the distanceto the manifold. for this particular input, or augment the training set with X ∗ i .This property can also be exploited to aid an adaptivity pro-cedure to augment the training set on the fly: if the distancefrom a certain new input X ∗ i to its manifold projection is toohigh, then we can add it (and its output from high-fidelitycomputation) to the training set. We have proposed a manifold learning approach to acceler-ate phase field fracture simulations in the RVE in the contextof the FE scheme. Considering a group of RVEs with thesame microstructure except for the microcracks, we use thephase field approach to represent such microcracks. We thenmake use of the LLE technique to construct a data manifoldthat contains a collection of similar cracked microstructures(RVEs). This LLE manifold can be used to efficiently andaccurately predict the phase field output as a function of the initial phase field, provided that all the analysis is done atthe same load applied to the RVE. The same approach canbe generalized to cases with more complicated RVEs suchas elastoplastic constitutive behavior.This new computational approach enjoys the followingfeatures:1. Only three hyperparameters need to be determined tolearn the manifold. And once the data manifold is con-structed, minimum computation is required to reconstructthe phase field output.2. There exists an indicator which can pre-estimate the re-construction error and pre-determine whether an inputdata is suitable to perform the reconstruction. We wouldlike to emphasize that this feature is very desirable, sincecompared with more popular machine-learning techniquessuch as neural networks – in many of those techniques, itis difficult to predict whether an interpolation is accurateor not without knowing the exact solution.3. A number of generalizations can be made, e.g., to threedimensions, and to the types of RVEs, boundary condi-tions, and outputs. In fact, the output can be of a highdimension, as long as there exists a continuous depen-dence of the output on the input, which is anyway a pre-requisite of a well-posed problem.4. The applicability of this approach is promising. The adap-tive algorithm makes efficient multiscale fracture simu-lation possible. Conflict of interest
The authors declare that they have no conflict of interest.
References
1. F. Feyel, Computational Materials Science (1), 344 (1999)2. V.B.C. Tan, K. Raju, H.P. Lee, Computer Methods in Applied Me-chanics and Engineering , 112694 (2020)3. B. Bourdin, G.A. Francfort, J.J. Marigo, Journal of the Mechanicsand Physics of Solids (4), 797 (2000)4. A.A. Griffith, Philosophical Transactions of the Royal Society ofLondon. Series A, Containing Papers of a Mathematical or Physi-cal Character (582–593), 163 (1921)5. W. Lai, J. Gao, Y. Li, M. Arroyo, Y. Shen, Computer Methods inApplied Mechanics and Engineering , 112787 (2020)6. F. Amiri, D. Milln, Y. Shen, T. Rabczuk, M. Arroyo, Theoreticaland Applied Fracture Mechanics , 102 (2014)7. P. Zhang, X. Hu, T.Q. Bui, W. Yao, International Journal of Me-chanical Sciences , 105008 (2019)8. P. Zhang, Y. Feng, T.Q. Bui, X. Hu, W. Yao, Composite Structures , 111551 (2020)9. T.T. Nguyen, D. Waldmann, T.Q. Bui, Computer Methods in Ap-plied Mechanics and Engineering , 1 (2019)10. T.T. Nguyen, D. Waldmann, T.Q. Bui, Journal of ComputationalPhysics , 585 (2019)11. M. Mollaali, V. Ziaei-Rad, Y. Shen, Journal of Natural Gas Sci-ence and Engineering , 102905 (2019)2 Liu, Weng, and Shen12. T. Heister, M.F. Wheeler, T. Wick, Computer Methods in AppliedMechanics and Engineering , 466 (2015)13. Y. Li, W. Lai, Y. Shen, International Journal of Fracture (1-2),83 (2019)14. V. Ziaei-Rad, Y. Shen, Computer Methods in Applied Mechanicsand Engineering , 224 (2016)15. T. Gerasimov, L.D. Lorenzis, Computer Methods in Applied Me-chanics and Engineering , 276 (2016)16. T. Wick, SIAM Journal on Scientific Computing (4), B589(2017)17. T. Wick, Computer Methods in Applied Mechanics and Engineer-ing , 577 (2017)18. P. Farrell, C. Maurini, International Journal for Numerical Meth-ods in Engineering (5), 648 (2017)19. J.Y. Wu, Y. Huang, V.P. Nguyen, Computer Methods in AppliedMechanics and Engineering , 112704 (2020)20. A. Kopaniˇc´akov´a, R. Krause, Computer Methods in Applied Me-chanics and Engineering , 112720 (2020)21. S. Bhattacharjee, K. Matou, Journal of Computational Physics , 635 (2016)22. D. Wirtz, N. Karajan, B. Haasdonk, International Journal for Nu-merical Methods in Engineering (1), 1 (2015)23. C. Wang, S. Mahadevan, in Twenty-Seventh AAAI Conference onArtificial Intelligence (2013)24. J. Yvonnet, Q.C. He, Journal of Computational Physics (1),341 (2007)25. K. Matou, M.G. Geers, V.G. Kouznetsova, A. Gillman, Journal ofComputational Physics , 192 (2017)26. S.T. Roweis, L.K. Saul, Science (5500), 2323 (2000)27. J. Chen, Z. Ma, International Journal of Pattern Recognition andArtificial Intelligence (07), 985 (2011)28. E. Lopez, D. Gonzalez, J.V. Aguado, E. Abisset-Chavanne,E. Cueto, C. Binetruy, F. Chinesta, Archives of ComputationalMethods in Engineering (1), 59 (2018)29. J. Guilleminot, J.E. Dolbow, Mechanics Research Communica-tions p. 103443 (2020)30. R. Iba˜nez, E. Abisset-Chavanne, J.V. Aguado, D. Gonzalez,E. Cueto, F. Chinesta, Archives of Computational Methods in En-gineering (1), 47 (2018)31. T. Furukawa, G. Yagawa, International Journal for NumericalMethods in Engineering (2), 195 (1998)32. J. Ghaboussi, D.A. Pecknold, M. Zhang, R.M. Haj-Ali, Interna-tional Journal for Numerical Methods in Engineering (1), 105(1998)33. Y.M.A. Hashash, S. Jung, J. Ghaboussi, International Journal forNumerical Methods in Engineering (7), 989 (2004)34. S. Jung, J. Ghaboussi, Computers & Structures (15), 955 (2006)35. Y. Sun, W.D. Zeng, Y.Q. Zhao, Y.L. Qi, X. Ma, Y.F. Han, Compu-tational Materials Science (3), 686 (2010)36. G. Ji, F. Li, Q. Li, H. Li, Z. Li, Materials Science and Engineering:A (13), 4774 (2011)37. R. Liu, T. Ruan, S. Song, Y. Lin, G. Jiang, Journal of Chromatog-raphy A , 13 (2015)38. R. Kondo, S. Yamakawa, Y. Masuoka, S. Tajima, R. Asahi, ActaMaterialia , 29 (2017)39. R. Cang, H. Li, H. Yao, Y. Jiao, Y. Ren, Computational MaterialsScience , 212 (2018)40. H. Wei, S. Zhao, Q. Rong, H. Bao, International Journal of Heatand Mass Transfer , 908 (2018)41. X. Li, Z. Liu, S. Cui, C. Luo, C. Li, Z. Zhuang, Computer Methodsin Applied Mechanics and Engineering , 735 (2019)42. H. Amor, J.J. Marigo, C. Maurini, Journal of the Mechanics andPhysics of Solids (8), 1209 (2009)43. C. Miehe, F. Welschinger, M. Hofacker, International Journal forNumerical Methods in Engineering (10), 1273 (2010)44. G.A. Francfort, J.J. Marigo, Journal of the Mechanics and Physicsof Solids (8), 1319 (1998) 45. Y. Shen, M. Mollaali, Y. Li, W. Ma, J. Jiang, Journal of ShanghaiJiao Tong University (Science) (1), 166 (2018)46. B. Sch¨olkopf, A. Smola, K.R. M¨uller, Neural Computation10