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 objects gain more and more interest in the context of the new design opportunities opened by additive manufacturing, mainly 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 V-rep models' inherent property 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 advanced materials that offer the possibility to exploit various desired physical properties within one component. FGMs allow manufacturing 'high-performance' and 'multi-functional' artifacts which can resist environmental 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. On the other hand, in FGM, 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 material deposition allows the fabrication of objects of almost arbitrary shape. AM is the method of choice for the manufacturing of FGM, as it can (i) resolve tiny structures, (ii) manufacture internal structures which could not be created with any other method, and (iii) the layer-wise 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]. Single-material 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 particular 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. Therefore, the goal of simulation supported development is to reduce testing to only calibrating data for functionally graded materials and then numerically analyze different shapes and compositions of artifacts. Within this paper's scope, 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 model's quality, 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 FGM description. B-rep models represent their volume implicitly by the boundary surfaces, which are modeled either with linear primitives (e.g., triangles and quads) or trimmed spline patches [19]. Consequently, B-rep models offer no possibility to represent a heterogeneous material distribution inside the body directly. 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. On the other hand, geometrically-independent functions 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 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 control points' dimension, 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 nonsingularity of trivariate B-splines, V-models are predestined for a subsequent simulation with 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 parameterizations -i.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 of controlling the interior of the volume, which can be used to model multimaterial FGM, the V-rep framework also offers the option 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 over-refined meshes, leading 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 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: Sections 2.1 and 2.2 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 Section 3.2. Section 3.3 presents and discusses several numerical examples before conclusions are drawn in Section 4.

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 Ω f ict forming an extended domain Ω ∪ , as illustrated in Figure 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 = ∞, Equation 1 recovers the standard weak form for Ω phy . In a finite element-like 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 complex geometry 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 ∈ Ω f ict . 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 Figure 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 Equation (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 V-models 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. 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 three-dimensional parameter space. It is typically represented as follows where V (u) is a point inside the volume and u = (u, v, w) T the corresponding three-dimensional parameter position in the parameter space u ∈ U × V × W ⊆ R 3 . B i,p denotes the i th one-dimensional B-spline basis function of polynomial degree p and P i,j,k ∈ R k are the l × m × n control points. The dimension of the control points is k = 3 + s, where k = 3 corresponds to the three geometric coordinates x T = [x, y, z]. Further information can be represented by additional dimensions s > 0.

V-rep primitives
Apart from the trivial case of a cuboid, the V-rep framework offers various high-level and simple primitives. Implemented are several high-level primitive constructors, all of which yield one trivariate patch (see Figure 2): 1. Extrusion: A surface is extruded along a vector.
2. Ruled solid: A volume is defined as a linear interpolation between two surfaces.
3. Solid of revolution: A volume is constructed by rotating a surface around an axis.
4. Boolean sum: A volume is created from six boundary surfaces [72].

Sweep/Loft: A sweep or loft interpolates several surfaces along a sweeping path.
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 Figure 3).

V-model construction
A trivariate B-spline is limited to a cuboid topology. In order to represent general volumetric shapes, so-called .., n}. These V-cells originate firstly from the primitives that constitute the CAD model. New V-cells occur due to the combination of the Boolean operations in the regions of overlapping primitives. Here, trivariate B-splines are trimmed at intersecting surfaces. Depending on the Boolean operation, the intersection volume is then remodeled from the trimming surfaces using the Boolean constructor (see Figure 4). Consequently, the V-cells of a V-model are non-intersecting  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 and noncoinciding 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 Section 3.2.1), which is ideally suited to model multi-material FGMs, or (b) an FGM can be created in a constructive manner in the form of a continuously changing microstructures, which corresponds to single-material FGMs (see Section 3.2.2).

