Acta Polytechnica CTU Proceedings doi:10.14311/APP.2017.13.0029 Acta Polytechnica CTU Proceedings 13:29–34, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/app SYNTHESIZED ENRICHMENT FUNCTIONS FOR EXTENDED FINITE ELEMENT ANALYSES WITH FULLY RESOLVED MICROSTRUCTURE Martin Doškářa,∗, Jan Nováka, Jan Zemana,b a Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7/2077, 166 29 Prague 6, Czech Republic b The Institute of Information Theory and Automation, Academy of Sciences of the Czech Republic, Pod Vodárenskou věží 4, 182 08 Prague 8, Czech Republic ∗ corresponding author: martin.doskar@fsv.cvut.cz Abstract. 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. Keywords: Wang tiling, microstructure synthesis, microstructure-informed enrichment functions, extended finite element method. 1. 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 ac- count and propagate knowledge of material composi- tion 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 character- istic length of the microstructure is comparable to the dimension of an upper-scale model; hence, the nested numerical homogenization, such as the FE2 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 regard- ing efficient modelling of materials with stochastic heterogeneous microstructure, e.g., [3, 4]. The mi- crostructural 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 re- sponses 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 fi- nite element procedures. Finally, Section 5 illustrates the methodology with a two-dimensional problem of steady-state thermal conduction in an L-shaped do- main. 2. Wang tiling concept Microstructure of random heterogeneous materials is usually represented by means of a periodic finite- size 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 FE2 simulations,e.g., [1], in which it is often accompanied with periodic boundary conditions imposed at the micro-scale level. How- ever, when used for microstructure synthesis, this rep- resentation introduces unrealistic, spurious periodic correlations. The formalism of Wang tiles represents a compro- mise between the amount of induced correlation and complexity of a microstructure representation. In- stead 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 follow- ing a simple stochastic algorithm: an empty grid of 29 http://dx.doi.org/10.14311/APP.2017.13.0029 http://ojs.cvut.cz/ojs/index.php/app M. Doškář, J. Novák, J. Zeman Acta Polytechnica CTU Proceedings a desired size is sequentially filled with instances of Wang tiles such that the adjacent tile edges are com- patible. At each position, a tile is randomly chosen from a subset of tiles compliant with previously placed tiles; hence, the algorithm generates stochastic real- izations, providing that set of tiles is rich enough such that multiple tiles are admissible at each assembly step; see Figure 1 for an illustration. Figure 1. The compressed microstructural geome- try represented with the set of Wang tiles (top) and an illustration of a step in the tiling algorithm with partially assembled microstructure (bottom). Procedures designed to compress a given microstruc- ture 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 com- putational cost of the optimization-based design [4]. Currently, a variety of material microstructures, rang- ing from particulate media to complex foam-like mi- crostructures, can be represented with the framework of Wang tiling. Individual realizations of the mi- crostructure are then generated in linear time with respect to the required sample size, which makes the tiling concept appealing for applications where multi- ple microstructure samples are required, e.g., investi- gation of the RVE size [8]. 3. 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 individ- ual tiles was optimized with respect to (i) meeting prescribed spatial statistics and (ii) minimizing trac- tion discrepancies among congruent edges in an auxil- iary 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 com- pression. We illustrate the methodology with the steady-state heat conduction problem, governed by ∇ · (−K(x) · ∇θ(x)) = 0 , x ∈ Ω , (1) 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 fluc- tuation part θ̃, caused by the presence of hetero- geneities, and a macroscopic part controlled with a macroscopic temperature gradient H, θ(x) = θ̃(x) + H · x , ∀x ∈ T(i),∀T(i) ∈ S, (2) where T(i) denotes the domain of the i-th tile and S ={ T(i), i = 1 . . . nT } represent the set of nT available tiles. Without loss of generality, each tile is assumed centred, i.e., T(i) = [−α, α]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 cor- responding 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 KS θ̃ = F (H) , (3) where KS = AT(i)∈S KT(i) 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). KS 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 KS by requirement 1 |S| ∫ S θ̃(x) dx = 0 . (4) The second constraint prevents the fluctuation field θ̃ from compensating for the loading induced through H . 30 vol. 13/2017 Synthesized microstructure-informed enrichments for XFEM analyses Inspired in numerical homogenization, e.g. [2], we introduce the second constraint enforcing 1 |Ω| ∫ Ω ∇θ̃(x) dx = 0 (5) in a chosen domain Ω loaded analogously to Eq. (2). Using Green’s first identity, Eq. (5) can be recast into a boundary integral∫ ∂Ω θ̃(x)n(x) dS = 0 , (6) 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 fluctu- ations at tile boundaries, θ̃|∂T(i) = 0 , ∀i = 1 . . . nT , (7) where ∂T(i) stands for the boundary of the i-th tile. (P) In the second type, the requirement (6) is posed tile-wise, i.e.,∫ ∂T(i) θ̃(x)n(x) dS = 0 , ∀i = 1 . . . nT . (8) (S) The last type imposes Eq. (6) over the whole set,∫ ∂S θ̃(x)n(x) dS = 0 , (9) where ∂S = { ∂T(i), i = 1 . . . nT } . Note that (K) corresponds directly to Kinematic Uni- form Boundary Conditions used in numerical homog- enization, while (P) and (S) mimic Periodic and Stat- ically Uniform Boundary Conditions, respectively. Instead of Eq. (3), the final fluctuations θ̃ follows from KS CT I CT II CI 0 0 CII 0 0  θ̃ λI λII  = F(H) 0 0  , (10) where CI and CII matrices represent the discretized versions of Eq. (4) and one of Eqs. (7), (8), or (9), re- spectively, and λ• stands for the corresponding vector of Lagrange multipliers. While CII straightforwardly eliminates edge un- knowns for the constraint type (K), constructing con- tributions of individual tile separately and localizing them into CII 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 compatible— hence they feature the same trace of θ̃—the corre- sponding component of ∇θ̃ average vanishes and so does the tile contribution to the relevant row in CII. In addition, the ill-posedness also stems from the rep- etition 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 CII for the least restrictive constraint (S), in which tile contributions are summed together. One way to treat the singularity of CII 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 con- straint, the regular part of the constraint can be di- rectly established as a projection of the original matrix CII. First, all nc combinations of edge codes pertinent to either horizontal or vertical edges are identified and enumerated. Next, an empty matrix P ∈ R2nT ×nc , for (P) type of the constraint, or P ∈ R2×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 cor- respond to vertical and horizontal edges, respectively. If P remains empty for (S), CII constraint is discarded entirely. Otherwise, the non-singular constraint ma- trix CII, replacing CII in Eq. (10), is obtained from a projection CII = PTCII. The outlined strategy allows to define one enrich- ment field θ̃ for a prescribed H and the selected type of the second constraint. For the linear problem consid- ered 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 CII, we can generate up to six enrichment fields. Note that un- like the nested numerical homogenization, where the choice of boundary conditions is based on their suit- ability in capturing effects of a surrounding medium, we combine different types of CII to control cardinality of the set of compressed fluctuations. Figure 2. A compressed fluctuation field obtained by prescribing (S) variant of the second constraint and setting H = {1, 0}T. 4. Numerical strategy The fields introduced above are employed as enrich- ments 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 upper- scale domain Ω and accompanied by boundary condi- tions θ(x) = θ̂(x) x ∈ Γθ 6= ∅ , (11) (−K(x)∇θ(x)) · n(x) = q̂(x) x ∈ Γq , (12) 31 M. Doškář, J. Novák, J. Zeman Acta Polytechnica CTU Proceedings 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 Find θ ∈ V such that a(θ, ϑ) = b(ϑ), ∀ϑ ∈ V0 . (13) and substitutes for the analytical solution that is usu- ally intractable. In Eq. (13), a(θ, ϑ) and b(ϑ) denote the bilinear and linear form pertinent to Eq. (1); V and V0 stand for the space of admissible temperature fields. Quality of the approximate solution is governed by suitability of a finite-dimensional subspace Vc ⊂ 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 applica- tions where microstructural geometry of the domain is involved. Besides local mesh refinements, the approximation space can be enhanced by incorporating prior knowl- edge of local character of the solution. Concretely, in the Extended Finite Element framework [5, 6], the approximate solution θc ∈ Vc takes the form θc(x) = nc∑ i=1 Ni(x) θi + nc∑ i=1 ne∑ j=1 Ni(x)ψj(x) θji , (14) where the first sum contains the standard finite el- ement shape functions Ni(x) and the second term adds ne enrichment functions ψj for each discretiza- tion node; θi denotes the regular Degree Of Freedom (DOF) associated with finite-element mesh node xi, θji are the additional DOFs. In traditional XFEM, only patches of the domain are usually enriched with ψj , suited for one geomet- rical 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 ma- terial microstructure. The original shape functions Ni(x) model the macroscopic response of the do- main, while products Ni(x)ψj(x) cover the fluctua- tions caused by material heterogeneity. A similar idea was introduced by Strouboulis et al. [10], who used nu- merical “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 pre- compute the fluctuations purely “off-line” at the level of microstructure compression. i.e., without any in- formation of the domain geometry or loading, and construct the enrichments as an assembly of the tile- defined fluctuation fields. In our approach, the upper scale discretization is constructed irrespectively of the underlying mi- crostructural geometry of Ω. Hence, the approxima- tion 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 Vf , assem- bled from the tile discretization used in the off-line phase for extracting fluctuation fields. The shape functions Ni(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 under- stood 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 projectedNi(x) and ψj(x) following Eq. (14). The range of Φ thus defines a subspace of Vf corresponding to Vc. Instead of solving the large, fine-discretization system Kfθf = Ff , we restrict the fine-discretization unknowns θf via θf = Φ a , (15) where a denotes the vector of unknown coefficients pertinent to the reduced modes. The final algebraic system then takes the form ΦTKfΦa = ΦTFf . (16) 5. Results We considered the material microstructure of mono- disperse 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 upper- scale 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, was investigated. The proximity of the reference (DNS) and XFEM solutions 32 vol. 13/2017 Synthesized microstructure-informed enrichments for XFEM analyses θ̂ = 20 10 tx 10 ty θ̂ = 0 10 tx 10 ty 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. was quantified via relative errors εL2 = ∥∥θDNS − θXFEM ∥∥ L2 ‖θDNS‖L2 and (17) εH1 = ∥∥θDNS − θXFEM ∥∥ H1 ‖θDNS‖H1 , (18) where ‖u‖2L2 = ∫ Ω u2(x) dx and (19) ‖u‖2H1 = ∫ Ω u2(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 Figure 4. 10!4 10!3 10!2 XFEM/DNS DOFs 10!2 10!1 0 0 L2 H1 Figure 4. Convergence of L2 and H1 errors with increasing number of DOFs in the coarse discretiza- tion. 6. 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 steady- state 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 acceler- ates only the solution of Eq. (16). For more significant time savings, a low-rank additional approximation of Kf , 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 de- composition, Eq. (2); higher order expansion of the macroscopic loading can be considered to further en- rich the approximation space. Acknowledgements M. Doškář and J. Novák acknowledge the endowment of the Ministry of Industry and Trade of the Czech Republic under the project No. FV10202. In addi- tion, M. Doškář was supported by the Grant Agency of the Czech Technical University in Prague, project No. SGS17/042/OHK1/1T/11 “Numerical Methods for Mod- eling Uncertainties in Civil Engineering”. References [1] F. Feyel, J.-L. Chaboche. FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Computer Methods in Applied Mechanics and Engineering 183(3-4):309–330, 2000. doi:10.1016/S0045-7825(99)00224-8. [2] M. G. D. Geers, V. G. Kouznetsova, K. Matouš, J. Yvonnet. Homogenization Methods and Multiscale Modeling: Nonlinear Problems. In E. Stein, R. de Borst, T. J. R. Hughes (eds.), Encyclopedia of Computational Mechanics Second Edition, pp. 1–34. John Wiley & Sons, Ltd, Chichester, UK, 2017. DOI: 10.1002/9781119176817.ecm107. [3] J. Novák, A. Kučerová, J. Zeman. Compressing random microstructures via stochastic Wang tilings. Physical Review E 86(4):040104(R), 2012. doi:10.1103/PhysRevE.86.040104. [4] M. Doškář, J. Novák, J. Zeman. Aperiodic compression and reconstruction of real-world material systems based on Wang tiles. Physical Review E 90(6):062118, 2014. doi:10.1103/PhysRevE.90.062118. [5] T.-P. Fries, T. Belytschko. The extended/generalized finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering 84(3):253–304, 2010. doi:10.1002/nme.2914. [6] T. Strouboulis, K. Copps, I. Babuška. The generalized finite element method. Computer Methods in Applied Mechanics and Engineering 190(32-33):4081–4193, 2001. doi:10.1016/S0045-7825(01)00188-8. 33 http://dx.doi.org/10.1016/S0045-7825(99)00224-8 http://dx.doi.org/10.1103/PhysRevE.86.040104 http://dx.doi.org/10.1103/PhysRevE.90.062118 http://dx.doi.org/10.1002/nme.2914 http://dx.doi.org/10.1016/S0045-7825(01)00188-8 M. Doškář, J. Novák, J. Zeman Acta Polytechnica CTU Proceedings [7] J. Yvonnet, Q.-C. He. The reduced model multiscale method (R3m) for the non-linear homogenization of hyperelastic media at finite strains. Journal of Computational Physics 223(1):341–368, 2007. doi:10.1016/j.jcp.2006.09.019. [8] M. Doškář, J. Novák. A jigsaw puzzle framework for homogenization of high porosity foams. Computers & Structures 166:33–41, 2016. doi:10.1016/j.compstruc.2016.01.003. [9] J. Novák, A. Kučerová, J. Zeman. Microstructural enrichment functions based on stochastic Wang tilings. Modelling and Simulation in Materials Science and Engineering 21(2):025014, 2013. doi:10.1088/0965-0393/21/2/025014. [10] T. Strouboulis, L. Zhang, I. Babuška. Generalized finite element method using mesh-based handbooks: application to problems in domains with many voids. Computer Methods in Applied Mechanics and Engineering 192(28-30):3109–3161, 2003. doi:10.1016/S0045-7825(03)00347-5. [11] J. Plews, C. Duarte. Bridging multiple structural scales with a generalized finite element method. International Journal for Numerical Methods in Engineering 102(3-4):180–201, 2015. doi:10.1002/nme.4703. [12] D. Ryckelynck. Hyper-reduction of mechanical models involving internal variables. International Journal for Numerical Methods in Engineering 77(1):75–89, 2009. doi:10.1002/nme.2406. 34 http://dx.doi.org/10.1016/j.jcp.2006.09.019 http://dx.doi.org/10.1016/j.compstruc.2016.01.003 http://dx.doi.org/10.1088/0965-0393/21/2/025014 http://dx.doi.org/10.1016/S0045-7825(03)00347-5 http://dx.doi.org/10.1002/nme.4703 http://dx.doi.org/10.1002/nme.2406 Acta Polytechnica CTU Proceedings 13:29–34, 2017 1 Introduction 2 Wang tiling concept 3 Synthesized enrichments fields 4 Numerical strategy 5 Results 6 Conclusions Acknowledgements References