Finite cell method for functionally graded materials based on V-models and homogenized microstructures

This paper proposes an extension of the finite cell method (FCM) to V-rep models, a novel geometric framework for volumetric representations. This combination of an embedded domain approach (FCM) and a new modeling framework (V-rep) forms the basis for an efficient and accurate simulation of mechanical artifacts, which are not only characterized by complex shapes but also by their non-standard interior structure. These types of objects gain more and more interest in the context of the new design opportunities opened by additive manufacturing, in particular when graded or micro-structured material is applied. Two different types of functionally graded materials (FGM) are considered: The first one, multi-material FGM is described using the inherent property of V-rep models to assign different properties throughout the interior of a domain. The second, single-material FGM—which is heterogeneously micro-structured—characterizes the effective material behavior of representative volume elements by homogenization and performs large-scale simulations using the embedded domain approach.


Introduction
Functionally graded materials (FGM) are a class of advanced materials that offer the possibility to exploit various desired physical properties within one component. This allows manufacturing 'high-performance' and 'multi-functional' artifacts which can resist physical exposures that could not be withstood by a single material [1]. The idea of combining different materials goes back more than 4000 years-the development of the composite bow-and has led to modern carbon fiber reinforced polymers. These composite materials change their material properties step-wise and are consequently prone to delamination. In FGM, on the other hand, material properties vary continuously inside the volume and avoid material interfaces [2]. Specific material properties are achieved by continuous changes in the microstructures, grain sizes, crystal structure, or composition of different materials such as metal, ceramics, polymers, or biological tissues [3,4]. Prototypes, especially for functionally graded microstructures, can be found in nature, such as in bones, seashells, skin, or wood [5] or obtained using topology optimization [6][7][8]. Fields of application are, amongst many others, corrosion resistance of chemically exposed components [9], bone-like lightweight porous medical implants [10], or heat resistance of load-bearing parts such as spacecraft thermal shielding, jet turbine blades, or nuclear reactors [3,11].
Additive manufacturing (AM) or 3D printing is a generic term for various production techniques in which an object is created by layer-wise material deposition. This allows the fabrication of objects of almost arbitrary shape. AM is the method of choice for the fabrication of FGM, as (i) it allows to resolve very fine structures, (ii) it can manufacture internal structures which could not be created with any other method, and (iii) the layerwise material deposition gives control over the composition of the processed material, as well as over the grain size [12,13]. With functionally graded additive manufacturing (FGAM), it is possible to create different single-and multi-material FGM [14]. Singlematerial FGM specimens consist only of one material that changes its properties due to an adaption of the microstructure, density, or grain size [15]. As AM allows the creation of free form structures, a single-material FGM in the form of a continuously changing microstructure can be fabricated with any printable material [16]. Multi-material FGMs, which blend two or more materials into each other within a volume, have recently been under intensive research [17]. A special focus was placed on metal-metal combinations, see e.g. [4], where steel and titanium-based combinations are investigated. More complex is the combination of materials of a different kind, such as ceramic-metal compositions [18]. However, these compositions might carry the most potential, as the underlying material properties are very distinct.
Material testing is the industry standard to determine the behavior of FGM components. Yet, physical test series are often elaborate and expensive. The goal of simulation supported development is therefore to reduce testing to only calibrating data for functionally graded materials and to then numerically analyze different shapes and compositions of artifacts. Within the scope of this paper, we present two distinct, novel approaches to perform numerical simulations on both, single-and multi-material FGMs, respectively. To this end, an analysis-suitable geometrical model needs to be provided which is naturally created with computer-aided design (CAD) and then transformed into a mesh. This transition process from CAD to an analysis-suitable mesh is error-prone. Depending on the quality of the model, manual work must be invested to heal the original geometry before mesh generation can be carried out successfully. Furthermore, the most used CAD representations, i.e. boundary representation (B-rep) or solid-based procedural models, are not well suited for an accurate description of FGM. B-rep models represent their volume implicitly by the boundary surfaces, which are modeled either with linear primitives (e.g. triangles and quads) or with trimmed spline patches [19]. Consequently, B-rep models offer no possibility to directly represent a heterogeneous material distribution inside the body. A workaround is to create vector functions that carry the material properties for each point. These functions can be classified into four different categories: (i) geometrically-independent, e.g., in Cartesian coordinates, (ii) distance-based, (iii) blending composition, and (iv) sweeping composition functions (for a detailed explanation refer to [20,21]). However, except (i), these functions only allow a smooth transition of material properties between the different surfaces, which is not suitable for all material distributions. Geometrically-independent functions, on the other hand, are cumbersome as they are not related to the object itself. CAD systems using solid-based procedural models follow the constructive solid geometry (CSG) idea [22]. Here, models are composed of simple primitives: spheres, cuboids, cylinders, etc. and more complex primitives: sweeps, lofts, extrusions, solid of revolution, etc. These primitives are combined with the classical Boolean operations: union, intersection, difference, negation, and their derivations: fillet, chamfer, holes, etc. Material properties can then easily be assigned to the respective primitives. Of course, this requires special treatment in regions with overlapping primitives [12]. Furthermore, as primitives are typically provided as implicit functions, they offer, similar to B-rep models, no possibility to further resolve the internal volume. Again, vector functions applied to the primitives are a possible workaround. Another possible geometrical representation is offered by spatial decomposition, such as voxelized models. Here, each voxel can carry its material properties. These voxel models mostly originate from CT scans (e.g. of bones) and provide only a coarse approximation while requiring an extensive amount of storage capacity. Nevertheless, voxel-based models have been used to resolve fine microstructures and quasi-continuous changes of the material properties [23,24].
Massarwi and Elber [25] recently proposed a novel volumetric representation technique (V-rep) for 3D models, which allow full control over the model's interior. V-reps consist of trimmed, trivariate B-spline patches which can be combined into V-models using Boolean operations. By extending the dimension of the control points, it is possible to assign material parameters directly to the model. This property can be used to model and simulate multi-material FGM. Potentially critical overlapping regions of the V-model are resolved by trimming the involved splines and creating new trivariate primitives for the respective overlapping volume. Due to the non-singularity of trivariate B-splines, V-models are predestined for a subsequent simulation using the isogeometric analysis (IGA) [26]. However, a direct application of IGA is often not feasible, since in overlapping regions the spline patches must be trimmed. Moreover, the respective spline parameterizationsi.e. the control point meshes, knot vectors, and polynomial degrees-do not coincide at adjacent faces. Hence, special techniques are required to glue them together, e.g. Mortar methods or T-splines [27,28]. By contrast, embedded domain methods require no special treatment of overlapping regions and pose far fewer requirements on the underlying geometric model.
Apart from the possibility to control the interior of the volume, which can be used to model multi-material FGM, the V-rep framework also offers the possibility to create single-material FGM, such as continuously changing microstructures. Although easy to fabricate with AM, these multiscale structures are critical from a simulation point of view. Due to the complexity of the underlying CAD models, the meshing becomes difficult. Additionally, attempts to resolve the structure sufficiently accurate may result in overrefined meshes -which in turn lead to an additional but unnecessary computational effort. This is where numerical homogenization provides an efficient tool to estimate an overall mechanical behavior of such structures. The basic idea of homogenization is to define a representative volume element (RVE) which is sufficiently large to represent the overall material behavior in the specific region [29][30][31]. In the case of periodic microstructures, a unit cell can be extracted for further material characterization. Periodic boundary conditions are then applied at their boundaries, which leads to the best possible estimate of the effective behavior [32,33]. The resulting material characterization can then be used to simulate a complete structure under complex loading. The computational cost is reduced considerably by 'smearing out' the detailed complex geometrical features of a microstructure and expressing them in terms of the effective behavior. Still, on the microscopic level of the RVE, the structure needs to be fully resolved in a boundary conforming fashion to account for all geometrical details. Here, embedded domain methods offer an elegant and reliable alternative over classical FEA also for non-periodic AM structures [34].
In this contribution, three novel methodologies are introduced: • The FCM is extended to V-models as a new CAD representation form.
• Based on the trivariate spline description of the V-models, a method for the simulation of multi-material FGM is introduced.
• Finally, a distinct approach is proposed that allows numerical analyses on large-scale continuously changing microstructures -i.e. single-material FGM-using homogenization.
The paper is structured as follows: "Finite cell method" and "Volumetric representation" sections provide a brief overview over the FCM and V-reps, respectively. The methodologies for the simulation of V-reps, single-and multi-material FGM are described in "Extension of the FCM to single-and multi-material FGMs" section. "Numerical examples" section presents and discusses several numerical examples before conclusions are drawn in "Conclusion" section.

