SYNTHESIZED ENRICHMENT FUNCTIONS FOR EXTENDED FINITE ELEMENT ANALYSES WITH FULLY RESOLVED MICROSTRUCTURE

Inspired by the first order numerical homogenization, we present a method for extracting continuous fluctuation fields from the Wang tile based compression of a material microstructure. The fluctuation fields are then used as enrichment basis in Extended Finite Element Method (XFEM) to reduce number of unknowns in problems with fully resolved microstructural geometry synthesized by means of the tiling concept. In addition, the XFEM basis functions are taken as reduced modes of a detailed discretization in order to circumvent the need for non-standard numerical quadratures. The methodology is illustrated with a scalar steady-state problem.


Introduction
All materials exhibit heterogeneous microstructure at certain level of detail.Besides the physical properties of microstructural constituents, the overall response is also dictated by their spatial arrangement.With the advance of numerical methods and computational power, there is an emphasis on devising approaches that take the microstructural characteristics into account and propagate knowledge of material composition into upper-scale models.
In our contribution, we focus on problems in which the underlying microstructural geometry is explicitly known.We consider problems where the characteristic length of the microstructure is comparable to the dimension of an upper-scale model; hence, the nested numerical homogenization, such as the FE 2 method [1], is not applicable because the underlying assumption of the separation of scales is violated [2, and references therein].
This work supplements our previous results regarding efficient modelling of materials with stochastic heterogeneous microstructure, e.g., [3,4].The microstructural geometry of a macro-scale model-as well as its finite element discretization-is synthesized using the Wang tiling concept, recalled in Section 2.
Here, we present a method to extract dominant responses of a compressed microstructure to prescribed macroscopic loading.We reuse the extracted responses to synthesize a microstructure-informed enrichment basis in order to accelerate the upper scale calculations.The numerical scheme, outlined in Section 4, combines eXtended Finite Element Framework (XFEM) [5, 6, and references therein], used for definition of local basis functions, with the reduced order modelling approach [7], allowing for utilization of standard finite element procedures.Finally, Section 5 illustrates the methodology with a two-dimensional problem of steady-state thermal conduction in an L-shaped domain.

Wang tiling concept
Microstructure of random heterogeneous materials is usually represented by means of a periodic finitesize model-appearing in the literature under various names such as Statistically Optimal Representative Unit Cell, Statistically Similar Representative Volume Element, or Statistically Equivalent Periodic Unit Cell (SEPUC)-that resembles the material in terms of various spatial statistics.This setting is particularly suitable for the family of FE 2 simulations,e.g., [1], in which it is often accompanied with periodic boundary conditions imposed at the micro-scale level.However, when used for microstructure synthesis, this representation introduces unrealistic, spurious periodic correlations.
The formalism of Wang tiles represents a compromise between the amount of induced correlation and complexity of a microstructure representation.Instead of a single cell, microstructural information is compressed into a set of equi-sized Wang tiles with pre-defined requirements on mutual compatibility.In the abstract setting, the compatibility constraints are represented by codes assigned to tile edges, illustrated with distinct colours in Figure 1.A realization of the compressed microstructure is then synthesized following a simple stochastic algorithm: an empty grid of a desired size is sequentially filled with instances of Wang tiles such that the adjacent tile edges are compatible.At each position, a tile is randomly chosen from a subset of tiles compliant with previously placed tiles; hence, the algorithm generates stochastic realizations, providing that set of tiles is rich enough such that multiple tiles are admissible at each assembly step; see Figure 1 for an illustration.Procedures designed to compress a given microstructure into a SEPUC can be straightforwardly extended to take into account generalized periodicity occurring in the tiling concept, e.g., the standard optimization procedure based on minimizing discrepancy in the two-point probability function was used in [3].We also adapted a sample-based approach originated in Computer Graphics in order to address the high computational cost of the optimization-based design [4].Currently, a variety of material microstructures, ranging from particulate media to complex foam-like microstructures, can be represented with the framework of Wang tiling.Individual realizations of the microstructure are then generated in linear time with respect to the required sample size, which makes the tiling concept appealing for applications where multiple microstructure samples are required, e.g., investigation of the RVE size [8].

