XIGA: An eXtended IsoGeometric analysis approach for multi-material problems

Abstract

Multi-material problems often exhibit complex geometries along with physical responses presenting large spatial gradients or discontinuities. In these cases, providing high-quality body-fitted finite element analysis meshes and obtaining accurate solutions remain challenging. Immersed boundary techniques provide elegant solutions for such problems. Enrichment methods alleviate the need for generating conforming analysis grids by capturing discontinuities within mesh elements. Additionally, increased accuracy of physical responses and geometry description can be achieved with higher-order approximation bases. In particular, using B-splines has become popular with the development of IsoGeometric Analysis. In this work, an eXtended IsoGeometric Analysis (XIGA) approach is proposed for multi-material problems. The computational domain geometry is described implicitly by level set functions. A novel generalized Heaviside enrichment strategy is employed to accommodate an arbitrary number of materials without artificially stiffening the physical response. Higher-order B-spline functions are used for both geometry representation and analysis. Boundary and interface conditions are enforced weakly via Nitsche’s method, and a new face-oriented ghost stabilization methodology is used to mitigate numerical instabilities arising from small material integration subdomains. Two- and three-dimensional heat transfer and elasticity problems are solved to validate the approach. Numerical studies provide insight into the ability to handle multiple materials considering sharp-edged and curved interfaces, as well as the impact of higher-order bases and stabilization on the solution accuracy and conditioning.

1 Introduction

Multi-material problems play an important role for a wide range of applications in engineering, such as problems involving composite or functionally graded materials, multiple phase interactions, or contact between components. However, efficiently and accurately predicting the physical response described by partial differential equations of multi-material problems remains challenging. Such problems often exhibit complex geometries, intricate material arrangements, or small features along with physical responses presenting large spatial gradients or discontinuities. Therefore, providing an accurate resolution for both the geometry and the physics around material interfaces is crucial. When performing analysis with traditional finite element methods (FEM), this can only be achieved by generating highly-refined body-fitted approximation meshes, which is a tedious task and is known to represent a substantial part of the analysis time, see Bazilevs et al. [5]. This is especially the case when the number of materials increases and when boundaries and material interfaces change in time due to physical phenomena, e.g., time evolving interfaces, see Kamensky et al. [40], or numerical processes, e.g., topology optimization, see Noël et al. [52].

Over the past decades, immersed boundary techniques have gained in popularity, providing an elegant solution for the multi-material problems characterized above. These methods can accommodate complex domain boundaries and interfaces without the need to construct a conforming body-fitted mesh. The first immersed boundary method was formally introduced by Peskin [56]. Since this first occurrence, several immersed finite element approaches have been presented in the literature. Fictitious domain methods, also known as embedded domain methods, circumvent the need to generate conforming analysis meshes by embedding the computational domain in a larger one and applying specific integration techniques. Variants of the method are based on different approaches to impose boundary and interface conditions with penalty methods by Ramière et al. [57], with Lagrange multipliers by Glowinski et al. [29], Glowinski and Kuznetsov [28], or Burman and Hansbo [9], or with Nitsche’s method by Hansbo and Hansbo [31], Dolbow and Harari [20], Burman and Hansbo [10], and Burman et al. [12]. Conversely, enrichment methods alleviate the need for generating conforming approximation meshes by capturing a priori known discontinuous behaviors within the mesh elements. Originally developed to represent moving fronts and crack propagation, enrichment based approaches have been extended to tackle various types of interface problems with strong and weak discontinuities. Among these techniques, a few noticeable ones are the Partition of Unity Method (PUM) proposed by Babuška and Melenk [3], the Generalized Extended Finite Element Method (GFEM) as introduced in Strouboulis et al. [68, 69], the eXtended Finite Element Method (XFEM) as proposed by Moës et al. [49] and Belytschko and Black [6], and the Interface enriched Generalized Finite Element Method (IGFEM) introduced by Soghrati et al. [65].

In most of the aforementioned papers, the implementation of immersed boundary techniques relies on low order approximation functions, in particular linear Lagrange basis functions, to represent both the geometry and the physics. Such a choice of basis functions suffers from several shortcomings in terms of geometry resolution, and accuracy of physical responses. Using p-version FEM, see Babuška [2] for FEM based on hierarchical basis functions or Patera [55] and Karniadakis and Sherwin [41] for spectral FEM based on higher-order spectral basis functions, enables improved accuracy of physical responses per degree of freedom (DOF) and leads to higher convergence rates with mesh refinement, i.e., h-refinement. Additionally, using higher-order functions for the geometry representation improves accuracy in the presence of curved interfaces and boundaries.

Several papers in the literature use higher-order bases, such as Lagrange or spectral ones, in combination with immersed boundary techniques. Parvizian et al. [54] and Düster et al. [22] introduced the Finite Cell Method (FCM). Similar to other fictitious domain approaches, the method extends the analysis domain to embed the physical one but makes use of higher-order Ansatz functions to approximate the extended variables. Numerous contributions focused on enrichment methods to accurately represent the geometry of and the physics around curved cracks. Wells et al. [73] studied the propagation of displacement discontinuities in strain-softening media with second order Lagrange bases. Working on crack propagation, Stazi et al. [67] used second order Lagrange bases for the finite element approximations, while Zi and Belytschko [75] extended this approach to higher-order enrichment functions. Tackling both strong and weak discontinuities, Cheng and Fries [15] resolved curved boundaries and interfaces by generating integration subcells with one curved side and by applying corrections to the enrichment formulation. Focusing on material interfaces, Dréau et al. [21] exploited the XFEM with a corrected enrichment and represented the geometry on sub-meshes finer than the one used for the mechanical fields. A similar approach was proposed in Legrain et al. [45]. Haasemann et al. [30] proposed a numerical integration strategy based on NURBS surfaces for higher-order XFEM and weak discontinuities. Lehrenfeld [46] used parametric mappings of the integration cells to reduce the interface representation error and the associated integration error. More recently, Saxby and Hazel [59] proposed a higher-order modified XFEM based on corrected basis functions for weak discontinuity problems. Working with higher-order spectral basis functions, Legay et al. [44] and more recently Chin and Sukumar [16] proposed a spectral XFEM approach to tackle weak discontinuity problems with curved interfaces.

Along with the development of IsoGeometric Analysis (IGA), using B-splines or NURBS as basis functions has become an increasingly popular approach, see Hughes et al. [37] and Cottrell et al. [17]. In IGA, both the geometry of a structure and its physical behavior are described using splines. From a geometry point of view, using B-splines and NURBS facilitates compatibility with Computer Aided Design (CAD) software. From an analysis point of view, using smooth and higher-order bases, such as quadratic and cubic B-splines, leads to more accurate physical responses per DOF than traditional $C^0$ finite element approaches, see Hughes et al. [35]; Evans et al. [25]; Hughes et al. [36].

Over the years, several research works have aimed at combining the advantages of immersed boundary techniques and smooth higher-order basis functions, such as B-splines, NURBS, and other variants. Based on fictitious domain approaches, Schillinger et al. [61] proposed a B-spline version of FCM. Kamensky et al. [40] further extended the concept to tackle fluid-solid interaction problems and introduced the term immersogeometric analysis. Höllig et al. [33] introduced the web-method using weighted extended B-splines as basis functions to solve Dirichlet problems. Modeling weak discontinuities and in particular material interfaces, Sanches et al. [58] developed an immersed boundary technique based on B-spline bases and a modified basis to locally interpolate the Dirichlet boundary conditions. Focusing on enrichment methods, Jiang et al. [39] proposed a robust Nitsche’s method to tackle interface problems with the XFEM based on B-spline basis functions. They used a separate locally refined mesh to improve the geometry representation. Jia et al. [38] solved curved material interface problems with XFEM based on NURBS and used curved integration elements for increased accuracy. To further resolve the interface geometry, Chen et al. [13] implemented the XFEM with locally refined B-splines to allow for adaptive local refinement around the interfaces. Recently, Elfverson et al. [23, 24] proposed a so-called cutIGA approach and a symmetric Nitsche’s method for imposing boundary conditions, as well as a drop of basis functions for improved stability.

To date, the scope and the applications of immersed boundary techniques with higher-order basis functions are rather limited in terms of the number of phases or materials, but also in terms of the geometric complexity considered. Tools to handle the geometric representation were investigated in numerous publications. Tran et al. [71] used several Level Set Functions (LSF) to accurately represent complex microstructures with multiple spatially close inclusions and avoid numerical artefacts using the XFEM. Moumnassi et al. [50] also used several LSFs to accurately represent sharp features and curved interface without mesh refinement. Xia et al. [74] proposed a matched interface and boundary method to tackle multi-material and triple junctions. Hou et al. [34] built specific approximations for interface elements presenting multiple material and triple-junctions. A similar approach has recently been proposed by Chen et al. [14], who treated triple-junction points with an immersed boundary technique through the construction of specific functions on interface elements. Soghrati [64] extended the IGFEM to handle multi-material interfaces by constructing special enrichment functions. However, most proposed frameworks lack versatility and do not offer a systematic approach to tackle multi-material problems in two and three dimensions.

In this paper, we propose a versatile XIGA approach to tackle multi-material problems in two and three dimensions. The geometry of the computational domain is represented implicitly by one or multiple LSFs. The LSFs are used to determine subregions of the computational domain that are associated to different phases and materials. This approach allows for handling of straight-edged and curved interfaces, as well as N-material junctions, in a systematic way. The governing equations are integrated separately on each material subdomain, and elements where multiple materials coexist are decomposed into single material integration subdomains. The finite element approximations for both the geometry and the mechanical fields use multi-variate B-splines. They are smooth higher-order basis functions and provide higher accuracy per DOF and higher convergence rates than traditional $C^0$ finite element bases. In this paper, a novel generalized Heaviside enrichment strategy is used with multiple enrichment levels to capture discontinuities at external boundaries and material interfaces. Boundary and interface conditions are weakly enforced using Nitsche’s method. Numerical instabilities associated with small material integration subdomains are mitigated by an adapted version of the face-oriented ghost stabilization. The combination of these ingredients results in a versatile and robust approach to tackle multi-material problems.

In most enrichment approaches, the approximation space is extended using different enrichments for each material domain. However, this approach can lead to an artificially stiffened physical response when one or more material domains are disconnected. To alleviate this issue, the approximation space was extended using different enrichments for each connected material subdomain in Terada et al. [70], and Hansbo and Hansbo [32]. Nonetheless, this approach can still lead to an artificially stiffened physical response when the intersection of a connected material subdomain with the support of a particular basis background function is disconnected, see Makhija and Maute [48]. This situation frequently arises when B-spline basis functions are employed rather than classical finite element basis functions. This is because B-spline basis functions have larger support regions than classical finite element basis functions. To tackle this issue, we build on previous work on Lagrange basis functions (see Makhija and Maute [48]) and enrich each individual basis function separately based on the topology of the material layout within the basis function support.

For multi-material problems, ghost stabilization procedures typically rely on the polynomial extension of function values within a material region of an element adjacent to a ghost facet to the ghost facet itself, see Burman and Hansbo [11]. However, with the enrichment strategy proposed in this work, such an extension is not well-defined when a material region within the element is disconnected. In particular, the polynomial extensions associated with different connected material subregions may differ. To overcome this issue, we introduce a ghost stabilization strategy that explicitly accounts for the topology of the material layout of elements adjacent to ghost facets.