Finite cell method
In the following, the basic concepts of the finite cell method are briefly summarized for linear elasticity. A detailed description of the FCM can, e.g., be found in [42]. The FCM embeds a physical domain phy into a fictitious domain fict forming an extended domain ∪ , as illustrated in Fig. 1 for two dimensions. The weak form of the equilibrium equation for the extended domain ∪ reads as follows where u is the unknown deflection, v is a test function, L is the kinematic differential operator and C is the constitutive material tensor. The body load and the prescribed tractions on the Neumann boundary N are denoted by b and t, respectively. To resolve the complex domain correctly, an indicator function α(x) is introduced which weights the material tensor C In the limit q = ∞, Eq. 1 recovers the standard weak form for phy . In a finite elementlike discretization, however, it leads to ill-conditioned systems. This can be avoided by choosing a finite q (in practice q = 6 . . . 10) in combination with a suitable preconditioning and/or orthogonalization of the shape functions [56]. This choice introduces a modeling error [57] but limits the conditioning number of the stiffness matrix. Further improvement on the conditioning can be obtained using preconditioning, orthogonalization of shape functions, and/or the increase of continuity between the cut cells [58]. The extended domain ∪ is of simple shape and can be easily meshed into regular cells, e.g. rectangles in 2D and cuboids in 3D, respectively. These cells can be locally refined into sub-cells or with respect to the order of the shape function [59,60].