Analyzing multi-material FGM with the finite cell method
A simulation of multi-material FGM using the FCM requires -apart from the point-inclusion statement -also the corresponding material properties. To this end, the spline-based description of the V-cells -as the smallest, non-intersecting building blocks -is extended to carry additional 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 Equation (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: P T i,j,k = [x, y, z, E, ν, κ] i,j,k . 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. Thus, its properties are inherited to the V-cell. Or some sort of blending scheme interpolates the material properties. 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 have material kinks, or even to be discontinuous. Nevertheless, due to the global influence of the knots' position and multiplicity, 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 Figure 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 Equation (3) for l, m, n). The least squares problem is then solved for each material function f σ m and the respective material 'coordinate' µ σ , σ ∈ [1, s] of the control mesh P i,j,k = x, y, z, µ 1 , ..., µ σ , ..., µ s T i,j,k . Matrix A ∈ R ν×(l·n·m) 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 Equation (3) Figure 5: 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.

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. 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. However, this infill is typically a repetitive lattice and is either not taken into account for the load transfer, or is assumed to be isotropic [74]. Nonetheless, recent approaches in 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 anisotropic microstructures with its tiling operation. With this, 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 built. As the resulting microstructure is composed of several V-cells, it is again a V-model. 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 Figure 6c, the point inclusion test can be carried out by inverse mapping, as described in Section 3.1. Yet in the 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 high-performance computers, or might even not be computable at all. In order 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. This method's basic idea is to approximate the solution of a macroscopic boundary value problem by solving less complicated microscopic problems [31]. This idea relies on the existence of a representative volume element (RVE). This microstructural domain 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 using the Hill-Mandel condition, 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 Figure 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 * T i are computed with a combination of numerical homogenization and the finite cell method [34] and stored in a look-up table (see Example 3

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 Figure 8a. The cuboid has assigned a constant Poisson ratio of ν = 0.3. The functionally graded Youngs 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 Equation (6)  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 accurately recover the cuboid's exact shape -similar to the smart octree [61]. In order 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 Figure 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.

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 high-temperature 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 Θ T i = 1.668 • C for titanium and Θ SiO2 = 1.710 • C for silica, which allows a fabrication with additive manufacturing using e.g. powder bed laser melting.
Particular 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 an 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 thermoelastic simulation. Consequently, the model will deform due to the different thermal expansion ratios. However, this deformation is 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 Figure 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 SiO2 and the thermal conductivity κ SiO2 must be adapted. This is implemented with the Gibson-Ashby criteria, which provide simple formulas to estimate the properties based 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 ν SiO2 and the thermal expansion α SiO2 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. In order 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. This kind of refinement introduces hanging nodes between refined and unrefined cells. The multi-level hp-method [60] resolves the resulting incompatibilities between the shape functions. For the discontinuous tile, both meshes are refined twice towards the coupling surface -meaning the individual 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 Figure 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  Figure 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. be seen in Figure 15, a stress concentration occurs at the coupling surface of the discontinuous tile. Figures 16  and 17 plot the temperature distribution, displacements, and stresses over the height at x = 5 cm and y = 25 cm. 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 stress concentration 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 Figure 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 Figure 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 V-model 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 V-cells' B-spline surfaces are extracted, and inner surfaces, between consecutive V-cells, are deleted. The resulting B-rep model consists of 8064 B-spline surfaces. With a B-rep CAD tool (Rhinoceros R ), the shell is added as a B-rep volume and combined with the microstructure using the Boolean union operation. Subsequently, the microstructure on the outer side of the shell is trimmed away using the trimming operation with the shell volume's outer surface. Finally, the computational model is extracted with a Boolean intersection with the computational domain. Figure 18 depicts the selection of the computational domain and the final model with the respective surfaces for the boundary conditions.  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 Section 3.2.2 especially because homogenization in the linear case allows the creation of a look-up table. However, the discussed fully resolved model can be used to verify the homogenization. Homogenization 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-voidmaterial' [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. A remedy is offered by homogenization. As explained in Section 3.2.2, for a functionally graded microstructure it is sufficient to compute the effective material tensors C * T i 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.3.3, 3.3.4 and 3.3.5, the diameter of the rod in x−direction and rotation angle around the z−axis. In order to compute the respective microscopic material behaviors, homogenization simulations are carried out for unit tiles with three different configurations of the diameter of the rod in x−direction (Ø 0.2 mm, Ø 0.3 mm, and Ø 0.4 mm), yielding the unrotated, effective material tensors C * T i . For the homogenization simulations, the material of the microstructure is considered to be steel with a Young's modulus of E = 210 GP a, 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 the Equations (19), (20) and (21), respectively. One can identify different material behaviors, which is expected due to the respective unit tiles' geometrical features. The orientation and the thickness of the rods have 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 T i 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 the Appendix, Section 6 presents the respective independent material tensor entries C ii of the three unit tiles for arbitrary rotations around the z−axis, following Equation (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 Figure 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. In contrast, 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.2 mm) towards the interface of the shell (Ø = 0.4 mm). A uni-axial compression state is achieved by applying a uniform displacement of ∆z = −1.0 on the top surface and restricting the displacements in the z−direction on the bottom surface. Three additional point-bearings block the rigid body motions.   The simulation uses 15 × 15 × 15 high-order Legendre finite cells with a polynomial degree of p = 4. For the integration, moment-fitting with the depth of an underlying octree of d = 4 is chosen. At the interface between shell and infill, one h−refinement step is carried out to capture the material discontinuity. As the unit tiles' homogenization 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 Equation (26).   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. In order to illustrate such a topological change, a hole is drilled through the structure (see Figure 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 hole walls.

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 tri-variate spline-formulation, the inversion -that is necessary for the point inclusion test -turns 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 V-cells' control points was increased to carry material information, as well. During the integration -apart from the pointinclusion 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 microstructures' parametric description, 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 microstructure's complexity and the amount of varying geometrical features, the number of representative unit tiles might become large. As for each of these unit tiles, an individual homogenization simulation needs to be carried out. Thus, these 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 microstructure results, 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. Considering that such kind of boundary layer effects usually vanish rapidly away from the interface, the error is not dominant.

Availability of data and materials
The geometric models simulated and analysed during the current study are either reproducible with the provided information, or available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.

Funding
We gratefully acknowledge the support of the German Research Foundation (DFG) under Grant No. Ra 624/22-2. We would also like to thank the German Research Foundation (DFG) for its support through the TUM International Graduate School of Science and Engineering (IGSSE), GSC 81. Furthermore, we gratefully acknowledge the support of the Transregional Collaborative Research Centre SFB/TRR 277 Additive Manufacturing in Construction. The Challenge of Large Scale, funded by the German Research Foundation (DFG).

Authors' contributions
BW was the corresponding author who wrote the central part of the paper, integrated the Irit geometry kernel into the Adhoc++ FCM framework, and carried out most of the simulations. NK wrote the section on homogenization and carried out the simulation and classification of the unit-tiles. SK was responsible for the content regarding the finite cell method and implementation issues. ER guided the paper's general structure and contents, cross-checked the results, and proposed most examples. GE provided the geometry kernel and assistance for its access. Additionally, he was responsible for the content of the V-reps. All authors read and approved the final manuscript.

Acknowledgements
We acknowledge the contributions of the research groups at the chair of Computation in Engineering regarding the development of the finite cell method framework Adhoc++ and at the Center for Graphics and Geometric Computing concerning the development of the geometry kernel Irit and CAD software GuIrit.

Appendix
Effective material tensors of unit tiles This section is an extension of Example 3.3.4 and provides material data, which is used in Example 3.3.5. The following polar diagrams depict the independent entries C ii of the three material tensors (of the unit tiles) for an arbitrary rotation around the z−axis. The values are computed with the Bond transformation matrices [86], according to Equation (23). Thus, at an angle of 0 • the value equals the corresponding entry of the respective unrotated material tensor C * T i of Example 3.3.4. Additionally, for a rotational degree of 45 • the results are numerically verified (see Figure 26). 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 concerning 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 Figure 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. Figure 29: Independent elastic constants for tile 3 under rotation around the z−axis.