Synthesized enrichments fields
The idea of compressing fluctuation fields using the Wang tiling concept has been introduced in Novák et al. [9], where the microstructural geometry of individual tiles was optimized with respect to (i) meeting prescribed spatial statistics and (ii) minimizing traction discrepancies among congruent edges in an auxiliary tiling.However, the resulting synthesized stress fields exhibited unavoidable discontinuities and the objective function required careful tuning of weight parameters because the two requirements turned to be hard to achieve simultaneously.In contrast, the present method generates continuous fields of the primal unknowns for any existing microstructure compression.
We illustrate the methodology with the steady-state heat conduction problem, governed by where θ denotes the unknown temperature field, K is a local heat conductivity tensor and no heat sources are considered.First, we assume the decomposition of θ into a fluctuation part θ, caused by the presence of heterogeneities, and a macroscopic part controlled with a macroscopic temperature gradient H, where T (i) denotes the domain of the i-th tile and S = T (i) , i = 1 . . .n T represent the set of n T available tiles.Without loss of generality, each tile is assumed centred, i.e., 2 , where 2α is the size of the tile edge.
In order to enforce continuity of the extracted fields, all tiles are loaded with the prescribed macroscopic gradient H and solved simultaneously.Providing that the tile geometries are discretized with a compatible triangulation-meaning that discretization of the corresponding tile edges is the same across all tiles-the continuity is ensured by associating the corresponding edge unknowns during the localization of an element contribution into the conductivity matrix of the whole set.The resulting system takes the standard form where is the conductivity matrix and F (H) is a H-dependent loading vector arising from the weak form after plugging Eq. (2) into Eq.( 1).K S has block-diagonal structure similar to primal domain decomposition approaches.
In addition, two types of constraints are imposed on the set system.First, we eliminate zero-energy modes of K S by requirement The second constraint prevents the fluctuation field θ from compensating for the loading induced through H.
Inspired in numerical homogenization, e.g.[2], we introduce the second constraint enforcing in a chosen domain Ω loaded analogously to Eq. ( 2).Using Green's first identity, Eq. ( 5) can be recast into a boundary integral where n(x) denotes the outer normal of the boundary.
In particular, we consider three mutually exclusive types of the second constraint: (K) The first type prescribes zero temperature fluctuations at tile boundaries, θ| where ∂T (i) stands for the boundary of the i-th tile.
Note that (K) corresponds directly to Kinematic Uniform Boundary Conditions used in numerical homogenization, while (P) and (S) mimic Periodic and Statically Uniform Boundary Conditions, respectively.Instead of Eq. ( 3), the final fluctuations θ follows from   where C I and C II matrices represent the discretized versions of Eq. ( 4) and one of Eqs. ( 7), (8), or (9), respectively, and λ • stands for the corresponding vector of Lagrange multipliers.While C II straightforwardly eliminates edge unknowns for the constraint type (K), constructing contributions of individual tile separately and localizing them into C II leads to an ill-posed problem for types (P) and (S) From Eq. ( 6), it clearly follows that if two opposite sides of a tile are required to be compatiblehence they feature the same trace of θ-the corresponding component of ∇ θ average vanishes and so does the tile contribution to the relevant row in C II .In addition, the ill-posedness also stems from the repetition of edge code pairs within the tile set (swapping the order of codes changes only the sign of the edge contribution to Eq. ( 5)).In certain cases, this can even result in zero matrix C II for the least restrictive constraint (S), in which tile contributions are summed together.
One way to treat the singularity of C II is to perform Singular Value Decomposition on the constraint and take into account only the non-zero singular values; however, such a procedure is sensitive to numerical precision.Thanks to the simple structure of the constraint, the regular part of the constraint can be directly established as a projection of the original matrix C II .First, all n c combinations of edge codes pertinent to either horizontal or vertical edges are identified and enumerated.Next, an empty matrix P ∈ R 2n T ×nc , for (P) type of the constraint, or P ∈ R 2×nc , for (S) type, is created and populated with 1 (or −1) if the tile contains an enumerated pair of codes (or a pair in the inverse order).Odd and even rows of the matrix correspond to vertical and horizontal edges, respectively.If P remains empty for (S), C II constraint is discarded entirely.Otherwise, the non-singular constraint matrix C II , replacing C II in Eq. ( 10), is obtained from a projection C II = P T C II .
The outlined strategy allows to define one enrichment field θ for a prescribed H and the selected type of the second constraint.For the linear problem considered in this work, two mutually orthonormal loading cases H = 1, 0 T and H = 0, 1 T cover all possible load cases.Combined with three types of C II , we can generate up to six enrichment fields.Note that unlike the nested numerical homogenization, where the choice of boundary conditions is based on their suitability in capturing effects of a surrounding medium, we combine different types of C II to control cardinality of the set of compressed fluctuations.