Geometry treatment
The FCM resolves the physical domain phy (i.e. the geometric model) by the discontinuous scalar field α(x), which is then queried during the integration of the system matrices and load vectors. Consequently, the resolution of the geometry's complexity is shifted from the discretization (conforming meshing) to the integration level. The only information the FCM requires from the geometry is an unambiguous point inclusion statement, i.e. it must be possible to decide for any point x whether x ∈ phy or x ∈ fict . Due to the discontinuity of α(x) on cut cells, the integration needs to be carried out using special quadrature rules. Common variants are composed integration on a space-tree reconstruction (see Fig. 1), smart quadtree/octree, or moment fitting [61][62][63]. Another approach uses dimensional reduction, i.e. the integration is not performed over the entire domain, but only along the boundary [64].

Boundary conditions
As the boundary of the physical domain phy typically does not coincide with the edges/faces of the finite cells, essential (Dirichlet) boundary conditions need to be applied in a weak sense. For this, several methods have been investigated-such as the Nitsche method, Lagrange multipliers, and the penalty method [65][66][67][68]. Natural (Neumann) boundary conditions are applied on N following Eq. (1). Homogeneous natural boundary conditions are automatically resolved by α(x) ≈ 0. Inhomogeneous natural and essential boundary conditions require an explicit integrable boundary description, which is either provided by the geometrical model or extracted directly from the volume using, e.g., the marching cubes algorithm, see e.g. [69].

Volumetric representation
The V-rep framework [25] provides methods and algorithms for the construction of Vmodels by combining simple (e.g. cylinder, sphere, etc.) or complex primitives (e.g. ruled primitives or solids of revolution) with the Boolean operations, thus following the idea of constructive solid modeling. Furthermore, it is possible to migrate spline-based B-rep models to V-rep models. The V-rep framework is embedded in the Irit geometry library [70], developed by Elber et al. Irit provides a vast amount of various geometric modeling and analysis functionalities, and it can be accessed as a C(++) library, via a scripting language, or graphically with the GuIrit CAD environment [71].

Trivariate B-splines
A trivariate B-spline is a parametric function that allows to span a volume over a threedimensional parameter space. It is typically represented as follows

V-rep primitives
Apart from the trivial case of a cuboid, the V-rep framework offers a variety of both highlevel and simple primitives. Implemented are several high-level primitive constructors, all of which yield one single trivariate patch (see Fig. 2 Simple primitives-such as spheres, cylinders, tori, and cones -can not be represented by a single trivariate patch without introducing singularities (e.g. at the mid axis of a sphere, the Jacobi matrix vanishes det(J V (r = 0)) = 0.) To this end, singular primitives are composed of several non-singular trivariate patches (see Fig. 3).

V-model construction
A trivariate B-spline is limited to a cuboid topology. To represent general volumetric shapes, so-called 3-manifold V-cells ν i C are introduced which correspond to trimmed . . , n}. These V-cells originate firstly from the primitives that constitute the CAD model. New Vcells occur due to the combination of the Boolean operations in the regions of overlapping primitives. Here, trivariate B-splines are trimmed at intersecting surfaces, and, depending on the Boolean operation, the intersection volume is then remodeled from the trimming surfaces using the Boolean constructor (see Fig. 4). Consequently, the V-cells of a V-model . . , n} and the parametrizations of the new 'intersection' V-cells are different from their parent primitives. This makes the use of IGA more complex. V-cells store additional topological and adjacency information, which allows an efficient model inquiry. Adjacent V-cells share common trimming/boundary surfaces. Analogously to B-Rep, the boundary of the V-model ∂V m forms a closed 2-manifold.

Extension of the FCM to V-reps
In the context of the finite cell method, at first, without considering functionally graded materials, the V-model only needs to provide a point inclusion test. To this end, an inverse mapping is carried out on each V-cell.
As splines can generally not be inverted analytically, the corresponding parameter position u must be determined iteratively using the Newton-Raphson algorithm. Yet, one should note that-since the splines are regular, i.e. the Jacobian never vanishes-a solution is always unique if one exists. In the case x ∩ ν i C = ∅ a parameter position can be found in the V-cell ν i C and the respective point x is inside the V-model. The number of required Newton-Raphson iterations for the inverse mapping can be substantially decreased providing a good guess as an initial value. This is efficiently exploited by the finite cell method as, due to the Cartesian grid-based data structure, consecutive integration points are very often geometrically adjacent. Therefore, the last inner point on each V-cell is cached and used as an initial guess for the next query.
Since, the underlying Irit library [70] offers already a robust point inclusion test, the extension of the FCM to V-models is straight-forward. It is noteworthy that-in contrast to IGA-trimmed splines, as well as non-coinciding spline parameters at adjacent faces, require no special treatment since the adaptive quadrature rules automatically recover the actual shape of the geometry.