The remainder of the paper is organized as follows. Sect. 2 focuses on the use of one or multiple LSFs to represent the geometry of external boundaries and material interfaces. The proposed XIGA formulation is detailed in Sect. 3. First, a brief discussion of B-splines as basis functions for finite element analysis is provided. Then, the immersed boundary technique, i.e., here the XFEM, is detailed in terms of enrichment strategy, creation of the model, formulation of the governing equations, stabilization via an adapted face-oriented ghost stabilization, and integration. Sect. 4 illustrates the capabilities and the robustness of the proposed XIGA approach by solving canonical two- and three-dimensional problems focusing on heat conduction and elasticity. Finally, Sect. 5 draws conclusions about the developed XIGA approach and proposes directions for future work.

2 Geometry representation

Although the proposed XIGA approach is not restricted to any particular geometry description method, the geometry of a computational domain is represented by one or multiple LSFs in this paper. This specific geometry representation was chosen to ease future work on evolving interfaces, in particular for design through topology optimization.

The Level Set Method (LSM) was developed by Osher and Sethian [53] to efficiently track front propagation. The method allows for the implicit representation of a geometry by a LSF. An iso-level $\phi _t$ of the LSF, generally chosen equal to 0, describes the interface $\varGamma _{\pm }$ between two regions $\varOmega _+$ and $\varOmega _-$ of the analysis domain $\varOmega$ via:

$$\begin{aligned} \begin{array}{ll} \phi ({\mathbf {x}}) < \phi _t, &{} \forall \, {\mathbf {x}} \in \varOmega _+,\\ \phi ({\mathbf {x}}) > \phi _t, &{} \forall \, {\mathbf {x}} \in \varOmega _-,\\ \phi ({\mathbf {x}}) = \phi _t, &{} \forall \, {\mathbf {x}} \in \varGamma _{\pm }.\\ \end{array} \end{aligned}$$

(1)

Fig. 1

Geometry description of a computational domain, made of three material domains $\varOmega ^{I}$, $\varOmega ^{II}$, and $\varOmega ^{III}$, using two LSFs, $\phi _1$ and $\phi _2$

An example of the type of multi-material problems addressed in this work is given in Fig. 1, where a computational domain $\varOmega$, made of three material domains $\varOmega ^{I}$, $\varOmega ^{II}$, and $\varOmega ^{III}$, is described using two LSFs, $\phi _1$ and $\phi _2$. In this paper, multiple LSFs $\phi _i({\mathbf {x}})$, $i=1, \dots , n$, are used to describe the external boundaries of and the material interfaces within the computational domain $\varOmega$. A multi-phase level set model is exploited, as introduced by Vese and Chan [72]. With n LSFs, a maximum of $2^n$ phases can be represented. In this paper, a phase represents a subregion of the analysis domain associated with a unique combination of positive or negative valued LSFs. A phase index ${\mathcal {P}}$ is assigned to each subregion based on the LSF signs. The phase assignment procedure is sequenced as follows. First, characteristic functions $f_i({\mathbf {x}})$, $i=1, \dots , n$, are used to characterize the point ${\mathbf {x}}$ with respect to an iso-level $\phi _t$ of the LSF $\phi _i({\mathbf {x}})$, i.e., whether ${\mathbf {x}}$ is inside, outside, or on the iso-level contour, as:

$$\begin{aligned} f_i({\mathbf {x}}) = \left\{ \begin{array}{ll} 0, &{} \phi _i({\mathbf {x}}) < \phi _t\\ 1, &{} \phi _i({\mathbf {x}}) > \phi _t\\ \text{ on } \text{ interface, } &{} \phi _i({\mathbf {x}}) = \phi _t\\ \end{array} \right. \end{aligned}$$

(2)

These characteristic functions $f_i({\mathbf {x}})$ are used to assign a unique index ${\mathcal {P}}({\mathbf {x}})$ to the point ${\mathbf {x}}$:

$$\begin{aligned} {\mathcal {P}}({\mathbf {x}}) = \sum _{j=1}^{n} 2^{j-1}\ f_{j}({\mathbf {x}}). \end{aligned}$$

(3)

Finally, a material describing the constitutive behavior is assigned to each phase. The phase indices ${\mathcal {P}}({\mathbf {x}})$ are associated with the corresponding material indices ${\mathcal {M}}({\mathbf {x}})$ through a map m following:

$$\begin{aligned} {\mathcal {M}}({\mathbf {x}}) = m({\mathcal {P}}({\mathbf {x}})). \end{aligned}$$

(4)

The phase and material assignment procedure is illustrated with the three-material problem in Fig. 2. First, the characteristic functions $f_i({\mathbf {x}})$ are evaluated based on the LSFs signs. In this picture, the minus sign indicates that $\phi _i({\mathbf {x}}) < \phi _t$ and thus, $f_i({\mathbf {x}}) = 0$. Then, the phase indices ${\mathcal {P}}({\mathbf {x}})$ are computed by Eq.(3) based on the characteristic functions $f_i({\mathbf {x}})$. Finally, a material is assigned to each phase, by associating a material index ${\mathcal {M}}$ to each phase index.

Fig. 2

Phase and material assignment procedure for a three-material problem described by two LSFs, $\phi _1$ and $\phi _2$

Each LSF $\phi _i({\mathbf {x}})$ is discretized on a mesh using basis functions $B_k({\mathbf {x}})$ as:

$$\begin{aligned} \phi _i^{h}({\mathbf {x}}) = \sum _{k} B_k({\mathbf {x}})\, \phi _{i}^{k}, \end{aligned}$$

(5)

where $\phi _{i}^{k}$ are the coefficients associated to the LSF $\phi _i({\mathbf {x}})$. In this work, B-spline basis functions are chosen. The coefficients and corresponding basis functions are used to evaluate nodal level set values on the background mesh generated for analysis, see Sect. 3. The LSF is interpolated linearly along the element edges to determine the intersection of the $\phi =\phi _t$ iso-contour with the element edges, see Subsect. 3.6. Using a linear interpolation inherently leads to a low order approximation of geometry, which might limit the accuracy of the physical response analysis. This issue can be alleviated by first refining the background mesh and then interpolating the LSF on this refined mesh.

3 XIGA formulation

This section focuses on the XIGA approach proposed in this paper. First, B-spline basis functions are briefly reviewed in Subsect. 3.1. Then, our novel enrichment strategy based on generalized Heaviside functions, used to accommodate multiple materials within a basis function support, is explained in Subsect. 3.2. Subsect. 3.3 summarizes the formulation of the governing equations considered in this work for heat conduction and elasticity. The techniques to enforce boundary and interface conditions and to stabilize the XIGA formulation are detailed in Subsects. 3.4 and 3.5 respectively. As several materials may coexist within the same element, special attention is required to perform numerical integration. This integration procedure is described in Subsect. 3.6.

3.1 B-splines for finite element analysis

In this paper, B-splines are used to approximate the level set and physics variable fields. This subsection briefly recalls basic concepts of constructing B-splines in one and multiple dimensions.

Considering a knot vector $\varXi = \{ \xi _{1},\xi _{2}, \dots ,\xi _{n+p+1} \}$, for which $\xi \in {\mathbb {R}}$ and $\xi _{1} \le \xi _{2} \le \dots \le \xi _{n+p+1}$, a univariate B-spline basis function $N_{i,p}(\xi )$ of degree p is constructed recursively starting from the piecewise constant basis function:

$$\begin{aligned} N_{i,0}(\xi ) = {\left\{ \begin{array}{ll} 1, &{} \text {if}\ \xi _{i} \le \xi \le \xi _{i+1},\\ 0, &{} \text {otherwise}. \end{array}\right. } \end{aligned}$$

(6)

The Cox de Boor recursion formula is used to obtain the basis functions for higher degrees $p>0$, see de Boor [18]:

$$\begin{aligned} N_{i,p}(\xi )= & {} \frac{\xi - \xi _{i}}{\xi _{i+p} - \xi _{i}}\ N_{i,p-1}(\xi )\nonumber \\&+ \frac{\xi _{i+p+1} - \xi }{\xi _{i+p+1} - \xi _{i+1}}\ N_{i+1,p-1}(\xi ). \end{aligned}$$

(7)

A knot is said to have a multiplicity k if it is repeated k times in the knot vector. The corresponding B-spline basis exhibits a $C^{p-k}$ continuity at that specific knot, while it is $C^{\infty }$ in between unique knots.

To tackle n-dimensional problems, multi-variate B-spline basis functions $B_{{\mathbf {i}},{\mathbf {p}}}(\varvec{\xi })$ are obtained by the tensor product of univariate B-spline basis functions. Denoting the parametric space dimension by $d_{p}$, a tensor-product B-spline basis is constructed starting from $d_{p}$ knot vectors $\varXi ^{m} = \{\xi _{1}^{m},\xi _{2}^{m},\dots ,\xi _{n_{m}+p_{m}+1}^{m}\}$ with $p_{m}$ being the polynomial degree and $n_{m}$ the number of basis functions in the parametric direction $m = 1, \dots , d_{p}$. A tensor-product B-spline basis function is generated from $d_{p}$ univariate B-splines $N_{i_{m},p_{m}}^{m}(\xi ^{m})$ in each parametric direction m using the formula:

$$\begin{aligned} B_{{\mathbf {i}},{\mathbf {p}}}(\varvec{\xi }) = \prod _{m=1}^{d_{p}} N_{i_{m},p_{m}}^{m}(\xi ^{m}), \end{aligned}$$

(8)

where the position in the tensor product structure is given by the index ${\mathbf {i}} =\{ i_{1}, \dots , i_{d_{p}} \}$, and the polynomial degree is denoted by ${\mathbf {p}} = \{ p_{1}, \dots , p_{d_{p}} \}$.

Fig. 3

Enrichment strategy for a basis function $B_k$ spanning the three material subdomains in a three-material problem

In this paper, Lagrange extraction, as introduced by Schillinger et al. [63], is used to facilitate a classical finite element implementation of the integration of the governing equations over the mesh of background elements, here defined as the tensor product of nonempty knot spans. This approach avoids the need to consider the non-elemental-locality of B-splines and simplifies the integration procedure for elements occupied by multiple materials.

3.2 Enrichment strategy

The XFEM was introduced by Moës et al. [49] and Belytschko and Black [6] to model crack propagation without remeshing. The method enables the prediction of discontinuous or singular behaviors within an element by adding specific enrichment functions to the classical finite element approximation. In this paper, we follow the work by Terada et al. [70], Hansbo and Hansbo [32], and Makhija and Maute [48], and use a generalized Heaviside enrichment strategy to introduce discontinuities along geometries and material interfaces. We further generalize the approach and enrich each basis function separately based on the material layout within the basis function support to ensure independent approximation on each connected material subregion.

The enrichment level selection procedure proposed in this paper is illustrated for a three-material problem in Fig. 3 for a basis function $B_k$ spanning the three material subdomains. The basis function support is delimited by a red dashed line. Within this support, $\text{ supp }(B_k)$, four separate connected material subregions $\varOmega _k^{\ell }$ exist, each occupied by one and only one material, such that $\text{ supp }(B_k) = \bigcup _{\ell =1}^{4} \varOmega _k^{\ell }$. Two subregions $\varOmega _k^{\ell =1}$ and $\varOmega _k^{\ell =2}$ are occupied by the same material I. Only one subregion $\varOmega _k^{\ell =3}$ is occupied by material II and one $\varOmega _k^{\ell =4}$ is occupied by material III. As four material subregions $\varOmega _k^{\ell }$ exist within the basis support, four enrichment levels $L_k=4$ are necessary.

Formally, considering a multi-material problem, a physics variable field ${\mathbf {u}}({\mathbf {x}})$ is approximated as:

$$\begin{aligned} {\mathbf {u}}^h({\mathbf {x}}) = \sum _{k=1}^{K} \sum _{\ell =1}^{L_k} \varphi _k^{\ell }({\mathbf {x}}) B_k({\mathbf {x}})\ u_{k}^{\ell }, \end{aligned}$$

(9)

where K is the number of background basis functions, and $L_k$ is the number of separate connected material subregions $\{ \varOmega _k^{\ell } \}_{\ell =1}^{L_k}$ in the support of background basis function $B_k$ for $k=1,\dots ,K$. The coefficient $u^{\ell }_{k}$ is the DOF associated with background basis function $B_k$ and material subregion $\varOmega _{k}^{\ell }$ for $k=1,\dots ,K$ and $\ell =1,\dots ,L_k$. The function $\varphi _k^{\ell }({\mathbf {x}})$ is an indicator function that determines whether a point ${\mathbf {x}}$ belongs to a material subregion $\varOmega _k^{\ell }$:

$$\begin{aligned} \varphi _k^{\ell }({\mathbf {x}}) = I_{\varOmega _k^{\ell }}({\mathbf {x}}) = \left\{ \begin{array}{ll} 1 &{} \displaystyle \text{ if }\ {\mathbf {x}} \in \varOmega _k^{\ell },\\ 0 &{} \displaystyle \text{ otherwise }. \end{array} \right. \end{aligned}$$

(10)

The set of functions

$${\mathcal {B}}:= \left\{ \varphi _k^{\ell } B_k : k \in \{ 1,\dots ,K\}\ \text{ and }\ \ell \in \{ 1,\dots ,L_k \} \right\} $$

possesses several useful properties. First, note that for each ${\mathbf {x}}\in \text{ supp }(N_k)$, there is one and only one $\ell$ for which $\varphi _k^{\ell }({\mathbf {x}})\ne 0$. Thus,

$$\begin{aligned} \sum \nolimits _{\ell } \varphi _k^{\ell }({\mathbf {x}}) B_k({\mathbf {x}}) = B_k({\mathbf {x}}), \end{aligned}$$

and as the background basis $\{ B_k \}_{k=1}^{K}$ forms a partition of unity, it follows that the functions in ${\mathcal {B}}$ do as well. The functions in ${\mathcal {B}}$ are also pointwise non-negative. Finally, as the background basis functions are locally linearly independent, so are the functions in ${\mathcal {B}}$. This indicates the functions in ${\mathcal {B}}$ form a basis. We refer to this basis as the enriched basis as it derives from the enrichment strategy.

3.3 Governing equations in discretized form

The proposed XIGA approach is not limited to any particular type of partial differential equation. However, in this paper, we restrict our attention to elliptic problems, namely multi-material linear elasticity and heat conduction problems. Either a linear elastic or linear diffusive material is assumed for each non-void domain.

In this work, the total residual ${\mathcal {R}}$, i.e., the discrete form of the governing equations, consists of four terms which are discussed subsequently. We solve for static equilibrium to enforce balance of linear momentum within each material domain $\varOmega ^I$, where I is the material index:

$$\begin{aligned} {\mathcal {R}}({\mathbf {u}},\delta {\mathbf {u}}) = {\mathcal {R}}_{Lin}^{{\mathbf {u}}} + {\mathcal {R}}_{D}^{{\mathbf {u}}} + {\mathcal {R}}_{Itf}^{{\mathbf {u}}} + {\mathcal {R}}_{Ghost}^{{\mathbf {u}}} = 0, \end{aligned}$$

(11)

where ${\mathbf {u}}$ and $\delta {\mathbf {u}}$ are the displacement field and the test function, respectively.

We solve for static equilibrium to enforce heat balance within each material domain $\varOmega ^I$, where I is the material index:

$$\begin{aligned} {\mathcal {R}}(\theta ,\delta \theta ) = {\mathcal {R}}_{Lin}^{\theta } + {\mathcal {R}}_{D}^{\theta } + {\mathcal {R}}_{Itf}^{\theta } + {\mathcal {R}}_{Ghost}^{\theta } = 0, \end{aligned}$$

(12)

where $\theta$ and $\delta \theta$ are the temperature field and the test function, respectively.

The first residual term ${\mathcal {R}}^{{\mathbf {u}}}_{Lin}$ for linear elasticity reads:

$$\begin{aligned} \begin{aligned} {\mathcal {R}}^{{\mathbf {u}}}_{Lin} = \sum _I \Bigg [&+ \int _{\varOmega ^I}\ \delta \varvec{\varepsilon } : \varvec{\sigma }\ d\varOmega \\&+ \int _{\varOmega ^I}\ \delta {\mathbf {u}} \cdot {\mathbf {b}}\ d\varOmega \\&- \int _{\varGamma _N^I} \delta {\mathbf {u}} \cdot {\mathbf {t}}_{N}\ d\varGamma \ \Bigg ], \end{aligned} \end{aligned}$$

(13)

where body loads, ${\mathbf {b}}$, are acting on the domain $\varOmega ^I$ and traction forces, ${\mathbf {t}}_N$, are applied on the Neumann boundary, $\varGamma _N^I$. The Cauchy stress tensor is denoted by $\varvec{\sigma } = {\mathbf {D}}\ \varvec{\varepsilon }$ and is obtained by multiplication of the infinitesimal strain tensor $\varvec{\varepsilon } = \frac{1}{2} \left( \nabla {\mathbf {u}} + \nabla {\mathbf {u}}^T \right)$ with the fourth order constitutive tensor ${\mathbf {D}}$, here for isotropic linear elasticity, expressed as a function of the Young’s modulus E and the Poisson ratio $\nu$ in Voigt notation in 2D as:

$$\begin{aligned} {\mathbf {D}} = \tilde{E}\left[ \begin{array}{cccccc} 1-\nu &{} \nu &{} 0\\ \nu &{} 1-\nu &{} 0\\ 0 &{} 0 &{} \frac{1-2\nu }{2}\\ \end{array} \right] , \end{aligned}$$

(14)

and in 3D as:

$$\begin{aligned} {\mathbf {D}} = \tilde{E}\left[ \begin{array}{cccccc} 1-\nu &{} \nu &{} \nu &{} 0 &{} 0 &{} 0\\ \nu &{} 1-\nu &{} \nu &{} 0 &{} 0 &{} 0\\ \nu &{} \nu &{} 1-\nu &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} \frac{1-2\nu }{2} &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} \frac{1-2\nu }{2} &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \frac{1-2\nu }{2}\\ \end{array} \right] , \end{aligned}$$

(15)

with $\displaystyle \tilde{E} = \frac{E}{(1+\nu )(1-2\nu )}$.

The first residual term ${\mathcal {R}}^{\theta }_{Lin}$ for heat conduction reads:

$$\begin{aligned} \begin{aligned} {\mathcal {R}}^{\theta }_{Lin} = \sum _I \Bigg [\&- \int _{\varOmega ^I}\ \delta \nabla \theta \cdot {\mathbf {q}}\ d\varOmega \\&+ \int _{\varOmega ^I}\ \delta \theta \ q_B\ d\varOmega \\&- \int _{\varGamma _N^I} \delta \theta \ q_{N}\ d\varGamma \ \Bigg ], \end{aligned} \end{aligned}$$

(16)

where body heat loads, $q_B$, are acting on the domain $\varOmega ^I$ and heat fluxes, $q_N$, are applied on the Neumann boundary, $\varGamma _N^I$. The heat flux, ${\mathbf {q}} = - \left( \varvec{\kappa } \cdot \nabla \theta \right)$, is obtained by multiplying the conductivity tensor $\varvec{\kappa } = \kappa {\mathbf {I}}$, here considering isotropic diffusion, by the temperature gradient $\nabla \theta$.

3.4 Weak enforcement of boundary and interface conditions

Boundary and interface conditions are imposed weakly via Nitsche’s formulation, see Nitsche [51]. To enforce prescribed displacements on Dirichlet boundaries, the static equilibrium in Eq. (13) is augmented with:

$$\begin{aligned} \begin{aligned} {\mathcal {R}}^{{\mathbf {u}}}_{D} = \sum _I \Bigg [\&- \int _{\varGamma ^{I}_{D}} \delta {\mathbf {u}} \cdot \left( \varvec{\sigma } \cdot {\mathbf {n}}_{\varGamma } \right) d\varGamma \\&+ \int _{\varGamma ^{I}_{D}} \delta \left( \varvec{\sigma } \cdot {\mathbf {n}}_{\varGamma } \right) \cdot \left( {\mathbf {u}} - {\mathbf {u}}_D \right) \ d\varGamma \\&+ \int _{\varGamma ^{I}_{D}} \gamma ^{{\mathbf {u}}}_{N}\ \delta {\mathbf {u}} \cdot \left( {\mathbf {u}} - {\mathbf {u}}_D \right) \ d\varGamma \ \Bigg ], \end{aligned} \end{aligned}$$

(17)

where a non-symmetric Nitsche formulation is considered, see for example Burman [8] and Schillinger et al. [62], and ${\mathbf {u}}_D$ is the displacement imposed on the Dirichlet boundary $\varGamma ^{I}_{D}$. The parameter $\gamma ^{{\mathbf {u}}}_{N}$ is chosen to achieve a desired accuracy in satisfying the boundary conditions and is a multiple of the ratio E/h, where E is the Young’s modulus of the considered material and h is the edge length of the intersected background element.

The same formulation is used to impose temperature on Dirichlet boundaries by augmenting the equilibrium in Eq. (16) with:

$$\begin{aligned} \begin{aligned} {\mathcal {R}}^{\theta }_{D} = \sum _I \Bigg [\&- \int _{\varGamma ^{I}_{D}} \delta \theta \ \left( {\mathbf {q}} \cdot {\mathbf {n}}_{\varGamma } \right) d\varGamma \\&+ \int _{\varGamma ^{I}_{D}} \delta \left( {\mathbf {q}} \cdot {\mathbf {n}}_{\varGamma } \right) \left( \theta - \theta _D \right) \ d\varGamma \\&+ \int _{\varGamma ^{I}_{D}} \gamma ^{\theta }_{N}\ \delta \theta \ \left( \theta - \theta _D \right) \ d\varGamma \ \Bigg ], \end{aligned} \end{aligned}$$

(18)

where a non-symmetric Nitsche formulation is again considered and $\theta _D$ is the imposed temperature on the Dirichlet boundary $\varGamma ^{I}_{D}$. The parameter $\gamma ^{\theta }_{N}$ is defined similarly to $\gamma ^{{\mathbf {u}}}_{N}$ as a multiple of the ratio $\kappa /h$, where $\kappa$ is the conductivity of the considered material.

Interface conditions between the displacements of materials I and J on $\varGamma ^{IJ} = \partial \varOmega ^{I} \cap \partial \varOmega ^{J}$ are imposed similarly using the non-symmetric formulation:

(19)

where the jump operator computes the difference in the considered quantity between material domains I and J as . The mean operator $\{ \bullet \}$ computes a weighted sum of the considered quantity over the materials I and J as $\{ \bullet \} = w^{I}\ \bullet ^{I} + w^{J}\ \bullet ^{J}$. The weights, following Dolbow and Harari [20] and Annavarapu et al. [1], are defined as:

$$\begin{aligned} w^{I} = \frac{\text{ meas }(\varOmega ^{I})/E^{I}}{\text{ meas }(\varOmega ^{I})/E^{I} + \text{ meas }(\varOmega ^{J})/E^{J}}, \end{aligned}$$

(20)

and:

$$\begin{aligned} w^{J} = \frac{\text{ meas }(\varOmega ^{J})/E^{J}}{\text{ meas }(\varOmega ^{I})/E^{I} + \text{ meas }(\varOmega ^{J})/E^{J}}, \end{aligned}$$

(21)

where $\text{ meas }(\varOmega ^{I})$ and $\text{ meas }(\varOmega ^{J})$ are the volume or the surface area of materials I and J within the intersected element in two and three dimensions respectively. The properties $E^{I}$ and $E^{J}$ are the Young’s moduli of materials I and J. The penalty parameter $\gamma ^{{\mathbf {u}}}_{Itf}$ is evaluated as:

$$\begin{aligned} \gamma ^{{\mathbf {u}}}_{Itf} = \frac{2\, E^{I}\, \text{ meas }(\varGamma ^{IJ})}{\text{ meas }(\varOmega ^{I})/E^{I} + \text{ meas }(\varOmega ^{J})/E^{J}}, \end{aligned}$$

(22)

where $\text{ meas }(\varGamma ^{IJ})$ is the surface area or the length of the interface within the intersected element in two or three dimensions, respectively.

Interface conditions on the temperature between material domains I and J are imposed as follows:

(23)

where the weights $w^{I}$, $w^{J}$ and the penalty parameter $\gamma ^{\theta }_{Itf}$ are evaluated similarly to the linear elastic case as described in Eqs. (20), (21), and (22), but by substituting the conductivity $\kappa$ for the Young’s modulus E of the considered material.

3.5 Face-oriented ghost stabilization

Numerical instabilities might arise when using enriched finite element techniques. This is the case when the contributions to the residual of one or several basis functions approximating the physics variable field vanish or when these contributions become linearly dependent. These issues typically occur when the level set field intersects elements such that small material subdomains emerge. This results in an ill-conditioning of the equation system and inaccurate prediction of solution field gradients along the interface.

Aside from mesh modification and update strategies, different techniques are available in the literature to mitigate this issue. Preconditioning schemes can be used to improve the condition number, e.g., geometric preconditioning (see Lang et al. [43] for the case of XFEM) and algebraic preconditioning (see de Prenter et al. [19] for the case of FCM). Other strategies focus on augmenting the weak form of the equations to promote a well-conditioned system, for example ghost stabilization (see Burman [7]), basis function removal (see Elfverson et al. [23] for the specific case of B-spline interpolation), and basis function aggregation (see Badia et al. [4]).

In this work, a generalized version of the face-oriented ghost stabilization, proposed by Burman and Hansbo [11], is used to fit the proposed enrichment strategy based on the basis functions supports. Despite an increase in the complexity of the implementation, this stabilization is selected as it is generally more robust than other strategies such as preconditioning, especially when working with higher-order approximation and nonlinear problems. Additionally, the ghost stabilization improves the prediction of heat flux and stress around the boundaries and interfaces when considering diffusion and elasticity problems, respectively. The proposed approach is described hereunder.

The domain $\varOmega$, made of the union of all material domains $\varOmega ^{I}$, is immersed in a background mesh. The set of all background elements in the mesh is denoted ${\mathcal {K}}$ and $K_{\varOmega }$ is the subset of background elements that have a non-empty intersection with $\varOmega$:

$$\begin{aligned} K_{\varOmega } :=\left\{ K \in {\mathcal {K}} : K \cap \varOmega \ne \emptyset \right\} , \end{aligned}$$

(24)

The set of interior facets of $K_{\varOmega }$ is denoted ${\mathcal {F}}_{int}$. Each interior facet $F \in {\mathcal {F}}_{int}$ is shared between two elements $\varOmega _{F}^{+}$ and $\varOmega _{F}^{-}$ of $K_{\varOmega }$. Finally ${\tilde{\varGamma }}$ is defined as the union of all material interfaces and geometric boundaries. The set of ghost facets is then taken to be:

$$\begin{aligned}&{\mathcal {F}}_{ghost} :=\Big \{ F \in {\mathcal {F}}_{int} : \varOmega _{F}^{+} \cap {\tilde{\varGamma }} \ne \emptyset ,\nonumber \\&\qquad \qquad \qquad \qquad \text{ or } \ \varOmega _{F}^{-} \cap {\tilde{\varGamma }} \ne \emptyset \Big \}. \end{aligned}$$

(25)

The set of interior facets ${\mathcal {F}}_{int}$ and the set of ghost facets ${\mathcal {F}}_{ghost}$ are illustrated on a three-material problem in Fig. 4. A three-material domain $\varOmega$ is immersed in a background mesh. The set $K_{\varOmega }$ of background elements with non-empty intersection with $\varOmega$ are hatched in grey. The facets that lie within $K_{\varOmega }$ form the set ${\mathcal {F}}_{int}$ and are drawn in red. Finally, the set of ghost facets is shown in yellow.

Fig. 4

Set of interior facets ${\mathcal {F}}_{int}$ and set of ghost facets ${\mathcal {F}}_{ghost}$ used for the face-oriented ghost stabilization for a three-material problem

Fig. 5

Description of ghost facet F shared between two adjacent background elements $\varOmega _{F}^{+}$ and $\varOmega _F^{-}$ for the face-oriented ghost stabilization for a three-material problem

Consider ghost facet $F \in {\mathcal {F}}_{ghost}$ shared between two adjacent background elements $\varOmega _{F}^{+}$ and $\varOmega _F^{-}$ as illustrated in Fig. 5. The outward facing normals to $\varOmega _{F}^{+}$ and $\varOmega _{F}^{-}$ along F are defined as ${\mathbf {n}}_{F}^{+}$ and ${\mathbf {n}}_{F}^{-}$, such that ${\mathbf {n}}_{F} = {\mathbf {n}}_{F}^{+} = - {\mathbf {n}}_{F}^{-}$. The material layout subdivides the element $\varOmega _{F}^{+}$ into $N_{F}^{+}$ connected subdomains $\{\varOmega _{F,i}^{+}\}_{i=1}^{N_{F}^{+}}$ and the element $\varOmega _{F}^{-}$ into $N_{F}^{-}$ connected subdomains $\{ \varOmega _{F}^{-} \}_{j=1}^{N_{F}^{-}}$ such that each subdomain is occupied by one and only one material.

The material index associated to $\varOmega _{F,i}^{+}$ for each $i=1,\dots ,N_{F}^{+}$ is denoted ${\mathcal {M}}_{F,i}^{+} = {\mathcal {M}}(\varOmega _{F,i}^{+})$ and the material index associated to $\varOmega _{F,j}^{-}$ for each $j=1,\dots ,N_{F}^{-}$ is ${\mathcal {M}}_{F,j}^{-} = {\mathcal {M}}(\varOmega _{F,j}^{-})$. Finally, the polynomial extension of the field ${\mathbf {u}}|_{\varOmega _{F,i}^{+}}$ to all of ${\mathbb {R}}^{d}$ for each $i=1,\dots ,N_F^{+}$ is defined as ${\mathbf {u}}_{F,i}^{+}$ and the polynomial extension of the field ${\mathbf {u}}|_{\varOmega _{F,j}^{-}}$ to all of ${\mathbb {R}}^{d}$ for each $j=1,\dots ,N_F^{-}$ as ${\mathbf {u}}_{F,j}^{-}$. The ghost stabilization ${\mathcal {G}}^{{\mathbf {u}}}_{F}$ for facet F, penalizing the jumps in the displacement gradients across the facet, is then taken to be:

(26)

where the set $J_{F,i}$ is defined so that:

$$\begin{aligned} J_{F,i}:= & {} \Big \{ j \in \{ 1,\dots , N_{F}^{-} \} : {\mathcal {M}}_{F,i}^{+}={\mathcal {M}}_{F,j}^{-}\ne 0,\nonumber \\&\qquad \qquad \text{ and }\ | \partial \varOmega _{F,i}^{+}\, \cap \, \partial \varOmega _{F,j}^{-} | \ne 0 \Big \}. \end{aligned}$$

(27)

The parameter $\tilde{k} = {2(k-1)+1}$ and is a jump operator such that:

(28)

and:

(29)

The operator $\partial _n^k(\bullet )$ is the $k^{th}$ order normal derivative operator and $\partial _n(\bullet ) = \nabla (\bullet ) \cdot {\mathbf {n}}_F$ with $\nabla (\bullet )$ the spatial derivative. The parameter p is the degree of the considered approximation. It should be noted that only the $p^{th}$ contribution is non-zero for maximally smooth splines, i.e., splines without repeating knots. The ghost penalty parameter $\gamma _{G}^{{\mathbf {u}}}$ is defined as a multiple of the Young’s modulus E of the considered material.

The formulation in Eq. (27) requires that:

$${\mathcal {M}}_{F,i}^{+}={\mathcal {M}}_{F,j}^{-}\ne 0,$$

assuming that the material index for void is zero, so that the ghost stabilization is only applied between ${\mathbf {u}}_{F,i}^{+}$ and ${\mathbf {u}}_{F,j}^{-}$ when $\varOmega _{F,i}^{+}$ and $\varOmega _{F,j}^{-}$ are occupied by the same non-void material. Additionally, the formulation requires that:

$$| \partial \varOmega _{F,i}^{+}\, \cap \, \partial \varOmega _{F,j}^{-} | \ne 0,$$

and the ghost stabilization is only applied between ${\mathbf {u}}_{F,i}^{+}$ and ${\mathbf {u}}_{F,j}^{-}$ when the boundaries of $\varOmega _{F,i}^{+}$ and $\varOmega _{F,j}^{-}$, $\partial \varOmega _{F,i}^{+}$ and $\partial \varOmega _{F,j}^{-}$, respectively, meet along a portion of the facet F with a non-zero measure, e.g., the boundaries meet along more than a point in two dimensions and along more than a line in three dimensions.

The subdivision of $\varOmega _{F}^{+}$ and $\varOmega _{F}^{-}$ into connected subdomains $\varOmega _{F,i}^{+}$ and $\varOmega _{F,j}^{-}$ and the associated material indices ${\mathcal {M}}_{F,i}^{+}$ and ${\mathcal {M}}_{F,j}^{-}$ is shown in Fig. 6 for two different material configurations. In the first case, marked by a red box, the background element $\varOmega _{F}^{-}$ is occupied by two connected material subdomains $\varOmega _{F,1}^{-}$ and $\varOmega _{F,2}^{-}$, while the background element $\varOmega _{F}^{+}$ is divided into three connected material subdomains $\varOmega _{F,1}^{+}$, $\varOmega _{F,2}^{+}$, and $\varOmega _{F,3}^{+}$. As $\varOmega _{F,1}^{-}$ and $\varOmega _{F,2}^{+}$ have the same material index ${\mathcal {M}}_{F,1}^{-} = {\mathcal {M}}_{F,2}^{+}$ and meet along the facet, the jump in the associated field gradients is penalized across the facet. The same holds for $\varOmega _{F,2}^{-}$ and $\varOmega _{F,1}^{+}$.

In the second case, marked by a yellow box, a different situation arises. While the background element $\varOmega _{F}^{-}$ is filled with only two different materials, the grey material lies within two connected subdomains within the element. The background element $\varOmega _{F}^{-}$ is thus divided into three connected material subdomains $\varOmega _{F,1}^{-}$, $\varOmega _{F,2}^{-}$, and $\varOmega _{F,3}^{-}$, while the background element $\varOmega _{F}^{+}$ is divided into two connected material subdomains $\varOmega _{F,1}^{+}$ and $\varOmega _{F,2}^{+}$. As $\varOmega _{F,1}^{-}$ and $\varOmega _{F,2}^{+}$ have the same material index, ${\mathcal {M}}_{F,1}^{-} = {\mathcal {M}}_{F,2}^{+}$, and meet along the facet, the jump in the associated field gradients is penalized across the facet. The same holds for the pairs $\varOmega _{F,3}^{-}$ and $\varOmega _{F,2}^{+}$, and $\varOmega _{F,2}^{-}$ and $\varOmega _{F,1}^{+}$.

Fig. 6

Subdivision of $\varOmega _{F}^{+}$ and $\varOmega _{F}^{+}$ into connected subdomains $\varOmega _{F,i}^{+},\ i=1, \dots , N_F^{+}$ and $\varOmega _{F,j}^{-},\ j=1, \dots , N_F^{-}$ for the face-oriented ghost stabilization for a three-material problem

The full contribution of the ghost stabilization to the residual equations is attained by summing over all ghost facets:

$$\begin{aligned} {\mathcal {R}}^{{\mathbf {u}}}_{Ghost} = \sum _{F \in {\mathcal {F}}_{ghost}} {\mathcal {G}}^{{\mathbf {u}}}_{F}({\mathbf {u}}, \delta {\mathbf {u}}). \end{aligned}$$

(30)

The ghost penalization for the temperature field is defined similarly as:

$$\begin{aligned} {\mathcal {R}}^{\theta }_{Ghost} = \sum _{F \in {\mathcal {F}}_{ghost}} {\mathcal {G}}^{\theta }_{F}(\theta , \delta \theta ), \end{aligned}$$

(31)

where the ghost stabilization ${\mathcal {G}}_{F}^{\theta }$ for facet F is:

(32)

with the ghost penalty parameter, $\gamma _{G}^{\theta }$, being defined as a multiple of the conductivity $\kappa$ for the considered material.

3.6 Numerical integration

Working with immersed boundary techniques and using the Heaviside enrichment, the weak form of the governing equations is integrated separately on each material subdomain. Elements occupied by more than one material are decomposed into conforming integration subdomains. In two (three) dimensions, a quadrangle (a hexahedron) is subdivided into a triangular (tetrahedral) integration mesh that conforms to the material interfaces. Gauss quadrature rules are used on the generated integration elements.

The subdivision strategy is illustrated with a two dimensional three-material problem in Fig. 7. First, to increase the accuracy of the interface detection, a primary subdivision is performed, and the background element is divided in four subtriangles. Then, a secondary subdivision is performed to construct a triangular mesh that conforms to the interfaces created by the LSFs. The LSFs are linearly interpolated along the element edges using ${\tilde{\phi }}_1$ and ${\tilde{\phi }}_2$. The intersections between the approximated LSFs and the element edges are determined, see black circles on the figure, and triangular integration elements are created. It should be noted that additional refinement of the background mesh can be carried out before constructing the integration mesh from the LSFs to achieve reduced geometric error, see Subsect. 4.4.

Fig. 7

Generation of conforming numerical integration mesh using a primary and a secondary subdivision for a three-material problem

4 Numerical examples

In this section, the versatility of the proposed XIGA approach and its ability to systematically and accurately address multi-material problems is demonstrated. First, the accuracy of the numerical predictions is investigated in Subsect. 4.1. The influence of the minimum size of material integration subdomains within a basis function and of the choice of the ghost penalty parameter on the conditioning of the system of equations is studied. Subsect. 4.2 focuses on the robustness of the approach with respect to the creation of small material integration subdomains by varying the location of the analysis domain within a fixed background mesh. In Subsect. 4.3, the approach is used to solve N-phase junction problems, and the convergence rates attained with h-refinement considering linear, quadratic, and cubic B-spline basis functions are investigated. The ability of the XIGA approach to handle non-planar interface configurations is assessed in Subsect. 4.4. Finally, a N-material problem is tackled in Subsect. 4.5.

In the following examples, the performances of the proposed XIGA approach are studied and characterized using three criteria: the system condition number, the relative $L^2$ error norm, and the relative $H^1$ error semi-norm. Ill-conditioning can affect the convergence of the solver of the system of equations. The condition number is used to assess the conditioning of the system of equations and is evaluated as:

$$\begin{aligned} \text{ cond }({\mathbf {A}}) = || {\mathbf {A}}^{-1} || \cdot || {\mathbf {A}} ||, \end{aligned}$$

(33)

where ${\mathbf {A}}$ is a matrix describing the system of equations to solve and $|| \bullet ||$ is the Frobenius norm.

For a generic physics variable field ${\mathbf {a}}$, the relative $L^2$ error norm is defined as:

$$\begin{aligned} L^2 = \sqrt{ \frac{\int _{\varOmega } \left( {\mathbf {a}}^h - {\mathbf {a}} \right) ^T\, \left( {\mathbf {a}}^h - {\mathbf {a}} \right) \ d \varOmega }{\int _{\varOmega } {\mathbf {a}}^T\, {\mathbf {a}}\ d \varOmega }}, \end{aligned}$$

(34)

where ${\mathbf {a}}$ is a reference solution, here chosen as either an analytical solution or a numerical solution computed on a sufficiently refined mesh, and ${\mathbf {a}}^h$ is the numerical solution evaluated with the XIGA approach.

For a generic physics variable field ${\mathbf {a}}$, the relative $H^1$ error semi-norm is evaluated as:

$$\begin{aligned} H^1 = \sqrt{ \frac{\int _{\varOmega } \left( \nabla {\mathbf {a}}^h - \nabla {\mathbf {a}} \right) ^T\, \left( \nabla {\mathbf {a}}^h - \nabla {\mathbf {a}} \right) \ d \varOmega }{\int _{\varOmega } \nabla {\mathbf {a}}^T\, \nabla {\mathbf {a}}\ d \varOmega }}, \end{aligned}$$

(35)

where $\nabla {\mathbf {a}}$ is the gradient of a reference solution, here chosen as either an analytical solution or a numerical solution computed on a sufficiently refined mesh, and $\nabla {\mathbf {a}}^h$ is the gradient of the numerical solution evaluated with the XIGA approach.

In all the following examples, the set of discretized governing equations is integrated using Gauss quadrature rules on each integration subelement depending on the order of the basis functions. In two dimensions, $2{\times }2$-, $3{\times }3$-, and $4{\times }4$-point integration rules are used for quadrangular integration elements and 7-, 12-, or 25-point integration rules are used for triangular integration elements for linear, quadratic, and cubic basis functions respectively. In three dimensions, $2{\times }2{\times }2$-, $3{\times }3{\times }3$-, and $4{\times }4{\times }4$-point integration rules are used for hexahedral integration elements and 11-, 35-, and 56-point integration rules are used for tetrahedral integration elements for linear, quadratic, and cubic basis functions respectively. The systems of discretized governing equations are built using an implementation of the XIGA approach within an in-house fully parallelized C++ code. They are solved by the direct solver PARDISO for 2D problems (see Kourounis et al. [42]), and by a GMRES algorithm for 3D problems, preconditioned by an algebraic multi-grid solver (see Gee et al. [27]). The choice of an iterative solver for 3D problems is guided by the larger number of equations and unknowns, and the increase in memory requirement, not by the nature of the considered problems. In some of the numerical examples, the systems of equations are poorly conditioned and condition numbers exceeding $10^{25}$ are observed. Despite these large condition numbers, the linear solve converged for all presented results.

4.1 Stability study with respect to material subdomains size and ghost penalty parameter

In this subsection, the accuracy of the XIGA approach is demonstrated by showing that low errors with respect to the analytical solution can be achieved. The conditioning of the system of equations with respect to the size of the created material integration subdomains and the value of the ghost penalty parameter is studied.

For this purpose, a single material bar is considered in three dimensions, with dimensions $L = 3 + \delta \, \text{ m }$, $l = 1\, \text{ m }$, and a cross-section area $A = l^2 = 1\, \text{ m}^2$. We consider a linear elastic problem defined on this geometry. The setup and boundary conditions are illustrated in Fig. 8. Considering linear elasticity, the material Young’s modulus is set to $E = 10.0\ \text{ N/m}^2$ and the Poisson ration is set to $\nu = 0.0$ to avoid any three-dimensional effect. The left side of the bar is clamped, i.e., ${\mathbf {u}}_D = [u_{Dx}\ u_{Dy}\ u_{Dz}]^T = [0.0\ 0.0\ 0.0]^T\, \text{ m }$, and the Nitsche’s penalty parameter is set to $\gamma _N = 100.0$. Three loading scenarios are considered: (i) a traction at the bar tip ${\mathbf {t}}_N = [t_{Nx}\ t_{Ny}\ t_{Nz}]^T = [5.0\ 0.0\ 0.0]\, \text{ N/m}^2$, (ii) a constant body load ${\mathbf {b}} = b_0 [1.0\ 0.0\ 0.0]\, \text{ N/m}^3$ with $b_0 = 2.0$ within the material, and (iii) a linear body load ${\mathbf {b}} = b_0 [x\ 0.0\ 0.0]\, \text{ N/m}^3$ with $b_0 = 2.0$ within the material. These three loading scenarios result in a one-dimensional linear elasticity problem and lead respectively to a linear, a quadratic, or a cubic displacement solution in x, that is the position in the horizontal direction:

$$\begin{aligned}&u_x(x) = \left\{ \begin{array}{lll} \displaystyle u_{Dx} + \frac{t_{Nx}}{EA} x,\\ \text{ for }\ t_{Nx} \ne 0.0,\, b_x = 0.0,\\ \displaystyle u_{Dx} + \frac{b_0}{2 EA} \left( 2Lx - x^2 \right) ,\\ \text{ for }\ t_{Nx} = 0.0,\, b_x = b_0, \\ \displaystyle u_{Dx} + \frac{b_0}{6 EA} \left( 3L^2x - x^3 \right) ,\\ \text{ for }\ t_{Nx} = 0.0,\, b_x = b_0 x, \end{array} \right. \nonumber \\&u_y = \quad 0,\nonumber \\&u_z = \quad 0. \end{aligned}$$

(36)

Fig. 8

Single material linear elastic bar with varying sliver size and ghost penalty parameter: problem setup, boundary conditions and solution for quadratic displacement solution using quadratic B-splines

Different locations of the end of the bar are considered generating different sliver sizes $\delta$ over the last element. The following sliver sizes are investigated:

$$\begin{aligned} \begin{aligned} \delta =\,&\big [ 0.001\ 0.002\ 0.0035\ 0.005\ 0.007\ 0.01\ 0.015\ 0.025\ \\&0.04\ 0.06\ 0.08\ 0.1\ 0.15\ 0.25\ 0.4\ 0.6\ 0.8\ 0.9 \big ] \times h, \end{aligned} \end{aligned}$$

where h is the background element size and is set to $h = 1\, \text{ m }$. Additionally, the following values are considered for the ghost penalty parameter:

$$\gamma _G = \left[ \, 0.0\ 10^{-9}\ 10^{-5}\ 10^{-4}\ 10^{-3}\ 10^{-2}\ 10^{-1}\ 1.0\, \right] .$$

The displacement fields are interpolated using linear, quadratic, and cubic B-spline basis functions.

The relative $L^2$ error norm and the relative $H^1$ error semi-norm are shown in Fig. 9 for the linear solution case, in Fig. 10 for the quadratic solution case, and in Fig. 11 for the cubic solution case. In these figures, each column presents a different interpolation order for the B-spline basis functions: linear, quadratic, and cubic. The rows display the relative $L^2$ error norm and the relative $H^1$ error semi-norm. Each ghost penalty value $\gamma _G$ is associated with a colored curve. Focusing on Fig. 9, the $L^2$ error norm and the $H^1$ error semi-norm values show that high accuracy can be achieved regardless of the B-spline order. This is expected as the exact solution is linear in x. In this case, the finite element solution is insensitive to the ghost penalty parameter choice or to the sliver size. Two additional observations are worth to be noted. A slight degradation of the $L^2$ and $H^1$ error occurs when the ghost penalty parameter is chosen too large, i.e., $\gamma _G \ge 1e-2$. In this case, the ghost stabilization acts as a coarsening operator and using high penalty values leads to errors similar to a coarsening of the mesh. When the ghost stabilization is turned off, i.e., $\gamma _G = 0$, or for small values of the penalty parameter, i.e.,$\gamma _G \le 1e-5$, a slight effect of the sliver size is observed on the $L^2$ and $H^1$ errors, namely it increases as the sliver vanishes as expected. Low ghost penalty leads to poorly conditioned systems, as further observed with the condition number in Fig. 12, and in turn to lower accuracy of the solution.

Fig. 9

Single material linear elastic bar undergoing a constant load at its tip, leading to a linear displacement solution

Fig. 10

Single material linear elastic bar undergoing a constant body load, leading to a quadratic displacement solution

Fig. 11

Single material linear elastic bar undergoing a linear body load, leading to a cubic displacement solution

The results associated with the second and third loading cases leading to a quadratic and a cubic displacement solution are presented in Fig. 10 and 11, respectively. The results support the observations made for the linear solution case. It should be noted that for the second case, linear B-splines are not sufficient to accurately represent the exact solution that is quadratic in x. For the third case, both the linear and quadratic B-splines are not able to capture accurately the exact solution that is cubic in x. These results suggest that our XIGA approach is accurate and leads to low errors with respect to the analytical solution if the basis function order is sufficient to represent the analytical solution.

The condition number $\text{ cond }({\mathbf {A}})$ is the same for all loading scenarios. The condition numbers are displayed in Fig. 12 for different interface configurations and values of the ghost penalty parameter $\gamma _G$. Figure 12 shows an increase in the condition number when higher-order bases are used. This is due to the increased number of basis functions that are supported on intersected elements, see de Prenter et al. [19]. Using ghost stabilization, the condition number can be significantly improved by using a penalty parameter value $\gamma _G \ge 1e-5$. This effect is particularly visible when using higher-order B-spline basis functions, i.e., quadratic and cubic.

Fig. 12

Single material linear elastic bar undergoing a constant load at its tip, leading to a linear displacement solution along the bar

4.2 Robustness study with respect to the creation of arbitrary integration subelements

To study the robustness of the method and the effectiveness of the stabilization with respect to different intersection configurations, a straight bar is immersed and rotated in a fixed background mesh. The problem is solved for linear elasticity with an imposed body load that is quadratic in $x_0$, the distance along the bar. The problem setup and boundary conditions are illustrated in Fig. 13 with $L = 1.0\, \text{ m }$, $l = 0.5\, \text{ m }$, and a cross-section area $A = l^2 = 0.25\, \text{ m}^2$. The bar is made of a single linear elastic material with a Young’s modulus $E = 10.0\, \text{ N/m}^2$ and a Poisson ratio set to $\nu = 0.0$ to avoid three-dimensional effects. The bar is clamped at its left extremity, and ${\mathbf {u}}_D = [u_{Dx_0}\ u_{Dy_0}\ u_{Dz_0}]^T = [0.0\ 0.0 \ 0.0]\, \text{ m }$. For all simulations, the Nitsche’s penalty parameter is set to $\gamma _N = 100.0$. When applying the ghost stabilization, the penalty parameter is fixed to $\gamma _G = 0.001$. A quadratic body load ${\textbf {b}} = [b_{x_0}\ b_{y_0}\ b_{z_0}]^T = b_0\, [x_0^2\ 0.0\ 0.0]\, \text{ N/m}^3$ with $b_{x_0} = 2.0$ is applied along the bar. This loading case yields a quartic displacement solution over the bar:

$$\begin{aligned} u_{x_0}(x)= & {} \displaystyle u_{Dx_0} + \frac{b_0}{12EA} \left( 4L^3x_0 - x_0^4\right) ,\nonumber \\&\qquad \qquad \text{ with }\quad b_{x_0} = b_0 x_0^2,\nonumber \\ u_{y_0}= & {} 0,\nonumber \\ u_{z_0}= & {} 0. \end{aligned}$$

(37)

Fig. 13

An immersed single elastic material bar rotated in a fixed 3D background mesh: problem setup, boundary conditions and solution using quadratic B-spline with rotation angles $\phi _y = 20^o$, $\phi _z = 20^o$ and a mesh size $h=0.125\, \text{ m }$

The bar is rotated by an angle $\phi _y$ around the y-axis and $\phi _z$ around the z-axis. The following angles with respect to the orientation of the background mesh are considered $\phi _y = \phi _z = [\, 10^o\ 20^o\ 30^o\ 40^o\ 50^o\ 60^o\ 70^o\ 80^o\, ]$, as illustrated in Fig. 14. The background mesh size is successively refined, and the following mesh sizes are considered: $h = [0.5\ 0.25\ 0.125\ 0.0625]\, \text{ m }$. Each setup is solved with linear, quadratic, and cubic B-splines. The system condition number, as defined in Eq. (33), is monitored. The accuracy of the evaluated physical responses is compared against the analytical solution given in Eq. (37) using the relative $L^2$ error norm and the relative $H^1$ error semi-norm, as defined in Eqs. (34, 35).

Fig. 14

Immersed bar rotated with an angle $\phi _y$ around the y-axis and an angle $\phi _z$ around the z-axis in a fixed three dimensional background mesh

For the configurations defined above, the condition number, the $L^2$ error norm, and the $H^1$ error semi-norm are given in Fig. 15. The first graph shows the evolution of the mean condition number averaged over the eight rotation angles in terms of the background mesh size. The second and third graphs present the mean relative $L^2$ error norm and the mean relative $H^1$ error semi-norm over the eight rotation angles with respect to the background mesh size. For each performance measure, the solid and dashed lines correspond to the application or the absence of the ghost penalty stabilization respectively. The error bars illustrate the range of $L^2$ and $H^1$ values when using the ghost penalty. The range is defined by the minimum and maximum values of the relative $L^2$ and $H^1$ errors. The use of linear, quadratic, or cubic B-splines is indicated by circle, triangle, or square marks.

Fig. 15

Robustness study on an immersed elastic bar in a fixed 3D background mesh: mean performance and error over eight rotation angles

When no stabilization is used (dashed lines in Fig. 15), an increase of the condition number is observed when using higher-order basis functions, and when using finer meshes. Introducing ghost penalty stabilization (solid lines in Fig. 15), the condition number can be significantly improved, especially when using higher-order basis functions, and the effect of the mesh refinement is mitigated. A higher accuracy is achieved when using higher order B-splines, as expected since the exact solution is quartic along the bar. The relative $L^2$ error norm and the relative $H^1$ error semi-norm values are similar with and without ghost stabilization. It should be noted that for cubic B-splines, ghost stabilization leads to a drastic improvement of the condition number at the price of a slight increase of the $L^2$ and $H^1$ errors. The error bars show that, using ghost stabilization, the finite element predictions are insensitive to the rotation of the bar. This suggests that the proposed XIGA approach exhibits robustness with respect to intersection configurations. Finally, the graphs show that the optimal convergence rates with respect to mesh refinement are recovered for all B-spline orders p, i.e., $p+1$ in the $L^2$ error norm and p in the $H^1$ error semi-norm.

Fig. 16

Problem setup, boundary conditions, and configurations for a linear elastic single material bar with N-phase junction

Fig. 17

Accuracy study on a linear elastic single material bar with N-phase junction configurations in three dimensions

4.3 Accuracy study for N-phase problems

In this section, the ability of the XIGA approach to handle N-phase junctions is investigated. Again, we consider a linear elastic problem with an imposed body load. The problem setup and boundary conditions are illustrated in Fig. 16 with $L = 1.0\, \text{ m }$, $l = 0.5\, \text{ m }$, and a cross-section area $A = l^2 = 0.25\, \text{ m}^2$. The bar problem presented in Subsect. 4.2 is reused, but with the geometric configurations illustrated in Fig. 16. The bar is made of four different phases filled with the same material, so that the numerical solution can be easily compared to an analytical one. The Young’s modulus of the material is set to $E = 10.0\, \text{ N/m}^2$ and a Poisson ratio $\nu = 0.0$ is chosen to avoid three-dimensional effects. The bar is clamped at its left extremity, ${\mathbf {u}}_D = [u_{Dx}\ u_{Dy}\ u_{Dz}]^T = [0.0\ 0.0 \ 0.0]\, \text{ m }$. For all simulations, the Nitsche’s penalty parameter is set to $\gamma _N = 100.0$ and the ghost penalty parameter is set to $\gamma _G = 0.001$. A quadratic body load ${\textbf {b}} = b_0\, [x^2\ 0.0\ 0.0]\, \text{ N/m}^3$ with $b_0 = 2.0$ is imposed, leading to a quartic displacement solution over the bar, as described hereunder:

$$\begin{aligned} \begin{array}{lll} u_x(x) = \displaystyle u_{Dx} + \frac{b_0}{12EA} \left( 4L^3x - x^4\right) ,\\ \text{ with }\quad b_x = b_0 x^2,\\ u_y = 0,\\ u_z = 0. \end{array} \end{aligned}$$

(38)

As illustrated in Fig. 16, four different configurations are considered: (1) a two-phase junction with the interface aligned with the y-axis, (2) a three-phase junction with the interfaces aligned with the x- and y-axes, (3) a four-phase junction with the interfaces aligned the x- and y-axes, and (4) a four-phase junction with the interfaces rotated by 45 degrees against the x- and y-axes. An h-refinement study is performed to investigate the convergence rates of linear, quadratic, and cubic B-spline basis functions. Five mesh sizes are considered, $h = [ 0.5\ 0.25\ 0.125\ 0.0625\ 0.03125]\, \text{ m }$.

The problem is solved for the four geometric configurations depicted in Fig. 16, successively refining the mesh. The performance in terms of system conditioning and solution accuracy is provided in Fig. 17. The graphs show the variations of the condition number, the relative $L^2$ error norm, and the relative $H^1$ error semi-norm for different mesh sizes. For each performance measure, the use of linear, quadratic, or cubic B-splines is indicated by circle, triangle, or square marks.

For all configurations, similar performance is observed in terms of condition number, $L^2$, and $H^1$ errors. This demonstrates the ability of the approach to handle N-phase junctions accurately. As expected, an increase in the condition number is observed for higher-order B-splines. In terms of accuracy, optimal convergence rates are recovered for both the $L^2$ and the $H^1$ error norms for all basis function orders. It should be noted that a slight increase in the condition number and in the errors is observed as the complexity of the geometric configurations is increased. Additionally, an increase in the condition number is observed for low refinement, due to the creation of small material integration subdomains with low volume ratio with respect to the background elements.

4.4 Accuracy study for two-material problems with curved interface

In the previous subsection, the surfaces and interfaces were planar and could be represented exactly with our level set approach, that inherently leads to a low order approximation of geometry, see Subsect. 2. In this subsection, the ability of the proposed XIGA approach to accurately resolve non-planar surfaces and interfaces is studied. A heated inclusion embedded in an infinite matrix problem is considered in two and three dimensions. The problem setup is illustrated in Fig. 18, where the dimensions are set to $L = 2.0\, \text{ m }$ and $a = 0.5\, \text{ m }$. The embedded inclusion is made of a material I with a conductivity $\kappa ^I = 1.0\, \text{ W/mK }$ and is undergoing a constant heat body load $q_B = 1.0\, \text{ W/m}^2$ or $\text{ W/m}^3$. The infinite medium is made of a material II with a conductivity $\kappa ^{II} = 0.125\, \text{ W/mK }$. The surrounding medium is not heated, i.e, $q_B^{II} = 0.0\, \text{ W/m}^2$ or $\text{ W/m}^3$. For the 2D case, a cylinder embedded in an infinite medium is considered. In cylindrical coordinates, the strong form of the heat conduction equation is given as:

$$\begin{aligned} \frac{1}{r} \frac{d}{dr} \left( r \frac{d \theta (r)}{dr} \right) + \frac{q_B}{\kappa } = 0, \end{aligned}$$

(39)

where r is the radius computed from the center of the cylinder and has the general solution:

$$\begin{aligned} \theta (r) = -\frac{q_B r^2}{4 \kappa } + C_1 \ln (r) + C_2, \end{aligned}$$

(40)

where $C_1$ and $C_2$ are integration constants that can be determined by considering the boundary conditions. A prescribed temperature $\theta _D = 0.375\, \text{ K }$ and a zero temperature gradient are prescribed at the origin. At the material interface, continuity is enforced for the temperature field and the heat flux. Finally, the temperature field over the computational domain is given as:

$$\begin{aligned} \theta (r) = \left\{ \begin{array}{lll} \displaystyle \theta _D -\frac{q_B^I r^2}{4 \kappa ^I}, &{}\text{ for }&{} r \le a,\\ \displaystyle \theta _D -\frac{q_B^I a^2}{4 \kappa ^I} - \frac{q_B^I a^2}{2 \kappa ^{II}} \ln \left( \frac{r}{a}\right) , &{}\text{ for }&{} r > a. \end{array} \right. \end{aligned}$$

(41)

For the 3D heat conduction problem, a sphere embedded in an infinite medium is considered. In spherical coordinates, the strong form of the heat conduction equation is given as:

$$\begin{aligned} \frac{1}{r^2}\frac{d}{dr} \left( r^2 \frac{d \theta (r)}{dr} \right) + \frac{q_B}{\kappa } = 0, \end{aligned}$$

(42)

where r is the radial coordinate. The general solution of Eq. (44) is:

$$\begin{aligned} \theta (r) = -\frac{q_B r^2}{6 \kappa } + \frac{C_1}{r} + C_2, \end{aligned}$$

(43)

where $C_1$ and $C_2$ are integration constants that can be determined by considering the boundary conditions. A prescribed temperature $\theta _D = 0.375\, \text{ K }$ and a zero temperature gradient are prescribed at the origin. At the material interface, continuity is enforced for the temperature field and the heat flux. Finally, the temperature field over the computational domain is given as:

$$\begin{aligned} \theta (r) = \left\{ \begin{array}{lll} \displaystyle \theta _D - \frac{q_B^I r^2}{6 \kappa ^I}, &{}\text{ for }&{} r \le a,\\ \displaystyle \theta _D - \frac{q_B^I a^2}{6 \kappa ^I} - \frac{q_B^I a^3}{3 \kappa ^{II}} \Big ( \frac{1}{a} - \frac{1}{r} \Big ), &{}\text{ for }&{} r > a. \end{array} \right. \end{aligned}$$

(44)

To simulate an infinite host domain, the analytical solutions, presented in Eqs. (41, 44), are applied to the outer faces of the host domain as a weak Dirichlet boundary condition. The Nitsche’s penalty parameter is set to $\gamma _N = 100.0$ for both the boundary and the interface conditions; the ghost penalty parameter is set to $\gamma _G = 0.001$. Simulations are performed for different mesh sizes $h = [0.5\ 0.25\ 0.125\ 0.06125\ 0.03125]\, \text{ m }$. The accuracy and convergence of the physical responses with mesh refinement are evaluated using the relative $L^2$ error norm and the relative $H^1$ error semi-norm, as defined in Eqs. (34, 35).

The temperature solution obtained for the five different mesh sizes with quadratic B-splines and a fixed integration mesh of size $h_{ \text{ int }} = 0.03125\, \text{ m }$ is shown in Fig. 19. Note that the geometry and the temperature fields are refined independently. The geometry is linearly interpolated and its accuracy can be increased by refining the integration mesh. A geometrical error, characterizing the accuracy of the cylinder or sphere representation, is evaluated as follows:

$$\begin{aligned} e_{ \text{ geo }} = \frac{V^h - V}{V}, \end{aligned}$$

(45)

where V is the reference volume or area in 2D or 3D, here chosen as $V = \pi a^2$ and $V = 4 \pi a^3/3$ in 2D and 3D respectively, and $V^h$, the numerical volume or area evaluated with the XIGA approach.

The results obtained in two and three dimensions, are shown in Figs. 20 and 21. In these figures, the rows present the relative $L^2$ error norm and the relative $H^1$ error semi-norm with respect to the h-refinement. The columns correspond to different levels of refinement of the integration mesh, i.e., $h_{ \text{ int }} = 3.125e^{-2}$, $7.813e^{-3}$, and $1.953e^{-3}\, \text{ m }$, that yield decreasing geometrical error $e_{ \text{ geo }}$, as shown on top of each graph. The highest level of refinement leads to a problem with around 12,200 integration elements and around 1,300 DOFs in 2D and around 11,500,000 integration elements and around 41,000 DOFs in 3D. Therefore, the cubic case with the highest refinement of the integration mesh is omitted in 3D, as the computational resources to tackle such a problem were not available.

Fig. 18

Problem setup and boundary conditions for a heated inclusion in an infinite host medium in two and three dimensions

Fig. 19

Heated inclusion in an infinite medium, refinement of the temperature field with a fixed geometry field $h_{ \text{ int }} = 0.03125\, \text{ m }$, solution using quadratic B-splines in three dimensions

Fig. 20

Heated cylindrical inclusion in an infinite medium, refinement of the temperature field for three different fixed geometry refinements $h_{ \text{ int }} = 3.125e^{-2}$, $7.813e^{-3}$, and $1.953e^{-3}\, \text{ m }$

Fig. 21

Heated spherical inclusion in an infinite medium, refinement of the temperature field for three different fixed geometry refinements $h_{ \text{ int }} = 3.125e^{-2}$, $7.813e^{-3}$, and $1.953e^{-3}\, \text{ m }$

In 2D and in 3D, the relative $L^2$ and $H^1$ error plots show that the ability of the proposed approach to represent the curved two-material interface problem solutions increases as the B-spline mesh is refined. However, optimal convergence rates in $L^2$ and $H^1$ errors are not fully recovered for all B-splines orders and geometry representations. When the integration mesh is not sufficiently refined, using higher-order bases provides only a slight improvement over using a linear basis. In this case, the lack of a sufficient geometry resolution limits the accuracy of the computed solution and prevents the recovery of optimal convergence rates of the $L^2$ and $H^1$ errors with h-refinement. As the integration mesh is further refined, the geometry error on the curved interface of the circle or the sphere drops and this issue is mitigated. Using an integration mesh size of $h_{ \text{ int }} =1.953e^{-3}\, \text{ m }$, optimal convergence rates in the relative $L^2$ and $H^1$ error norm is fully recovered for the quadratic B-splines and partially recovered for the cubic ones.

4.5 Accuracy study for N-material problems

The ability to accurately tackle multi-material problems in two and three dimensions is addressed in this subsection. In two dimensions, a square inclusion embedded in a four-material host matrix is considered. In three dimensions, a cubic inclusion is embedded in an eight-material matrix. The setup of the problem is illustrated in Fig. 22, where the dimensions are set to $L = 2.0\, \text{ m }$ and $a = 0.5\, \text{ m }$. A fixed temperature $\theta = 0\, \, \text{ K }$ is imposed on the outer faces of the host domain as a weak Dirichlet boundary condition. The ghost penalty parameter is set to $\gamma _G = 0.001$ and the Nitsche’s penalty parameter to $\gamma _N = 50.0$ for both the boundary and the interface conditions. The entire domain, inclusion and matrix, is undergoing a varying heat body load $q_B = \sin (2\pi x) \sin (2\pi y)\, \text{ W/m}^2$ in 2D and $\sin (2\pi x) \sin (2\pi y) \sin (2\pi z)\, \text{ W/m}^3$ in 3D. Two different material configurations are considered: a single material one and a multi-material one. In the multi-material setting, singularities arise at the sharp corners of the inclusion. As this is not the case with a single material problem, the two configurations are considered to assess the effect of such singularities on the accuracy of the numerical solution. In the single material case, the inclusion and the four, in 2D, or eight, in 3D, matrix materials share the same conductivity $\kappa = 1.0\, \text{ W/mK }$. In the multi-material case, the embedded inclusion is made of a material I with a conductivity $\kappa ^I = 1.0\, \text{ W/mK }$ and the host medium is made of four, in 2D, or eight, in 3D, different materials with a conductivity $\kappa ^{i} = i \times 0.125\, \text{ W/mK }$.

Simulations are performed using linear, quadratic, and cubic B-spline basis functions. The accuracy of the physical responses is evaluated by performing a h-refinement with $h = [ 0.5\ 0.25\ 0.125\ 0.0625\ 0.03125]\, \text{ m }$ and comparing the obtained results in terms of the relative $L^2$ and $H^1$ error norms, as defined in Eqs. (34, 35), with the response evaluated with a higher resolution $h = 0.015625\,\text{ m }$.

The inclusion presents sharp corners which are within a single background element. Using the same basis function to approximate the temperature field around the corner might degrade the accuracy of the solution. To investigate this situation, two schemes for the phase and material assignments are considered for the two-dimensional case. Each region in the matrix associated with identical material properties can be divided into one or three connected subregions leading to a total number of five or thirteen materials, as illustrated in Fig. 23. As detailed in Subsect. 3.2, the two schemes lead to different enrichment. When considering five materials, a single enrichment level $\ell =1$ is necessary to approximate the temperature field around the inclusion corner; while for the thirteen-material case, three enrichment levels are used $\ell =1,\dots ,3$.

Fig. 22

Setup and boundary conditions for a heated multi-material medium in two and three dimensions

Fig. 23

Schemes for phase and material assignment for a heated single or multi-material medium

For the two-dimensional case, the simulations use linear, quadratic, and cubic B-splines basis functions. The relative $L^2$ error norm and the relative $H^1$ error semi-norm with respect to the h-refinement are presented in Fig. 24. The highest level of refinement leads to a problem with 30,500 integration elements and about 22,000 DOFs. In the figure, the first column corresponds to the single material case and the second column to the multi-material case. Additionally, the results are presented for the five-material and the thirteen-material settings in solid and dashed line respectively.

Fig. 24

Accuracy study on a heated single or multi-material medium in two dimensions

In 2D for the single material setting, optimal convergence rates in the relative $L^2$ and $H^1$ error norms are recovered for linear, quadratic, and cubic B-spline basis functions. For lower level of refinement, the solution obtained using cubic B-splines is less accurate than for quadratic B-splines. This lower accuracy results from the effective coarsening introduced by the ghost stabilization. For the multi-material setting, optimal convergence rates are only recovered for linear and quadratic B-splines basis functions. Using cubic basis functions provides a slight improvement in terms of accuracy, but not in terms of convergence rates with respect to the quadratic basis functions. This is due to singularities at the sharp corners of the inclusion. Introducing auxiliary material domains and enhancing the approximation around the inclusion sharp corner allows for the recovery of optimal convergence rates.

For the three-dimensional case, only the scheme with the auxiliary material domains around the corners and edges is considered. The resulting setup has a total number of 57 materials; each area in the host medium is made of seven material domains with same properties. The simulations are carried out for quadratic B-splines basis functions only. The temperature solution obtained for the multi-material case and for the five different mesh sizes with quadratic B-spline is illustrated in Fig. 25. It should be noted that for a coarse mesh size as shown in Fig. 25(a) with $h=0.5\, \text{ m }$, the nature of the solution is not representable by the spline space.

Fig. 25

Multi-material heated medium with refinement of the temperature field solution using quadratic B-splines in three dimensions

The relative $L^2$ error norm and the relative $H^1$ error semi-norm with respect to the h-refinement are presented in Fig. 26. The highest level of refinement leads to a problem with around 15,400,000 integration elements and around 3,000,000 DOFs. Similarly to the 2D results, optimal convergence rates in the $L^2$ and $H^1$ error norms are recovered for the singe material case and close to optimal ones for the multi-material case, as singularities are present at the sharp corners of the inclusion.

Fig. 26

Accuracy study on a heated single or multi-material medium in three dimensions

This study suggests that our XIGA approach can deal with multi-material problems and provides the flexibility to enhance the approximation near singularities.

4.6 Complex geometries with curved interfaces

Finally, the ability of the proposed XIGA approach to analyze engineering problems with complex geometries, including multiple curved interfaces, is studied in this subsection. The elastic response of a sample of a fiber-reinforced composite is investigated. The geometry of the composite sample is generated using TexGen, an open source software for modeling the geometry of textile structures, see Long and Brown [47]. The setup of the problem is illustrated in Fig. 27 and is generated from the polyester.tg3 example, provided with TexGen v3.12.2. The outer dimensions of the sample are set to $L = 4.4\, \text{ mm }$ and $l = 8.5\, \text{ mm }$. The fiber yarns have an elliptical cross-section with the semi-axes set to $a = 1.8\, \text{ mm }$ and $b = 0.4\, \text{ mm }$. The composite is made of an epoxy matrix characterized by an elastic modulus $E^m = 2.5{\times }10^{3}\, \text{ N/mm}^2$ and a Poisson ratio $\nu ^m = 0.3$, and is reinforced by woven yarns made of carbon fibers with an elastic modulus $E^f = 250{\times }10^{3}\, \text{ N/mm}^2$ and a Poisson ratio $\nu ^f = 0.3$. To demonstrate the analysis capabilities of the XIGA approach, the following simple boundary conditions are considered: the composite sample is clamped on its back side with ${\mathbf {u}} = {\mathbf {u}}_D = {\mathbf {0}}$ and a constant traction is applied on its front side with $\mathbf {t}_N = [2500\, 0\, 0]$ N/mm$^2$. For this example, the three-dimensional simulation is carried out with the Nitsche’s penalty parameter set to $\gamma _N = 10.0$ and the ghost penalty parameter set to $\gamma _G = 0.01$. The displacements are approximated using quadratic B-spline basis functions.

Fig. 27

Setup and boundary conditions for a carbon fiber-reinforced composite sample

Fig. 28

Geometry details and integration mesh for a carbon fiber-reinforced composite sample

Fig. 29

Von Mises stress in the matrix and the fibers for a carbon fiber-reinforced composite sample

As illustrated in Fig. 28, the geometry of the composite sample presents multiple curved interfaces. The entire composite sample consisting of six fiber yarns and a matrix is immersed into a background mesh following the approach presented in Sect. 2. The level set fields are constructed as signed distance fields for each fiber. Close-ups on specific areas of the sample display details of the three-dimensional woven arrangement of the fibers and of the wrapping of the fibers in the matrix. The sample is made of five different materials, including the matrix in grey and four groups of fibers in red, blue, cyan, and yellow respectively. The same properties are assigned to all the fibers. Similarly to Subsect. 4.5, the use of distinct materials with identical properties allows modeling each fiber individually and prevents artificially merging them when the fibers are in contact. The background mesh size is set to $h = 0.05\, \text{ mm }$, leading to a model with around 1,800,000 integration elements and a system of equations with around 862,000 DOFs.

The displacement solution obtained is provided in Fig. 27. Figure 29 shows the resulting Von Mises stress contours in the matrix and the fibers separately. It should be noted that different color scales are used for the matrix and the fibers for visualization purposes. In the epoxy matrix, higher stresses are observed under the load as expected. For the fibers, larger stress values occur when they are loaded in tension. It should also be noted that local stress concentrations at the fiber surfaces are observed when the fibers are in contact. Note that in this analysis prefect bonding is assumed.

This example shows that the XIGA approach is applicable to engineering problems with multiple materials and complex interface geometries.

5 Conclusions

In this paper, an XIGA approach is proposed to simulate multi-material problems. To achieve a high resolution of both the geometry and the physical response around the material interfaces, the proposed approach combines an immersed boundary technique, namely the XFEM, and the use of smooth and higher-order bases, here B-splines. Although not restricted to this approach, the geometry of the computational domains is defined using one or multiple LSFs in this work. The iso-levels of the LSF define subregions within the computational domain that are associated with materials. To provide additional flexibility, the computational domains are fully immersed in a background mesh. The physical responses of the systems are evaluated with the proposed XIGA approach using a novel generalized Heaviside enrichment strategy, where each basis function is enriched separately based on the material layout within the basis function support. Boundary and interface conditions are imposed weakly through Nitsche’s formulation. Instabilities related to the creation of small material integration subdomains are mitigated using a generalized formulation of the face-oriented ghost penalty stabilization methodology adapted to our enrichment strategy.

The performance and the versatility of the proposed XIGA approach are studied through numerical experiments. In particular, canonical heat conduction and linear elastic problems in two and three dimensions are considered. The stability, the accuracy, and the robustness of the evaluated solutions are measured with the relative $L^2$ error norm, the relative $H^1$ error semi-norm, and the condition number.

Numerical experiments show that accurate solutions with optimal convergence rate with h-refinement can be recovered. Additionally, the ghost penalty stabilization mitigates the effect of small material subdomains on the conditioning of the system, but also on the achieved accuracy of the solution. This remark holds regardless of the order of the B-spline basis functions, i.e., linear, quadratic, and cubic. Additionally, by studying the effect of the size of the created material integration subdomains and of the value of the ghost penalty parameter, a valid range for the ghost penalty parameter is suggested for the considered types of problems.

Numerical examples study the application of our XIGA approach to problems with planar and curved interfaces. A slight degradation of the condition number and the solution accuracy is observed when increasing the geometric complexity. Optimal convergence rates with mesh refinement are observed for the $L^2$ and $H^1$ errors using linear, quadratic, and cubic B-splines. It should be noted that, when dealing with curved interfaces, the resolution of the geometry limits the accuracy of the evaluated physics variable field. Individual refinement of the geometry and solution field can be used to mitigate this issue, see Dréau et al. [21] and Legrain et al. [45]. Alternatively one can use a curved integration mesh, see Cheng and Fries [15], Haasemann et al. [30], Lehrenfeld [46], and Stavrev et al. [66], or Green’s theorem based integration schemes, see Saye [60].

Moreover, numerical simulations demonstrate that N-phase and N-material problems can be resolved accurately. The approach allows for a large flexibility in handling and assigning phases and materials to different subregions of the computational domain, even when dealing with complex geometries exhibiting multiple curved interfaces. For both multi-phase and multi-material problems, optimal convergence rates with mesh refinement are recovered for the $L^2$ and $H^1$ errors using linear, quadratic, and cubic B-splines. It should be noted that singularities arise at interfaces between more than two materials, and suboptimal convergence rates are observed. To tackle such singularities, the approximation space could be extended with specific enrichment functions if the explicit form of the singularity is known a priori, see for example Hou et al. [34] and Chen et al. [14].

Several follow-up research developments are foreseen. In particular, the XIGA approach can be studied for other classes of problems, such as nonlinear and multi-physics problems. Furthermore, hierarchical B-splines, LR-splines, or T-splines can be used in place of B-splines to carry out local mesh refinement, see Schillinger et al. [61] or Garau and Vàzquez [26]. This procedure allows for further enhancement of the resolution of both the geometry and the physics around the interface and the accurate modelling of systems across several length scales.