Numerical strategy
The fields introduced above are employed as enrichments in a upper-scale model whose microstructural geometry is assembled by means of tiling from the set depicted in Figure 1.Assume a upper-scale problem governed by the same differential equation as in Eq. ( 1), valid in a upperscale domain Ω and accompanied by boundary conditions where θ and q are the given temperature and normal heat flux profiles, respectively, prescribed at parts Γ θ and Γ q of the domain boundary ∂Ω.An approximate numerical solution arises from the weak form of the problem and substitutes for the analytical solution that is usually intractable.In Eq. ( 13), a(θ, ϑ) and b(ϑ) denote the bilinear and linear form pertinent to Eq. ( 1); V and V 0 stand for the space of admissible temperature fields.Quality of the approximate solution is governed by suitability of a finite-dimensional subspace V c ⊂ V in which the solution is sought for.In particular, domain discretization must accurately resolve all geometrical details in the standard Finite Element (FE) setting.This requirement leads to very fine meshes in applications where microstructural geometry of the domain is involved.
Besides local mesh refinements, the approximation space can be enhanced by incorporating prior knowledge of local character of the solution.Concretely, in the Extended Finite Element framework [5,6], the approximate solution θ c ∈ V c takes the form where the first sum contains the standard finite element shape functions N i (x) and the second term adds n e enrichment functions ψ j for each discretization node; θ i denotes the regular Degree Of Freedom (DOF) associated with finite-element mesh node x i , θ j i are the additional DOFs.
In traditional XFEM, only patches of the domain are usually enriched with ψ j , suited for one geometrical feature.Typically, the enrichment provides asymptotic solution near a crack tip or introduces strong/weak discontinuities [6].Here, we add global enrichments that capture collective response of a material microstructure.The original shape functions N i (x) model the macroscopic response of the domain, while products N i (x)ψ j (x) cover the fluctuations caused by material heterogeneity.A similar idea was introduced by Strouboulis et al. [10], who used numerical "handbook" functions for assemblies of closely packed inclusions, or Plews and Duarte [11].However, both approaches solve local boundary value problems defined on subdomains first and subsequently run the analysis of the whole domain.In contrast, we precompute the fluctuations purely "off-line" at the level of microstructure compression.i.e., without any information of the domain geometry or loading, and construct the enrichments as an assembly of the tiledefined fluctuation fields.
In our approach, the upper scale discretization is constructed irrespectively of the underlying microstructural geometry of Ω.Hence, the approximation basis functions in the form of Eq. ( 14) together with fluctuations in material parameters due to the presence of microstructural geometry preclude efficient use of standard numerical procedures for evaluation of the bi-linear form in Eq. ( 13).
In order to circumvent this drawback, we construct another, finer domain discretization space V f , assembled from the tile discretization used in the off-line phase for extracting fluctuation fields.The shape functions N i (x), defined at the coarser discretization, are projected onto the fine mesh; enrichments ψ j (x) are already defined within the fine discretization by construction.Consequently, Eq. ( 14) can be understood as a definition of reduced modes for the fine discretization.
Let matrix Φ comprises individual reduced modes as its columns, computed as an element-wise product of the projected N i (x) and ψ j (x) following Eq.( 14).The range of Φ thus defines a subspace of V f corresponding to V c .Instead of solving the large, fine-discretization system K f θ f = F f , we restrict the fine-discretization unknowns θ f via where a denotes the vector of unknown coefficients pertinent to the reduced modes.The final algebraic system then takes the form

Results
We considered the material microstructure of monodisperse elliptical inclusions, shown in Figure 1.Both phases were assumed isotropic, with the inclusion material being more conductive (K(x) = 100 I) than the matrix phase (K(x) = 1 I).Six fluctuation fields in total, discretized with linear triangular elements, were extracted, see Figure 2. The microstructure was superimposed to an upperscale problem represented by an L-shaped domain with uniform temperature profiles prescribed at the bottom and right-hand side edges.The initial coarse discretization of the domain is depicted in Figure 3.The linear Lagrange basis functions were assumed at the coarse level, whilst the fine space comprised the quadratic Lagrange approximation arising from the product of linear coarse shape functions and the enrichments.
The solution θ XFEM of Eq. ( 16) was compared against Direct Numerical Simulation (DNS) θ DNS , which was taken as the reference solution.A sequence of five uniformly refined coarse discretizations, the first one shown in Figure 3

Conclusions
We have presented a method for extracting continuous fluctuation fields from a microstructure compressed by means of Wang tiles.We have also demonstrated utilization of these fields as enrichments in XFEM.
In order to avoid non-standard quadratures, XFEM basis functions were projected onto a discretization arising from assembly of the tile discretization used for computing the fluctuation fields.
The proposed methodology was applied to a steadystate heat conduction problem defined in an L-shaped domain with microstructural geometry provided by the tiling concept.Only 1.6 % of the original degrees of freedom were sufficient to obtain 2 % relative error compared to DNS.
On the other hand, the reduction in DOFs accelerates only the solution of Eq. ( 16).For more significant time savings, a low-rank additional approximation of K f , e.g. using the idea of hyper-reduction [12], has to be introduced and is the focus of our current work.We have also restricted ourself to the first-order decomposition, Eq. ( 2); higher order expansion of the macroscopic loading can be considered to further enrich the approximation space.

Figure 1 .
Figure 1.The compressed microstructural geometry represented with the set of Wang tiles (top) and an illustration of a step in the tiling algorithm with partially assembled microstructure (bottom).

Figure 2 .
Figure 2. A compressed fluctuation field obtained by prescribing (S) variant of the second constraint and setting H = {1, 0} T .

Figure 3 .
Figure 3.A scheme of the problem definition.The edge lengths tx and ty of the tiles were set as tx = ty = 6.All quantities are in consistent units.

u 2 H1 = Ω u 2
(x) + ∇u(x) • ∇u(x) dx .(20)Mean and standard deviations of the errors computed for five coarse discretizations and six different tiling realizations are given in Figure4.

Figure 4 .
Figure 4. Convergence of L2 and H1 errors with increasing number of DOFs in the coarse discretization.