Extension of the FCM to single-and multi-material FGMs
The V-Rep framework provides two different ways to realize functionally graded materials which can be produced by additive manufacturing techniques: (a) the material properties can either be encoded directly into the volume of the V-cells (see "Analyzing multimaterial FGM with the finite cell method" section), which is perfectly suited to model multi-material FGMs, or (b) an FGM can be created in a constructive manner in the form of a continuously changing mircostructures, which corresponds to single-material FGMs (see "Analyzing multi-material FGM with the finite cell method" section).

Analyzing multi-material FGM with the finite cell method
A simulation of multi-material FGM using the FCM requires-apart from the pointinclusion statement-also the corresponding material properties at any location. To this end, the spline-based description of the V-cells-as the smallest, non-intersecting building blocks-is extended to also carry material information. V-Rep material representation Material properties such as Young's modulus, Poisson's ratio, thermal conductivity, density, etc. can easily be represented on the V-cells by simply extending the dimension of the control points R 3+s , with s > 0 being the additional material parameters (see Eq. (3)). Consequently, evaluating the V-cell yields, in addition to the geometric coordinates, also the respective material values As an example, consider a control point that carries additional material properties for the Young's modulus E, Poisson's ratio ν, and thermal conductivity κ as needed for Example 3.3.2: The material properties of a V-cell, created from the overlap of two or more trivariate B-splines carrying different material information, require additional handling. Either one of the initial trivariate B-spline can be set prevailing and, thus, its properties are inherited to the V-cell, or the material properties are interpolated by some sort of blending scheme. For detailed information refer to [25]. Spline based material approximation Inside a patch, splines are typical of higher continuity, which renders them perfectly suitable for modeling smooth geometries. However, this restricts the material function to be of the same continuity. A remedy to also represent C 0 or discontinuous material distributions is given by knot-insertion, as the continuity depends on the multiplicity of the knots C p−m , where p is the polynomial degree and m the number of multiple knots. Naturally, knot-insertion also reduces the potential continuity of the geometry. However, the original higher continuity is preserved in a geometrical sense. 1 Hence, the model keeps its geometrical shape, whereas the material is allowed to One-dimensional least squares approximation of a hypothetical sinusoidal material function m σ (x) = sin(2πn p x), with n p = 2.5 being the number of periods, yields the material 'coordinates' μ σ i . Note that the rather large deviation between the curves comes from the fact that the locations (i.e. x−coordinates) of the material control points as well as the knot vector are fixed have material kinks, or even to be discontinuous. Nevertheless, due to the global influence of the position and multiplicity of the knots, splines are not the method of choice to represent highly discontinuous material distributions, as e.g. underlying voxel data provided by CT-scans.
Given a sufficiently smooth material distribution, the material 'coordinates' of the control points can be obtained using least-squares approximation (see Fig. 5). For each material property, the least-squares problem reads where n LS is the number of sample points and μ σ = μ σ i,j,k ∈ R l·n·m are the minimization variables (see Eq. (3) for l, m, n). The least squares problem is then solved for each material function f σ m and the respective material 'coordinate' contains the spline basis functions. The sample points are evaluated in the parameter space Consequently, the material function needs to be evaluated in the same space (see Eq. (3))

Analyzing single-material FGM with the finite cell method
Single-material FGM structures change their material properties due to adaptions in the microstructure, density, grain size, etc. A prominent example in nature is the trabecular bone, where the size and alignment of thin rods and plates of bone tissue create stiffness trajectories that follow the principal stresses for the most common load cases [73]. Today, additive manufacturing (AM) offers the possibility to create similarly complex structures. To this end, AM uses porous infill structures to support the outer hull. How- 6 Functionally graded microstructure: b Three different anisotropic tiles, with a changing stiffer direction, are used to tile a a rotating ruled volume. c The entire resulting microstructure exhibits a continuously changing anisotropic stiffness ever, this infill is typically a repetitive lattice and is either not taken into account for the load transfer, or it is assumed to be isotropic [74]. Nonetheless, recent approaches in the field of topology optimization try to exploit the contribution of the infill to the load transfer [75]. Problem-fitted complex 3D anisotropic microstructures can reduce the printing time and material consumption substantially and at the same time improve the load-carrying properties and buckling behavior. Gradually changing microstructure The V-rep framework offers the possibility to create complex aniso-tropic microstructures with its tiling operation. Hereby, copies of a unit structure are consecutively created inside a base volume. Following the shape of the base volume and by using layers of different unit cells, a complex constructive FGM can be created. As the resulting microstructure is composed of several V-cells, it is again a Vmodel. Naturally, each V-cell can again represent a heterogeneous material distribution within its volume. Even for complex tile-based structures, like the example shown in Fig. 6c, the point inclusion test can be carried out by inverse mapping as described in "Extension of the FCM to V-reps" section. Yet in case of single-material structures, it turns out that a conversion into an auxiliary B-rep and a consecutive ray tracing based test (see [76]) is computationally more efficient. In our implementation, the B-rep surface is subdivided into a fine triangular mesh and stored in a kd−tree [77]. Certainly, in contrast to the inverse mapping on trivariate B-splines, the surface triangulation causes a further approximation error, which can yet be controlled by refining the surface subdivision.
Simulation of large-scale single-material FGM with the FCM Detailed geometrical features of microstructures require a fine numerical resolution to achieve reliable simulation results. Hence, large-scale structures necessitate the use of high-performance computers, or might even not be computable at all. To reduce the computational cost and thus allowing the computation of large microstructures, a numerical homogenization can be used to evaluate a macroscopic mechanical behavior under specified loadings. The basic idea of this method is to approximate the solution of a macroscopic boundary value problem by solving less demanding microscopic problems [31]. This idea relies on the existence of a representative volume element (RVE), which is a microstructural domain that is large enough to represent macroscopic behavior and small enough to ensure the scale separation. The mechanical quantities can then be transferred from the micro-to the macro-scale by using the Hill-Mandel condition, which is also called 'macro-homogeneity condition'. This mean-field numerical homogenization provides reliable estimates for the effective mechanical behavior if appropriate boundary conditions are chosen. For the herein considered microstructures that are created with Irit's tiling operation, periodic boundary conditions provide the best effective material properties.
Certainly, a functionally graded microstructure cannot be represented by one single RVE. However, since the parameter-based construction plan of the FGM is known apriori, it is sufficient to compute the effective material tensors C * i for several 'representative' RVEs (see Fig. 7). At any realization in-between, the material is then be interpolated from corresponding adjacent representative tensors C * i . For the microstructures, considered in this paper, the RVEs correspond to the constituting unit tiles, which can have different properties, for instance, a stiffer direction, a rotation around some axis, or a material composition. All these properties are defined using suitable construction parameters. The following approach then allows the efficient computation of large-scale functionally graded microstructures with the FCM: • Using the parametric description of the microstructure, several representative unit tiles are selected in different configurations. • For the unit tiles, the effective material tensors C * Ti are computed with a combination of numerical homogenization and the finite cell method [34] and stored in a look-up

Numerical examples
To demonstrate the variety of simulatable functionally graded materials using a combination of V-reps and the FCM, five examples are presented. The first example serves as a verification of the extension of the FCM to multi-material FGMs. To this end, a linear elastic simulation of a simple cuboid with a prescribed material distribution is performed. The second example, a coupled heat, thermo-elastic simulation of a curved thermal protection tile, underlines the applicability to examples of engineering relevance. The third example shows a simulation of a fully resolved single-material FGM-i.e. a continuously changing microstructure. In the fourth example, the underlying tiles of the third example are evaluated in terms of a homogenization, which are then used in the fifth example to perform a simulation on a large-scale homogenized single-material FGM.

Example 1: Cuboid with sinusoidal material distribution
As a benchmark problem, the cuboid with varying material distribution in z−direction is chosen. 2 The cuboid is a trivariate B-spline and is created with GuIrit [71]. As the spline basis functions are initially linear in each direction, they are not able to represent the material function E(z). For this reason, a degree elevation to r = 3 and subsequent multiple knot insertions in z−direction were carried out, yielding a knot-vector of W = [0, 0, 0, 0, 0.2, 0.4, 0.6, 0.8, 1, 1, 1, 1]. The control points in z−direction are depicted in Fig. 8a. The cuboid has assigned a constant Poisson ratio of ν = 0.3. The functionally graded Young's modulus is given as an analytical function E(z) = 10 6 + 5 · 10 4 · sin(zπ) .
The material 'coordinates' μ E i of the control points are computed using least squares with n LS = 100 sample points, according to Eq. To prove the validity of the FCM for multi-material FGM, a convergence study is carried out. The polynomial degrees of the Legendre ansatz function are increased from p = 1 . . . 8 and the corresponding strain energies are compared to a reference solution U ref , which was computed by a boundary-conforming p−FEM analysis. To minimize integration errors, a composed integration is used, which can exactly recover the exact shape of the cuboid -similar to the smart octree [61]. To compare the convergence behavior of the FCM with the standard FEM, two additional convergence studies using h−refinement are carried out on boundary conforming FEM discretizations, with polynomial degrees of p = 1 and p = 2, respectively.
As depicted in Fig. 9b, the FCM shows a pre-asymptotic exponential convergence until it reaches the numerical precision of the underlying Irit library at p = 4, whereas the h−studies show algebraic convergence-as expected. 3 Obviously, in terms of degrees of freedom, the FCM outperforms classical h−versions. Figure 10 shows the displacements and the von Mises stresses of the deformed cuboid. As expected, the regions of lower stiffness are undergoing a larger deformation.

Example 2: Curved thermal shielding tile
The second example demonstrates the applicability for practical applications. To this end, three curved thermal shielding tiles are simulated. Such tiles are needed for hightemperature applications, such as re-entrance shielding for spacecraft or the inner coating of fusion power plants. The tiles consist of a load-carrying zone made of titanium Ti and an insulating zone made of porous silica SiO 2 with a porosity of 70%. Both materials have similar melting points of Ti = 1.668 • C for titanium and SiO 2 = 1.710 • C for silica, which allows a fabrication with additive manufacturing using e.g. powder bed laser melting.
Special focus is laid on the continuity of the transition zone between these materials. The first discontinuous tile consists of two distinct domains where both domains are assumed to be homogeneous titanium and silica, respectively, i.e. there is no transition zone. Hence, the first tile is not a FGM, but a composite material. The material is changed C 0 −continuously in the second tile, and C 1 −continuously in the third tile. To evaluate the stresses under a heat load, a coupled simulation is carried out. An initial thermal simulation provides the temperature distribution, which is then used to apply thermal strains for the subsequent thermo-elastic simulation. Consequently, the model will deform due to the different thermal expansion ratios. This deformation is, however, hindered by the different Young's moduli in the transition zone, then leading to internal stresses.
The underlying V-model consists of one V-cell and was generated by extruding a curved two-dimensional B-spline surface 5 cm in z−direction. The control point mesh of the curved surface is defined as follows 4 The resulting material distributions are depicted in Fig. 12 exemplary for the Young's modulus. The other material properties are distributed similarly. The parameters for the B-splines were chosen such that the integral of the material over the thickness is equal for all three tiles. Figure 11 shows the outer dimensions of the tiles in cm.
To perform the coupled simulation, four different material parameters are required for both materials (see Table 1). The properties were taken from AZO Materials and averaged if necessary [78]. Due to the porosity of the silica, the respective Young's modulus E SiO 2 and the thermal conductivity κ SiO 2 must be adapted. This is implemented with the Gibson-Ashby criteria, which provide simple formulas to estimate the properties based of the porosity [79,80] where φ is the porosity (in this example φ = 0.7) and κ r and E r are the weighting factors for the thermal conductivity and Young's modulus, respectively. In contrast, the Poisson's ratio ν SiO 2 and the thermal expansion α SiO 2 require no adjustment [81].  The simulation uses 16 × 23 × 9 finite cells with a polynomial degree of p = 3 and an integration subdivision depth of n = 3. For the preceding heat simulations, Dirichlet boundary conditions are applied with a prescribed heat of 1000 • C on the top surface and 20 • C on the bottom surface. The resulting temperature inside the tiles is then transferred as a body strain to perform a thermo-elastic simulation. Additionally, the tiles are clamped at the bottom surface. Since the higher-order shape functions are not able to represent jumps in the material distribution, the simulations of the tile with the discontinuous material distribution are carried out on two separate meshes -one for each domain -, which are 'glued' together in a weak sense along their coupling surface, using the penalty method [82]. Both meshes are equally discretized with 16 × 23 × 9 finite cells. Thus, on each mesh, only material jumps from the physical into the fictitious domain appear, which can be decently represented by the FCM.
To resolve the critical regions, the finite cell mesh is refined using h−refinement. One h−refinement step yields eight subcells for each (refined) finite cell, which can then be further refined in a subsequent refinement step. Certainly, this kind of refinement introduces hanging nodes between refined and unrefined cells. The resulting incompatibilities between the shape functions are resolved by the multi-level hp-method [60]. For the discontinuous tile, both meshes are refined twice towards the coupling surface -meaning the respective finite cells are refined with a minimum of 15 and a maximum of 64 subcells. For the continuous tiles, all finite cells that are intersecting the respective transition zones are refined once (see Fig. 13).
To visualize the results inside the tiles, a cut through the model is investigated at x = 5 cm. Figure 14 shows the temperature distribution and displacements of the C 0 −continuous tile. The temperature and the displacement distributions are almost identical for all tiles. More relevant are the stress distributions. As can be seen in Fig. 15, a  Fig. 13 Discretizations of a the discontinuous tile: The mesh is refined twice around the coupling surface (yellow), which divides the upper (light blue) and lower domain (purple). b The C 0 −continuous tile: The FCM mesh is refined once in the transition zone (cells in blue are unrefined, and cells in red are refined once). The grey mesh in the background corresponds to the octree for the integration. The C 1 −continuous tile is meshed and refined analogously The discontinuous material distribution yields a C 0 −continuous heat and displacement distribution, which then entails a discontinuous stress distribution with a maximum peak at the interface region. This is critical as it will potentially cause delamination. The C 0 −continuous material distribution, on the other hand, ensures a continuous and much smaller stress distribution throughout the entire domain. This effect can be augmented further by using a C 1 −continuous material distribution. Continuous materials, on the other hand, involve a larger heat flux. For the 1D case, the thermal resistance is reduced to approximately 86% for the C 0 −continuous and approximately 75% for the C 1 −continuous material with respect to the discontinuous material distribution.

Example 3: Anisotropic microstructure
The third example addresses the second kind of functionally graded materials-namely single-material FGM. For this, a linear-elastic simulation of the continuously changing microstructure, depicted in Fig. 6, is carried out. It resembles a porous, foam-like microstructure stiffened by an outer shell. To generate this model, a continuously changing microstructure is created with GuIrit. Different unit tiles-each composed of seven trivariate B-splines-are used to tile a parametrically described ruled body (see Fig. 2). The unit tiles have a growing stiffness from bottom to top, realized by an increasing diameter of the rod in x−direction. 5 The resulting microstructure consists of 6 × 6 × 9 unit tiles and an overall number of 2268 trivariate B-splines, or V-cells. A direct simulation on this Vmodel leads to unreasonably high runtimes, due to the complexity of the inverse mapping. However, since in this example, the FGM is not modeled within the individual V-cells, but as a single-material continuously changing microstructure, it is possible to carry out a simulation significantly faster on an auxiliary B-rep model. To this end, the B-spline surfaces  Figures 19 and 20 show the displacements and the von Mises stresses. Certainly, such a fully resolved simulation is slower than the numerical homogenization presented in "Analyzing single-material FGM with the finite cell method" section-especially because homogenization in the linear case allows the creation of a lookup table. However, the discussed fully resolved model can be used to verify the homogenization. This is addressed in the following Examples 3.3.4, and 3.3.5. Note, since the shape functions are badly suited to represent holes inside one finite cell, meaning 'material-void-material' [83], the microstructure needs to be resolved with many finite cells. A remedy can be local enrichment as presented in [84].

Example 4: Material characterization database for unit tiles
A fully resolved numerical simulation of a microstructure -as presented in Example 3.3.3-is computationally very demanding in both memory consumption and simulation time. For large-scale microstructures (as in Example 3.3.5), fully-resolved computations need to be carried out on a high-performance computer, or might even be not applicable at all. A remedy is offered by homogenization. As explained in "Analyzing single-material FGM with the finite cell method" section, for a functionally graded microstructure it is sufficient to compute the effective material tensors C * Ti only for several representative unit tiles, and interpolate the material properties in-between, according to the parametrization of microstructure.
Two parameters are used to characterize the unit tiles in the Examples 3.  For the homogenization simulations, the material of the microstructure is considered to be steel with a Young's modulus of E = 210 GPa, and a Poisson's ratio of ν = 0.3. Each tile is discretized with 11 × 11 × 11 finite cells of polynomial degree p = 5. For the domain integration, the moment-fitting approach [62] with the depth of an underlying octree of d = 6 is chosen. As the structures under consideration are, in good approximation, geometrically periodic, periodic boundary conditions are the natural choice for transferring the macroscopic quantities to the microscopic unit cells. Figure 21 shows the displacement fields under shear load for the unit tiles in the unrotated configuration. The resulting homogenized material tensors for the tiles 1, 2 and 3 are summarized in Eqs. (19), (20) and (21), respectively. One can identify different material behaviors, which is expected due to the geometrical features of the respective unit tiles. The orientation and the thickness of the rods has an important effect on the final material behavior. Tile 1 shows a cubic macroscopic material symmetry with three independent elasticity coefficients [85], namely C 11 , C 12 and C 44 Due to the stiffer direction in x−direction, tile 2 and 3 show a tetragonal effective material symmetry with C 11 , C 22 , C 44 , C 55 , C 12 and C 23 as independent entries: The material tensors C Ti for the second changing parameter-the rotation around the z−axis-can be computed by a coordinate transformation, and thus require no homogenization simulations. The Bond-Transformation matrices [86] can be used to rotate the effective elasticity tensor by a matrix-matrix multiplication. Assume the following ordering of the macroscopic stresses σ M ij and strains ε M ij in the Voigt notation Then, the transformation of the effective elastic tensor reads as follows where C * is the effective elasticity tensor, C is the effective elasticity tensor in rotated coordinates, and M and N are the Bond-stress and the Bond-strain transformation matrices, respectively. For the rotation around the z−axis, the Bond strain and stress matrices are defined as follows In "Appendix" section presents the respective independent material tensor entries C ii of the three unit tiles for arbitrary rotations around the z−axis, following Eq. (23). Given a set of different (an-)isotropic unit tiles that can be used to construct such microstructures, it is possible to create a look-up table of homogenized materials, which can then be used to simulate different macroscopic load cases. Table 2 is a snippet of such a look-up table, and it shows the effective elasticity tensors for the two varying material properties. The material properties in-between can be interpolated. This Table 2, will be used in the following Example 3.3.5 to compute a large-scale microstructure with interpolated homogenized material properties.

Example 5: Homogenized microstructure
Consider the model of Example 3.3.3 to be a part of a larger structure (see Fig. 22). Based on the material database for the homogenized unit tiles (see Table 2) it is possible to simulate such a structure with a homogenized material. Similar to Example 3.3.3, the corresponding geometric parts are modeled as B-rep models. For the simulation, the model is subdivided into an outer shell and an infill. The shell is considered to be of solid isotropic material, whereas the infill is a homogenized microstructure which continuously changes the two known properties: the rotation angle ψ around the z−axis varies from 0 • at the bottom to 90 • at the top and the thickness of the rod Ø increases from the center z−axis of the infill (Ø = 0.  unit tiles was carried out with periodic boundary conditions, the behavior at the interface between shell and infill is not captured precisely. However, the affected domain is small compared to the overall structure, thus, the introduced error is negligible. If, however, the microscopic stress state at the transition from the micro-tiles to the shell is of interest, then a geometrically resolving simulation as in Example 3.3.3 can be performed. A total of 13 independent material coefficients are required to evaluate the material tensor of the continuously changing microstructure. To this end, the material coefficients that were computed in Example 3.3.4 and that are stored in a look-up table (see Table 2) are interpolated using spline fitting. Figure 23 exemplary shows the interpolation for the material coefficients C 11 and C 22 of the homogenized material tensor shown in Eq. (26).   Figure 24 shows the displacements in x−direction and von Mises stresses of the structure under uni-axial compression z−direction. The load is mainly transferred through the stiffer shell, yet the contribution of the infill cannot be neglected. Due to the uni-axial compression, the rotation angle ψ of the microstructure has only little influence. The thickness of the rod Ø, on the other hand, can be deducted directly from the stress field of the infill. It should be noted that a geometrical change does not influence the overall workflow. Even a topological change does not lead to a re-meshing as it would be required for a simulation with classical FEM or IGA. To illustrate such a topological change, a hole is drilled through the structure (see Fig. 25). In the context of the FCM, a cylinder is subtracted with a Boolean difference. As can be seen, the infill contributes less to the load transfer, and high stress concentrations appear at the walls of the hole.

Conclusions
In this paper, three novel methodologies were presented: (a) At first, the FCM was extended to V-models, as novel CAD representation form. As V-rep is based on a trivariate spline-formulation, the inversion-that is necessary for the point inclusion testturns out to be costly, in particular in cases where due to the geometric complexity of the model a large number of integration points has to be used. In these cases the definition of an auxiliary B-rep model using ray tracing for the point membership test turns out to improve the computational performance significantly. (b) Secondly, the FCM was extended to multi-material FGM. For this, the dimension of the control points of the V-cells was increased to carry material information, as well. During the integration-apart from the point-inclusion test-also the material properties are retrieved. The spline-based description of the V-cells renders the V-rep framework perfectly suitable to model smooth material distributions. Yet, also rapidly changing materials can be represented using knot-insertion.
(c) Finally, an efficient method for the simulation of large-scale single-material FGM-in this case continuously changing microstructures-was presented. Using the parametric description of the microstructures, representative unit tiles can be selected on which homogenization simulations provide effective material properties. Material properties for adjacent parameter sets are then interpolated using these values. Although this approach allows the efficient simulation of large-scale microstructures, two problems arise. Firstly, depending on the complexity of the microstructure and the amount of varying geometrical features, the number of representative unit tiles might become large. Since for each of these unit tiles, an individual homogenization simulation needs to be carried out, such structures can become demanding in memory consumption as well as in computational time. And secondly, the homogenization simulations with periodic boundary conditions provide only precise results within the microstructure, yet not at the interface to another material or a free surface. However, provided this interface or surface area is small compared to the overall domain and considering that such kind of boundary layer effects usually vanish rapidly away from the interface, the error is not dominant. A rotation of tile 1 around the z−axis does not influence the third, fifth, and sixth columns, neither on the respective rows of the effective tensor. The coefficient C 11 equals C 22 due to the geometrical symmetry in x− and y−direction. C 14 and C 24 are of equal magnitude but have opposite signs. Figure 27 shows the remaining independent material constants with respect to the rotational angle. The results of the numerical simulation at 45 • are indicated with red crosses. For tile 2, only the coefficient C 33 which corresponds to the stiffness in z−direction remains unchanged under rotation around the z−axis. All other entries are affected by the altered symmetry. Considering a rotation angle of 90 • , it is noteworthy that the coefficients C 11 and C 22 are switched with regard to the initial position. The same holds for the coefficient pairs C 55 -C 66 , and C 13 -C 23 . The rest of the independent material parameters are depicted in Fig. 28. Again, the results of the numerical simulation at 45 • are marked with red crosses. Tile 3 exhibits similar material symmetries as the second tile. Figure 29 shows the material coefficients. Again, the results of the numerical simulation at 45 • are marked with red crosses.