Thispaperpresentsanewcomputationalmethodforsolvingtwo-phaseengineeringproblems. The novelty of this approach lies in its ability to perform the analysis over the exact two-phase geometry. The geometry is defined by a spline-based representation and constructed based on the level-set method. The level-set-based geometry is mapped into a spline-based geometry using an untrimming technique that is newly developed in this paper. The spline-based representation is defined over unstructured meshes enabling precise representations of geometrically complex domains. An analysis-suitable mechanical model that accurately replicates the geometry is con- structed and the governing equations are solved using the isogeometric analysis method. The interface between phases is explicitly represented by cubic B-spline curves, allowing for the eval- uation of the solution field and various physical measures over its exact geometry. Additionally, algorithms that allow for changing the continuity of the solution field between the phases while precisely preserving the discontinuous two-phase domain are developed. Consistent analytical sensitivity analysis is formulated for generating the geometrical and mechanical models to allow for solving engineering problems with moving boundaries following the gradient-based approach. The proposed approach is tested for solving Poisson’s equation, linear elasticity, topology op- timization, and two-phase flow problems. The results expose the capability of the proposed approach to solve two-phase engineering problems while accurately representing the interface. 1. Introduction Various engineering problems involve multiple material phases to describe their physical be- havior. Examples include elasticity problems with multiple material properties, contact between multiple objects, and two-phase flow problems. In some cases, a static interface exists between the material phases, such as in steady-state fluid-structure interaction and static linear elastic- ity. Conversely, the shape and topology of the interface can be dynamic variables as in topology Email address: emad.shakur@ricam.oeaw.ac.at (Emad Shakur) Preprint submitted to November 16, 2023 optimization and dynamic fluid mechanics. Multiphase engineering problems are typically sen- sitive to the quality of representing the interface between the different phases. An inaccurate representation of the interface can lead to undesired inaccuracies in the solution of the physical problem. While the geometrical models for representing the multiphase domains are typically precisely defined by using phase-field [1, 2] or level-set-based methods [3, 4], the models that are used to numerically solve the physical problem – the mechanical models – are constructed using various approximation techniques. The key difference between these techniques lies in dealing with the interface representation in the mechanical model. With the advancement of the finite element analysis approach, various methods have been developed aiming to increase the quality of representing the interface in the mechanical model. These methods could be distinguished as diffusive interface methods and sharp interface methods. In diffusive interface methods, the exact interface is smeared to constitute a transition region of non-zero width that separates the domains of the different phases. The transition region is typically smooth and considered an intermediate phase that averages the physical properties of the adjacent phases. In two-phase flow problems, the intermediate phase is occasionally interpreted as an artificial mixing between the two fluids [5], whereas, in other multiphase engineering problems, it is considered as an approximation of the actual physical behavior. Following this class of methods, interactions between two phases, such as surface tension in two-phase flow problems, are usually formulated over the smeared interface region and not along the actual interface geometry. Diffusive interface methods have been used in various multiphase engineering problems such as topology optimization [6, 7], two-phase flow problems [8] and crack-propagation [9]. Despite the popularity and computational efficiency of diffusive interface methods in solving various engineering problems, the interfaces according to these methods are not well defined, making it challenging to evaluate physical and geometrical measures in proximity to the interface regions. In contrast, sharp interface methods aim for a crisp separation between the different material phases in the mechanical model. Sharp interface methods include the major class of immersed boundary techniques. A well-known immersed boundary techniques are the extended finite element method (X-FEM) [10, 11] and the CutFEM methods [12, 13, 14, 15]. In these methods, elements along the crossing interface are decomposed into sub-domains representing the different material phases. The decomposition is performed by enriching the finite element shape functions inthevicinityoftheinterfaceusingprojectionfunctions. X-FEMandCutFEMhavebeenusedto solvevariousmultiphaseengineeringproblemsincludingtwo-phaseflow[16,17], contact[18], and multiphase linear elasticity problems [19, 20, 21]. In other immersed boundary techniques, the finite elements mesh can be constantly modified to approximate the exact geometrical interface, such as in arbitrary Lagrangian-Eulerian methods [22]. A main difference between the various methods for constructing mechanical models in multi- phase engineering problems is their ability to approximate the actual geometry of the interface. The quality of the approximation impacts the accuracy of evaluating physical measures near 2 interfaces, such as interactions between the phases, and consequently, it might influence the accuracy of the analysis results. In general, inconsistencies between the mechanical model and the exact geometry emanate from two main reasons. The first reason is that the interfaces in the geometrical models are typically implicitly defined, making them hard to replicate in the mechanical model. The second reason is the difficulty of the traditionally adopted finite element method to precisely describe free-formed curved interfaces. In principle, the primary objective of the current paper is to solve this problem by developing a methodology characterized by a unified representation of the geometrical and mechanical models to eliminate inconsistencies between them. Additionally, to develop a representation that has an explicit and precise de- scription of the interfaces between the phases to allow for more accurate evaluations of physical and geometrical measures in the proximity of interface regions. The need to eliminate inconsistencies between the geometrical and mechanical models is not new, while several approaches have been developed aiming for a unified representation for both models [23, 24] leading to the development of isogeometric analysis (IGA) [25, 26]. The basic idea of IGA is that it adopts a spline-based basis to represent the mechanical model aiming to match the spline-based representation that is commonly used to represent geometrical models in Computer-Aided Design systems (CAD). This approach enables reaching a unified and explicit representation for the geometrical and mechanical models, rendering it an attractive choice for addressing multiphase engineering problems that require accurate interface representation. One of the main challenges in solving two-phase engineering problems is the inclusion of the trimmed patches at the interface regions in the mechanical model. Numerous studies have been conducted addressing this challenge using IGA, aiming to increase the accuracy of the analysis results. Kimetal.[27]introducedamethodthatdealswithtrimmedpatchesinB-splinesurfaces. Their method involves replacing the trimmed patches with triangular patches, each of which incorporates curved edges that accurately replicate given B-spline trimming curves. Various methods have been developed with the aim of improving the accuracy of physical results near trimmed patches bydeveloping specific quadrature rules [28, 29, 30]. Wei et al. [31] introduces an approach that combines immersed boundary techniques with conformal discretization near the interface aiming for more accurate physical responses. Lieu and Lee [32] solved multi-material topology optimization problems constructed over a B-spline domain, employing k-refinement techniques to enhance mesh resolution near the interface. Hosseini et al. [33] adopted IGA to solve two-phase flow problems where the geometry of the fluids is governed by the Cahn–Hilliard phase-field approach. In Yan et al. [34] high-order nonuniform rational B-splines (NURBS) were used to model multiphase flow problems with geometry based on the level-set method. IGA has been integrated with enrichment techniques to the basis functions to make physical results near interfaces more accurate. This led to the development of the Extended Isogeometric Analysis (X-IGA). Noël et al. [35] introduced a new generalized Heaviside enrichment strategy that addresses numerical instabilities that characterize the X-IGA approach to solving problems overdomainsthataccommodateanarbitrarynumberofmaterials. Thispaper’sconcepthasbeen furtherenhancedinSchmidtetal.[36]byutilizingtherefinabilityofspline-basedrepresentations, 3 resulting in an X-IGA procedure that is based on high-order hierarchical B-spline domains. It is worth noting that high-order hierarchical B-spline representation was previously used by Noël et al. [37] for geometric design. Additionally, X-IGA has been successfully adopted to solve fracture propagation [38, 39] and topology optimization problems [40, 41]. Most of the conducted research that uses IGA for solving multiphase engineering problems adopts implicit representation for the geometrical design and structured meshes with NURBS to represent the mechanical model. This obstructs the achievement of a unified representation for both the geometrical and mechanical models. While NURBS offer high polynomial degrees and continuity compared to finite element methods, their structured configuration limits their ability to accurately represent geometrically complex domains. This requires using projection procedures or other approximation techniques in order to capture the physical behavior near the interface between the material phases. Other research addresses this issue by using local refinement techniques to increase mesh resolution near the interface. This indicates that the full potential of IGA for integration over exact domains of complex geometries has yet to be realized. To take a step toward achieving this goal, the approach in the current paper adopts an explicit spline-based representation of the geometry, enabling a precise replication of the mechanical model. Nowadays, the topic of performing IGA over exact geometrically complex domains is on active research with a focus on domains that are defined over unstructured meshes [42, 43, 44]. A recent paper by the author introduces a new approach that uses IGA on unstructured meshes to solve single-phase topology optimization problems [45]. The main advantage of this approach is the adoption of untrimming techniques to generate unstructured meshes that can represent complex geometries in a precise and explicit manner. Additionally, the approach uti- lizes the ability of IGA to perform integrations over the exact domain of the geometry. This approach was used to solve stress-constrained topology optimization problems in [46]. The work in [46] emphasizes the importance of precisely representing the boundaries of the mechan- ical model in engineering problems where the physical behavior near the boundaries is crucial. However, this approach is developed mainly to solve single-phase linear elasticity problems and topology optimization and it is incapable of representing multiphase domains. Inspired by the work in [45], in this paper, a new method for solving two-phase engineering problems is developed. The novelty of this method with respect to current research on solving two-phase engineering problems is that it allows for a precise and explicit spline-based represen- tation for both the geometrical and the mechanical models. The main idea of this method is that it considers the implicit level-set-based representation as an intermediate step rather than a final representation of the geometry. Instead, the final geometry is defined by a spline-based rep- resentation, constructed using a linear mapping procedure that is applied to the level-set-based representation. The resulting spline-based representation possesses several analysis advantages and desired geometrical features: 1. Crisp phase separation: The proposed method utilizes a new untrimming technique that allows for a crisp separation between the phases in both the geometrical and mechanical models. 4 2. Accurate interface representation: The interface between the phases is accurately rep- resented as a set of cubic B-spline curves in both the geometrical and the mechanical models. This precise representation enables the evaluation of various physical and geomet- rical measures along the exact interface, including curvature, tangent vector, and values of the solution field, while also ensuring compatibility with CAD systems. 3. Conforming mesh: This approach provides a conforming mesh between different phases, eliminating the requirement for special treatments to achieve a tight bond between phases. 4. Precise geometry replication by the mechanical model: The mechanical model according to this approach precisely replicates the physical domain of the geometry, eliminating the need for approximation procedures between the two models. 5. Controlled continuity: The approach provides control over the continuity between different phases in the solution field without affecting the interface’s geometry. This feature is particularly useful when solving multi-physics engineering problems where multiple partial differential equations (PDEs) with varying continuity requirements within the same two- phase domain are solved. 6. Analytical sensitivity analysis: All stages of this method are differentiable, allowing for the development of a consistent and analytical sensitivity analysis for the entire solving proce- dure. This facilitates the integration of the proposed approach into engineering problems involving moving interfaces. Additionally, it is worth noting that this approach relies on the level-set method for the generation of the two-phase domain, allowing it to efficiently handle topological changes that mayoccurinengineeringproblemswithmovinginterfaces. Tothebestoftheauthor’sknowledge, these advantages have not yet been integrated into a single procedure in the literature for solving two-phase engineering problems. The remainder of this paper is organized as follows: Section 2 presents the procedure for generatingthegeometricalandmechanicalmodels. Algorithmsforchangingthecontinuityofthe solution field between the phases are given in Section 3. The sensitivity analysis for generating the geometrical and mechanical models is developed in Section 4. The ability of the proposed approach to solve two-phase engineering problems is demonstrated in Section 5. A discussion and concluding remarks are given in Sections 6 and 7, respectively. Finally, Appendices A and B are attached to the paper to provide further insights into the proposed approach. 2. Generating the geometrical and mechanical models This section outlines the proposed procedure for generating the geometrical and mechanical models for solving two-phase engineering problems. The procedure assumes a set of given trim- ming curves, ζ, that are used to construct the interfaces in the geometrical model. The trimming curves in the current paper are generated following the level-set method, assuming a bi-cubic 5 B-spline surface in R3as the level-set function. The z-coordinates of the B-spline surface served as control parameters of the level-set function and are denoted by α. It should be noted that the curves in ζdo not directly define the interfaces in the geometrical model. Instead, they are mapped into cubic B-spline curves that define the interfaces. This mapping is performed using untrimming techniques, resulting in a spline-based model defined over an unstructured mesh. The procedure for generating the geometrical and mechanical models will be demonstrated using the example in Figures 1, 2, and 5. Figure 1 displays the level-set function that is used for the illustrative example. The zero contour line of the level-set function is considered as the trimming curves ζ. The level-set function is defined by a B-spline surface constructed over a mesh of 11×11control points and open uniform knots vectors. Figure 1: The level-set function that is used to demonstrate the procedure for generating the geometrical and mechanical models. The zero contour level of the level-set function is marked with red lines. 2.1. Generating the geometrical model To generate the geometrical model using the given trimming curves, ζ, the following steps are performed: Step1. Define a bi-cubic B-spline surface, S0, that covers the entire physical domain of the problem. S0is considered as an initial geometry for the untrimming procedure. The knots vector and the planar locations of the control points of S0are chosen to be equal to those that define the level-set function. Figure 2a shows the initial B-spline surface, while Figure 2b shows its control polygon with black and red lines representing unit and zero knot intervals, respectively. Step2. Draw the curves in ζon top of S0and its control polygon. Classify the trimmed faces in the control polygon according to the trimming manner, as shown in the first column of Table 1. The trimming manner for each patch is determined based on the values of the level-set function at its corner points. These trimming manners can be represented as cases 1, 2, and 3 in the table. Additionally, case 4 is introduced to handle zero-measure faces in the control polygons and by that to ensure proper Bézier end conditions for the trimmed patches along the boundary of S0. Note that each case in the table can be rotated by 90 degrees to cover all the possible trimming cases. While cases 1, 2, 6 (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 2: The steps involved in generating the geometrical model: (a) A bi-cubic B-spline surface S0is initially created. (b) The control polygon that corresponds to S0. (c) The curves in ζdrawn on top of S0. (d) The curves in ζdrawn on top of the control polygon in (b). (e) The unstructured control polygon after modifying the trimmed faces according to the trimming manner presented in Table 1. (f) The control polygon after computing the locations of the control points PG. (g) The physical domain corresponds to the control polygon in (f). (h) The control polygon after performing Laplacian smoothing. (i) The physical domain corresponds to the control polygon in (h). and 4 can be directly determined by the values of the level-set function at the corners of the Bézier patches, case 3 requires additional knowledge of the contour path inside the domain. For the same level-set function values at the corners associated with case 3, two possible different configurations exist: the first is the one presented in the third row of the second column in the table, and the second is its rotation by 90 degrees. The correct configuration is determined by analyzing the contour path within the domain using square marching algorithms. For further details on this topic, readers are referred to the works in [47, 48]. Figures 2c and 2d demonstrate the curves in ζdrawn on top ofS0and its control polygon, respectively. Step3. Replace all the trimmed faces according to their corresponding trimming manner as illustrated in the second column in Table 1. The blue dots in the same column represent the control points of the control polygon after the untrimming procedure. The filled dots 7 Trimming manner Replacement C a s e 1 ? > 0 ? > 0? < 0 ? > 0 → C a s e 2 ? > 0 ? > 0? < 0 ? < 0 → C a s e 3 ? < 0 ? > 0? < 0 ? > 0 → C a s e 4 ? > 0? < 0 → Table 1: Modifying the trimmed faces in the control polygon of S0according to the trimming manner. First column: trimmed faces with the possible contours configurations drawn in dashed lines. Second column: the corresponding modifications that should be applied to the trimmed faces. represent control points already existing in S0, whereas the hollow dots represent newly added control points. The positions of the new control points are computed according to the subsequent step. The thick lines indicate the boundary edges of the different phases. The newly added edges that intersect boundary edges are assigned with zero knot intervals while all the other edges are assigned with unit knot intervals. In the table, black and red lines mark unit and zero knot intervals, respectively. The control polygon after this step is shown in Figure 2e. Note that the mesh after this step becomes unstructured due to the existence of extraordinary points, which are interior vertices in the control polygon that are not T-junctions and the number of their connected edges 8 is not 4. Step4. Compute the control points of the untrimmed geometry with respect to the control parameters α. It is important to note that the computation is not unique and can be chosen freely by the user. However, it is recommended to choose a smooth and differ- entiable computation method to facilitate the evaluation of sensitivities when dealing with problems with moving interfaces. We denote the control points of the mesh after this step by PG. These control points are a function of the control parameters of ζ, i.e.PG(α). A suitable computation procedure for positioning the newly added control points is provided in Appendix A. This Appendix refers to the 4 cases presented in Table 1. The relation between PGandαin the provided computation is linear and can be represented by: PG=Aα, (1) whereAis the correlation matrix. This computation procedure served to compute PG for the illustrative example in this section and the examples in Section 5. The resulting control polygon after this step is displayed in Figure 2f, while the corresponding physical domain is shown in Figure 2g. Step5. Perform Laplacian smoothing to the control points in the domain. This step is highly beneficial for avoiding interface layers with very small-sized patches, as they could affect the conditionality of the linear system. Additionally, implementing this step leads to designs with reduced interface curvature. This is particularly advantageous in applica- tions sensitive to high curvature values along the interface, such as problems related to stress concentrations. While this step successfully addresses these issues in the examples in the current paper, further in-depth studies on this matter are encouraged. The modified location of a certain control point PG iafter one Laplacian smoothing iteration is obtained by: PGℓ i=PNlp j=1wjPG jPNlp j=1wj(2) where PGℓ iis the modified location of PG i;PG jrepresents a control point that is con- nected to PG i;Nlpis the number of control points that are connected to PG i; and wj is the Laplacian smoothing weight that is associated with control point PG j. All the weights are set to 1at this step. Corner control points that have two zero-measure edges are excluded from the Laplacian smoothing procedure and their locations are not modified. The resulting control polygon and its corresponding physical domain after a single Laplacian smoothing iteration are displayed in Figures 2h and 2i, respectively. The modified control points PGℓare linearly correlated with the control points before 9 Laplacian smoothing, PG. This correlation can be expressed as: PGℓ=LGPG(3) whereLGis the correlation matrix. The physical domains corresponding to the unstructured meshes in the present work are constructed using the Bézier extraction procedure described in [49]. We note that the proposed approach distinguishes itself from other level-set-based methods in that the actual geometry of the two-phase domain is defined by the geometrical model presented in Figure 2g (or Figure 2i if Step 5 is applied), rather than the level-set-based geometry whose boundaries correspond to the zero contour level of the level-set function. The level-set-based geometry following the presented approach serves only as a preliminary stage in constructing the spline-based geometry. Meanwhile, the primary objective of adopting the level-set method is to govern the topological changesthatmayoccurinthetwo-phasedomainwhensolvingengineeringproblemswithmoving interfaces. Thus, it is worth clarifying that the geometries depicted in Figures 2g and 2i are neither an approximation of the level-set-based topology nor a precise replication of it, rather, they are the outcome of the mapping procedure introduced in this section. The main aim of this mapping procedure is to generate geometrical models with topologies that follow the ones of the level-set models not to replicate or approximate them. It should be noted, that the suggested mapping procedure is differentiable with derivatives presented in Section 4. These derivatives are taken into account when solving the examples in Section 5. The continuity of the spline-based model is derived from the adopted Bézier extraction pro- cedure. According to the Bézier extraction procedure outlined in [49], the continuity between Bézier patches that belong to the same phase and are at least two patches away from an extraor- dinary point is C2, while the continuity between patches in the proximity of an extraordinary point can be either C0orC1. Increased continuity between patches in the proximity of extraor- dinary points can be achieved by applying smoothing algorithms, such as in [43]. Moreover, the continuity between the different phases is determined by the values of the knot intervals that separate them. According to the knot intervals in Table 1, the different phases are separated with two layers of zero knot intervals. These values result in a continuity of C0between the different phases, i.e., in the direction perpendicular to the interface. Nevertheless, the continuity between the phases can be adjusted to C−1,C1orC2according to the requirements of the physical problem as will be discussed in Section 3. An advantage of the proposed approach is the ability to represent interfaces between phases as smooth cubic B-spline curves, enabling precise evaluation of physical measures at exact points along the interfaces. Examples of these measures include stresses, pressure, and velocities, which may vary depending on the engineering problem that is solved. The procedure for extracting the interfaces from the unstructured control polygon is elaborated in Appendix B. The interfaces in the geometrical model are highlighted in green lines in Figures 2g and 2i. The geometrical model generated so far is not necessarily analysis suitable. As discussed in [49] and [43], extraordinary points should be at least two patches away from each other to 10 avoid linear dependencies between the basis functions. This condition is not guaranteed in the geometrical model following the proposed approach. For this purpose, a mechanical model that ensures the analysis suitability of the mesh is constructed. It should be emphasized that the shape and topology of the B-spline curves that represent the interfaces are precisely preserved in the mechanical model. 2.2. Generating the mechanical model An analysis suitable mechanical model is generated by applying refinement and Laplacian smoothing procedures to the geometrical model. These procedures result in a mechanical model that is constructed over an unstructured mesh which includes T-junctions. The basis func- tions corresponding to the T-junctions’ vertices are evaluated following the theory of T-splines [50], while Bézier extraction to the patches near T-junctions is performed following the Bézier extraction algorithm in [51]. To generate the mechanical model, apply the following steps: Step1. Perform h-refinement to certain faces in the control polygon to ensure the analysis suit- ability of the mesh. The refinement is performed primarily to avoid linear dependencies between the basis functions. The refinement is performed as follows: – All faces that touch an extraordinary point are subjected to 2 h-refinement levels. This ensures that the extraordinary points are separated from each other by at least 4 patches. We denote these faces by f1. – Allzero-measurefacesinthegeometricalmodelaresubjectedto2h-refinementlevels. This refinement is performed to ensure that all extraordinary points are at least 3 patches away from the interface. As a result, the reduced continuity lines in the vicinity of an extraordinary point are pushed away from the interface. We denote the refined faces in this step by f2. – All faces that share an edge with f1andf2are subjected to one h-refinement level. The aim of applying this refinement is to avoid significant differences in the sizes of neighboring patches in the physical domain. – The refinement steps so far lead to the existence of T-junctions in the domain. To ensure the analysis suitability of the mesh in the regions of T-junctions, additional edges and vertices are added to the mesh. These edges and vertices are added to avoid intersections between T-junctions extensions, as elaborated in [52]. To ensure the interfaces of the geometrical model are accurately replicated in the me- chanical model, the newly added control points along interface edges in the mesh are positioned using the B-spline curve refinement algorithm outlined in [26]. All other newly added control points are equally spaced between the vertices of the faces. The control points of the mesh after this step are denoted by PMand are linearly correlated withPGℓ. This correlation can be expressed by: PM=TPGℓ, (4) 11 where Tis the correlation matrix between PMandPGℓ. Figures3and4illustratetherefinementprocedureinthevicinityofextraordinarypoints and zero-measure faces, respectively. Furthermore, these figures depict the formation of T-junctions after the refinement procedure. The refined mesh after applying the refinement procedure to the mesh in Figure 2f is demonstrated in Figure 5a. A zoom- in to the interface region is demonstrated in Figure 5b. The face extensions of the T-junction according to [52] are represented by green lines in both figures. Step2. To further enhance the quality of the mesh, a Laplacian smoothing procedure is per- formed on all internal control points in the control polygon. The interface control points remain fixed in place to preserve the shape and topology of the interface from the pre- vious step. This step is again linear and can be expressed as: PMℓ=LMPM, (5) where PMℓrepresents the control points after Laplacian smoothing and LMis the correlation matrix. Figure 5c shows the control polygon after the Laplacian smoothing procedure. A zoom- in to the interface region is presented in Figure 5d. This refined mesh is considered to be the final control polygon that defines the mechanical model. The physical domain of the mechanical model is shown in Figure 5e, and a zoom-in to the interface region is shown in Figure 5f. (a) (b) Figure 3: Refinement in the vicinity of extraordinary points. (a) The control mesh before the refinement. f1 faces are marked with cross markers. (b) The refined mesh. f1faces are subjected to two h-refinement levels while neighboring faces are subjected to a single h-refinement level. 12 I n t e r f a c e + 1 + 2(a) I n t e r f a c e + 1 + 2 (b) Figure 4: Refinement in the vicinity of zero-measure faces. (a) The control mesh before the refinement. f2faces are marked with cross markers. (b) The refined mesh. f2faces are subjected to two h-refinement levels while neighboring faces are subjected to a single h-refinement level. Notice that interface control points in the refined mesh are located to precisely replicate the interface geometry. 13 (a) (b) (c) (d) (e) (f) Figure 5: The procedure for generating the mechanical model. (a) The refined mesh according to Step 1 in Section 2.2. Unit and zero knot intervals are marked with black and red lines, respectively. The green lines represent the face extensions of the T-junctions in the mesh. (b) A zoom-in to the interface region in (a). (c) The refined mesh after applying the Laplacian smoothing procedure according to Step 2 in Section 2.2. (d) A zoom-in to the interface region in (c). (e) The physical domain of the mechanical model. (f) A zoom-in to the interface region in (e). The extraordinary points in this subfigure are highlighted with green dots, while t-junctions are highlighted by green cross markers. 14 3. Continuity control between the phases One of the advantages of the proposed approach is its ability to control the continuity of the solution field between the different phases, in the direction perpendicular to the interface, without affecting the shape or topology of the interfaces in the geometrical and mechanical models. As elaborated in Section 2, the continuity between the phases that is obtained by the procedure presented so far is C0. However, for certain multiphase engineering problems, different continuity conditions may be required between the phases. For example, in multiphase flow problems, both C−1andC0continuity conditions are required to represent the pressure and the velocity fields, respectively. For the degree 3polynomials that are used in the present paper, four different continuity conditions can be imposed between the phases: C−1,C0,C1, and C2. Accordingly, we denote the meshes that correspond to the different continuity conditions byMC−1,MC0,MC1, andMC2. To achieve the required continuity condition between different phases, the knot intervals and control points near the interface are modified. This process is demonstrated using the example shown in Figure 6a, which depicts the interface region in an MC0mesh that could be obtained according to Section 2. The corresponding physical domain is demonstrated in Figure 6b. The control points lying on the interface between the phases in MC0are denoted by P0, while the control points in the adjacent layers on Ω1andΩ2are denoted by P1 Ω1andP1 Ω2, respectively. The control points in the following layers are denoted in a similar pattern by P2 Ω1,P2 Ω2,P3 Ω1, andP3 Ω2, as shown in Figure 6a. Starting from MC0, the meshes MC−1,MC1, andMC2are obtained as follows: •MC−1: This algorithm replicates all the control points P0inMC0, creating two separate sets of control points P0 Ω1andP0 Ω2(Lines 1-2). Then, it associates each set with its respective phase, leading to two disconnected domains, Ω1andΩ2(Lines 3-4). Algorithm 1: The algorithm for modifying MC0toMC−1. Data:MC0 Result: MC−1 1P0 Ω1=P0; 2P0 Ω2=P0; 3Assign P0 Ω1toΩ1; 4Assign P0 Ω2toΩ2; The resulting MC−1mesh is presented in Figure 6c. The control points of Ω1are kept blue, and the control points of Ω2are marked in red and slightly decreased in size to demonstrate the multiple control points along the interface. The physical domain of MC−1 is demonstrated in Figure 6d. •MC1: This algorithm modifies MC0to have a single layer zero knot intervals separating the different phases. In order to preserve the interface’s geometry, the control points in theP1 Ω1andP1 Ω2layers are repositioned so that their average location matches the control points in the P0layer (lines 5-6, 8-9). This repositioning is based on the distances between 15 the control points in P0,P2 Ω1, andP2 Ω2to prevent overlapping faces in the modified mesh (lines 2-4). Finally, the algorithm removes the control points in the P0layer (Line 10) and connects the control points in the P1 Ω1andP1 Ω2layers with zero knot intervals (lines 10-11). Algorithm 2: The algorithm for modifying MC0toMC1. Data:MC0 Result: MC1 1foreach point in P0do 2 d1=dist(P0,P2 Ω1); 3 d2=dist(P0,P2 Ω2); 4ifd1< d 2then 5 P1 Ω1=mean (P0,P2 Ω1); 6 P1 Ω2= 2P0−P1 Ω1; 7else 8 P1 Ω2=mean (P0,P2 Ω2); 9 P1 Ω1= 2P0−P1 Ω2; 10Remove P0; 11Connect P1 Ω2andP1 Ω2with zero knot intervals. TheMC1mesh after applying this procedure and its corresponding physical domain are demonstrated in Figures 6e and 6f, respectively. In this case, the zero-measure faces between the phases in Figure 6e are used to compute the basis functions that are associated with both phases. Note that the single zero knot intervals along the interface leads to a multiplicity of m= 2in the local knots of the basis functions in the interface region. Following the relation between the continuity, the polynomial degree, p, and the knots multiplicity, m, the obtained continuity along the interface is Cp−m=C3−2=C1. •MC2: This algorithm transforms MC0toMC2by eliminating zero knot intervals in the interface region. In order to preserve the interface’s geometry, the control points in layer P1 Ω2are first positioned at the same locations as the control points in layer P0(line 2). Then, the control points in the layers P2 Ω1andP2 Ω2are repositioned so that their average locationmatchesthecontrolpointsinlayer P0(lines6-7, 9-10). Thisrepositioningisbased on the distances between the control points in P0,P3 Ω1, and P3 Ω2to prevent overlapping faces in the modified mesh (lines 3-4). Finally, the two layers of control points P0andP1 Ω2 are removed, and the control points in layers P2 Ω1andP1 Ω2are connected with a non-zero knot intervals (lines 11-12). 16 Algorithm 3: The algorithm for modifying MC0toMC2. Data:MC0 Result: MC2 1foreach point in P0do 2P1 Ω2=P0; 3 d1=dist(P0,P3 Ω1); 4 d2=dist(P0,P3 Ω2); 5ifd1< d 2then 6 P2 Ω1=mean (P0,P3 Ω1); 7 P2 Ω2= 2P0−P2 Ω1; 8else 9 P2 Ω2=mean (P0,P3 Ω2); 10 P2 Ω1= 2P0−P2 Ω2; 11Remove P0andP1 Ω1; 12Connect P2 Ω1andP1 Ω2with non-zero knot interval. The mesh after applying this procedure and its corresponding physical domain are demonstratedinFigures6gand6h, respectively. Notethatinordertoobtainanadmissible mesh, all the knot intervals between P2 Ω1andP1 Ω2should have the same non-zero value. The primary objective of these algorithms is to control the continuity between phases while preserving the interface geometry. It is worth noting, however, that Algorithms 2 and 3 for generating the MC1andMC2meshes may alter the shape of internal edges in the mesh in the proximity of the interface. Nevertheless, this paper focuses on 2D planar geometries, and the shape preservation of the internal edges is not required. 17 (a) (b) (c) (d) (e) (f) (g) (h) Figure 6: (a) A demonstration of MC0mesh that could be obtained after applying the algorithm in Section 2. (b) The physical domain of the mesh in (a). (c) The MC−1mesh. The replicated control points could be seen along the interface as coinciding blue and red dots. (d) The physical domain of (c). (e) The MC1mesh. The locations of the removed P0layer are marked with blue cross-marks. (f) The physical domain of (e). (g) The MC2mesh. (h) The physical domain of (g). A superscript □Modis added to the modified control points. 18 4. Sensitivity analysis In this section, we present the sensitivity of the mechanical model with respect to the control parameters of the multiphase domain. Since the level-set method is adopted in this paper to construct the multiphase domain, these control parameters are, in fact, the control parameters of the level-set function that are denoted by αin Section 2. For a clear presentation of the sen- sitivities, we partition the procedure for generating the mechanical model into five consecutively related differentiable parameters: 1.α– The control parameters of the level-set function, 2.PG(α)– The control points of the mesh after applying the untrimming procedure in Step 4 in Section 2.1, 3.PGℓ(PG)– The control points of the mesh after applying the Laplacian smoothing in Step 5 in Section 2.1, 4.PM(PGℓ)– The control points of the refined model after applying Step 1 in Section 2.2, 5.PMℓ(PM)– The control points of the final mechanical model after applying Step 2 in Section 2.2. By this partition, the control points that represent the final mechanical model can be repre- sented as: PMℓ(PM(PGℓ(PG(α)))), (6) Correspondingly, the sensitivity of the mechanical model with respect to the control design variables is expressed by: ∂PMℓ ∂α=∂PMℓ ∂PM∂PM ∂PGℓ∂PGℓ ∂PG∂PG ∂α. (7) The first term of Eq. (7) represents the sensitivity due to the Laplacian smoothing procedure in Step 2 in Section 2.2. The second term represents the sensitivity due to the refinement procedure in Step 1 in Section 2.2. The third term represents the sensitivity of the Laplacian smoothing procedure in Step 5 in Section 2.1. The fourth term represents the sensitivity of the untrimmed mesh after Step 4 in Section 2.1 with respect to the control parameters of the level-set function. These four terms are stated mathematically in Equations (5), (4), (3) and (1) and their sensitivities are given by: ∂PMℓ ∂RM=LM, (8) ∂PM ∂PGℓ=T, (9) ∂PGℓ ∂PG=LG,and (10) ∂PG ∂α=A, (11) 19 respectively. Developing these sensitivities is essential for solving engineering problems with moving inter- faces using the proposed approach. Engineering problems that require these sensitivities include topology optimization, multiphase flow, and phase transition (such as solidification and evapora- tion of liquids). Specifically for topology optimization, the sensitivity of the geometrical model, ∂PGℓ ∂α=∂PGℓ ∂PG∂PG ∂α, might be required when dealing with physical or geometrical measures along the interface, such as limiting the maximal curvature of the interface in the optimized design. 5. Examples This section is dedicated to demonstrating the ability of the proposed approach to solve discontinuous two-phase engineering problems. The examples include two-phase problems from different engineering fields. We begin with solving Poisson’s equation, where the convergence and the effect of controlling the continuity between the phases are investigated. In the second example, a linear elasticity problem is solved both for a problem with an analytical solution and for a problem with a free-formed domain. In the third and fourth examples, the ability of the proposed approach to solve engineering problems with a moving interface is examined, where topology optimization and two-phase flow problems are solved. In the last two examples, the analytical sensitivities discussed in Section 4 are utilized to update the evolving interface. In all examples, the domain of the problem and its boundary are denoted by Ωand∂Ω, respectively; The two different phases occupy the two open subdomains Ω1andΩ2such that Ω1∪Ω2=ΩandΩ1∩Ω2=∅; The Ω1phase corresponds to regions with positive values of the level-set function, indicated by the color red in all illustrations, while the Ω2phase corresponds to regions with negative values of the level-set function and colored blue. These phases are separated by the interface Γ =Ω1∩Ω2, as depicted in Figure 7. Figure 7: An illustration of a two-phase domain. 20 5.1.Example 1: Poisson’s equation with different continuity conditions between phases In the first example, we examine the convergence of the proposed approach for different continuity conditions between the phases. For this purpose, we solve the following Poisson’s equation: −∆φ=8kiπ2 l2sin2πx l sin2πy l inΩ, (12) where zero Dirichlet boundary conditions are imposed on the external boundary of the domain, ∂Ω;kirepresents the material characteristics within Ωi; and lis a geometrical parameter de- pending on the dimensions of the domain Ω. In this particular example, Ωhas a square shape, andlis selected to be equal to the square’s side length. The examination is conducted on a model of a square plate with a quarter-circle inclusion made of a different material in one of its corners. The side length of the plate is 2 units, while the radius of the inclusion is 1 unit. The control mesh of the mechanical model is presented in Figure 8a, with extraordinary points marked by green dots and the extensions of T-junctions indicated by green lines. The corresponding physical domain is presented in 8b. The development so far assumes a unified weight for all the basis functions, so the hole is not precisely circular. Instead, it will be a degree 3 B-spline whose control points are chosen to replicate as closely as possible the circular shape. This mechanical model consists of 520Bézier patches and 699control points. (a) (b) Figure 8: The mechanical model for a square plate with a quarter-circle inclusion. (a) The control mesh. (b) The physical domain. First, to investigate convergence under different continuity conditions in the perpendicular direction to the interface, and to compare with an analytical solution, we assume that the two materials have the same characteristics, k1=k2=k. The corresponding analytical solution in this case is: φa=ksin2πx l sin2πy l . (13) The numerical solution, φn, fork= 2with C0continuity between the phases is shown in Figure 21 9a. To enhance distinguishing between the different phases, the solution is presented once more using the original red-blue color scheme for the phases in Figure 9b. (a) (b) Figure 9: Numerical solution of Eq.(12) on the mechanical model shown in Figure 8 for k= 2with C0continuity between phases. (a) Color-scaled representation. (b) Representation with the original red-blue color scheme for the phases. To generate refined meshes, the T-junctions are initially eliminated from the model through local refinement in their vicinity. Subsequently, refined meshes are generated through a series of single h-refinement steps. The numerical solutions over these meshes for k= 2with C0 continuity between the phases are displayed in Figure 10. The numbers of degrees of freedom for all the refinement levels under different continuity conditions between phases are summarized in Table 2. the first row in the Table corresponds to the meshes presented in Figure 10. Note that the number of degrees of freedom in this example is equal to the number of control points in the mechanical model. The corresponding convergence graphs for the L2norm and the H1 semi-norm with respect to the number of degrees of freedom are displayed in Figures 11a and 11b, respectively. RL0RL1RL2RL3 MC0741 2857 11217 44449 MC1716 2808 11120 44256 MC2691 2759 11023 44063 Table 2: The numbers of degrees of freedom for all the refinement levels under different continuity conditions between phases. The first row in this table corresponds to the meshes presented in Figure 10. The abbreviation RLicorresponds to - Refinement Level i. The convergence graphs indicate the convergence of the proposed approach for different continuity conditions between the phases. Even though the continuity between the phases in the analytical solution is ∞, increasing the continuity between the phases in the numerical solution led to a larger error. This is due to the reduction of the number of degrees of freedom that is accompanied with increasing the continuity. Note that as the continuity increases between the different phases, the number of degrees of freedom decreases. This is the result of applying the suggested algorithms in Section 3. It should be noted that developing algorithms that increase the number of degrees of freedom as continuity between phases increases could potentially reduce 22 (a) (b) (c) (d) Figure 10: The numerical solutions over meshes with different refinement levels for k= 2with C0continuity between the phases. Sub-figures (a) to (d) correspond to refinement levels 0 to 3, respectively. The numbers of degrees of freedom associated with these meshes are presented in the first row of Table 2. 1 031 041 05 N u m b e r o f d e g r e e s o f f r e e d o m1 0! 51 0! 41 0! 3j j e j j L2 ( + ) 11 . 2 7MC0 MC1 MC2 (a) 1 031 041 05 N u m b e r o f d e g r e e s o f f r e e d o m1 0! 31 0! 21 0! 1j e j H1 ( + ) 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 10 . 7 6MC0 MC1 MC2 (b) Figure 11: Convergence graphs for the Poisson’s equation with k1=k2= 2and different continuity conditions between the phases. (a) L2norm. (b) H1semi-norm. the error, but this is beyond the scope of the current research as the error values for different continuity conditions are of the same order of magnitude and possess very close convergence rates. Furthermore, the convergence rates illustrated in both convergence graphs align closely with the theoretical expectations detailed in [53]. In mechanical models incorporating extraordinary points, theconvergenceratesdependnotonlyuponpolynomialdegreesbutalsouponthenumber and valence of the extraordinary points within the model. As discussed in Section 2.2, the continuity along the spoke edges of an extraordinary point is reduced from C2toC0. With 23 mesh refinement, the physical region supported by extraordinary points shrinks, causing the physical length of these C0continuity lines to decrease. And this has a direct effect on the convergence rates. For more details on convergence rates in mechanical models that incorporate extraordinary points, the reader is referred to [53]. Convergence for the case where the phases have different material characteristics is conducted in the subsequent example. Finally, we display the solution of Eq. (12) over the mechanical model depicted in Figure 5e. The material parameters are k1= 2andk2= 4and continuity of C0between the phases. This model is a square with a side length of 8 units and consists of 2,516Bézier patches and 3,007 control points. The numerical solution is presented in Figure 12. (a) (b) Figure12: NumericalsolutionofEq.(12)onthemechanicalmodelshowninFigure5efor k= 2withC0continuity between phases. (a) Color-scaled representation. (b) Representation with the original red-blue color scheme for the phases. 24 5.2.Example 2: Static linear elasticity This example demonstrates the effectiveness of the proposed approach in solving static lin- ear elasticity problems within two-phase domains. Additionally, it showcases the approach’s capability to ensure smooth and well-defined stresses along the explicitly defined interface. The example starts with a theoretical overview, followed by a showcase example that includes a com- parison with an analytical solution. Finally, an illustrative example of a free-formed geometry is presented. 5.2.1. Theoretical overview This example deals with static linear elasticity problems while neglecting body forces. The corresponding partial differential equation is: ∇ ·σ=0,inΩ, (14) where σrepresents the stress tensor, which characterizes the internal forces within the material. Foreachmaterialphase, thestress-strainrelationshipisgivenbyHooke’sLaw: σ=Ci:ε, where εis the strain tensor and Ciis elasticity tensor, describing the material’s mechanical properties inΩi. The strain tensor εis related to the displacement field uas follows: ε=1 2 ∇u+ (∇u)T . Here, ucorresponds to the displacement vector field, representing the material’s deformation, and serves as the solution field for Eq. (14). 5.2.2. A disk plate with a circular inclusion Toassesstheeffectivenessoftheproposedapproachinsolvingstaticlinearelasticityproblems within two-phase domains, we first consider a problem with an analytical solution: a circular inclusion of a different material centered within a disk plate. The disk plate is subjected to a radial displacement of ur(r) =r, where rrepresents the radius. The analytical solution for the plain strain case is provided in [54]. The inner and outer radii are set to r1= 1andr2= 3, respectively. By exploiting the tangential symmetry of the problem, we model only a quarter of the disk plate, while the numerical model encompasses a bi-symmetric square region, as depicted in Figure 13. This numerical model covers a square with a side length of a= 2. Symmetric boundary conditions are applied along the two symmetry axes. Precise tractions, based on the analytical solution in [54], are applied to the remaining edges of the model. The Young’s modulus and Poisson’s ratio for the material in Ω1areE1= 2andν1= 0.2, respectively. For the material in Ω2, these values areE2= 1andν2= 0.25. The mechanical model in this example is identical to the one in the previous example, with its control mesh and physical domain presented in Figures 8a and 8b, respectively. To examine the proposed approach, we initially present the numerical solutions for the radial and tangential stresses, σrrandσθθ. These solutions are illustrated in Figures 14a and 14b, respectively. Recall that the solution to this physical problem herein is the displacement field, and its continuity in the perpendicular direction to the interface is C0. This leads to C−1 continuity in the numerical solutions of both σrrandσθθ, in the direction perpendicular to the 25 Figure 13: Problem settings and numerical model for the disk plate with inclusion example. interface. However, the analytical solutions exhibit C0andC−1continuity in the same direction forσrrandσθθ, respectively. Therefore, despite the C−1continuity of the numerical σrrin the direction perpendicular to the interface, we expect to obtain equal values of σrrfrom both sides of the interface. We consider this as an accuracy indicator of this method. In Figure 15, we present numerical solutions of σrralong the interface within the domain 30◦< θ < 60◦, obtained through varying mesh refinement levels. The analytical σrris constant along the interface and it equals 2.09302. This solution is marked with a black line in the Figure. Dashed and continuous colored lines correspond to solutions within Ω1andΩ2respectively. From this graph, we can see that as the mesh is refined the difference between the values of σrrfrom opposing sides of the interface converge. Moreover, with increasing mesh refinement, these values approach proximity to the analytical solution. Moreover, due to the C2displacement field along the tangential direction of the interface, both numerical stresses exhibit C1continuity in that same direction. (a) (b) Figure 14: Numerical stresses over the disk plate with a centered circular inclusion example. (a) Radial stress, σrr; (b) Tangential Stress, σθθ. Additionally, Figure 14 demonstrates the reduced continuity lines near an extraordinary point. As stated in Section 2, the construction of the basis functions in this paper results in a solution field of C0continuity along the spoke edges of an extraordinary point. Consequently, this leads to C−1continuity in the stress field along these edges. However, it is essential to 26 30°40°50°60° 32.09272.09282.09292.0932.0931< r rR L0- 6 9 9 C P s R L1- 2 8 5 7 C P s R L2- 1 1 2 1 7 C P s R L3- 4 4 4 4 9 C P sFigure 15: σrralong the circular interface within the domain 30◦< θ < 60◦. The analytical solution is marked by a solid black line, while dashed colored lines represent the numerical solution within Ω1, and continuous colored lines represent the numerical solution within Ω2. Here, RL and CPs correspond to refinement levels and control points, respectively. emphasize that the extraordinary points are situated at least 3patches away from the interface due to the refinement procedure outlined in Section 2. As a result, the reduced continuity lines surrounding an extraordinary point do not impact the smoothness along the interfaces. To enhance the continuity of the solution field near extraordinary points, the algorithm proposed in [49] can be used where G1continuity for the basis functions in the vicinity of extraordinary points can be obtained. The convergence graph for the L2norm of the solution field in this example with respect to the number of degrees of freedom is displayed in Figure 16. This graph indicates the convergence of the proposed approach, while the convergence rate aligns with their theoretical expectations as presented in [53] and discussed in the previous example. 1 031 041 05 N u m b e r o f d e g r e e s o f f r e e d o m1 0! 71 0! 61 0! 5j j e j j L2 ( + ) 11 . 2 6 Figure 16: Convergence graph of the L2norm for the linear elasticity example. 5.2.3. General free-formed geometry In this subsection, the proposed approach is tested for solving linear elasticity over a free- formed geometry. The examined structure is a simply supported beam with a uniformly dis- tributed load where only the symmetric half part is modeled. The problem settings are depicted 27 in Figure 17. The elasticity moduli for the material in Ω1andΩ2areE1= 10andE2= 5, respectively, while the Poisson’s ratio is set to 0.25. The geometrical model in this example is constructed based on the level-set function shown in Figure 18. The level-set function is defined over a mesh of 15×60bi-cubic Bézier patches and18×63control points. The corresponding geometrical and mechanical models according to Section 2 are presented in Figures 19a and 19b , respectively. A zoomed-in views to the interface regions in the green boxes in Figures 19a and 19b are presented in Figures 19c and 19d, respectively. The geometrical model consists of 1,416control points and 991Bézier patches while the mechanical model consists of 11,243control points and 10,527Bézier patches. Figure 17: Linear elasticity over a free-formed geometry - Problem settings. Figure 18: The level-set function for generating the geometrical model of the simply supported beam. The stress field σxis demonstrated in Figure 20a. One of the benefits of the proposed approach is demonstrated in Figure 20b, which displays a zoomed-in perspective view of the green box in the stress field shown in Figure 20a. This figure indicates to the smoothness of the stresses along the interface in both phases. Furthermore, as previously noted, the solution field maintains C2continuity in the tangential direction and C0continuity in the perpendicular directiontotheinterface. Consequently, thecontinuityofthestressinthosedirectionsis C−1and C1,respectively. Inaddition,Figure20bshowsthereducedcontinuitylinesnearanextraordinary point. To further emphasize the benefits of the proposed approach for solving static linear elasticity problems, we evaluate σxalong the exact interface between the phases. For this purpose, the internal interface on the right-hand side in Figure 19b is selected. The stresses in both phases are presented in the graph in Figure 21. According to Section 2, the interface curve is represented as a periodic cubic B-spline curve. The B-spline curve is extracted according to Appendix B, and 28 (a) (b) (c) (d) Figure 19: (a) The geometrical model of the simply supported beam corresponding to the level-set function in Figure 18. (b) The Mechanical model corresponding to (a). (c) A zoom-in to green square in (a). (d) A zoom-in to green square in (b). (a) (b) Figure 20: (a) The stress field σxover the domain of the beam.. (b) A zoom-in to the green box in Figure 20a. it is marked with a magenta color in the figure. Five points, A-E, along the curve, are chosen to be located on the graphs. We note that the continuity of the displacement field along the 29 interface is C2. This leads to stresses of C1continuity along the interface. 0 5 1 0 1 5 2 0 2 5 3 0 3 5 4 0 9- 0 . 4- 0 . 200 . 20 . 4< x ABCD EAA B CDE AA Figure 21: The stress σxevaluated along the exact interface in both phases. The stress field is C1continuous in the tangential direction and C−1continuous in the perpendicular direction to the interface. 5.3.Example 3: Topology optimization for static linear elasticity In this example, a topology optimization problem that involves static linear elasticity over a two-phase domain is solved. This example shows the ability of the proposed approach to deal with engineering problems that involve moving boundaries and topological changes. The problem settings are depicted in Figure 22. The two-phase domain consists of two materials with different properties. The optimization problem is solved to minimize the compliance of the structure while constraining the amount of the stiffer material, with the larger elasticity modulus. This optimization problem is commonly used to maximize the rigidity of the structure while controlling the amount of the constituent materials. The formulation of the optimization problem is as follows: minimize: αJ(α,u), s.t.: V1(α)/V∗ 1−1≤0, with: K(α)u(α) =f,(15) where Jis the compliance, J=Z Ωϵ(u)Cϵ(u)dΩ. (16) αisthevectoroftheoptimizationdesignvariables,namelythez-coordinatesofthecontrolpoints parametrizing the level-set function; the expression ϵ(u)Cϵ(u)represents the strain energy, whose minimization is equivalent to compliance minimization; V1(α)is the volume of Ω1, and computed by V1(α) =Z Ω11dΩ1; (17) V∗ 1is the maximum allowed volume of Ω1; and K(α)andfare the global stiffness matrix and the global external force vector, respectively. Without loss of generality, we assume for the current demonstration that the force vector fdoes not depend on the design variables α. 30 The sensitivities of the compliance and the volume expression with respect to the control design variables αare analytically computed by: ∂J ∂α=∂J ∂PΓ∂PΓ ∂α,and (18) ∂V1 ∂α=∂V1 ∂PΓ∂PΓ ∂α, (19) respectively, where PΓrepresent the control points along the interface. The derivatives∂J ∂PΓare computed as suggested by Allaire et al. [55], where the compliance derivative is calculated using the adjoint approach. The design variables are updated using the Method of Moving Asymptotes (MMA) as sug- gested in [56]. These variables are bounded within the range of −1to1. External move limits of0.1are enforced on MMA updates to the design variables. The optimization process termi- nates when the absolute difference in both the objective function and the constraint over three consecutive optimization cycles is less than 10−3. Thestructuresolvedinthisexampleisathree-pointbendingbeamwhereonlyonesymmetric half is modeled and optimized. The initial B-spline surface consists of 25×60Bézier patches and is constructed over a 28×63control points mesh leading to 1764optimization design variables. The maximum allowed volume, V∗ 1, is set to 900, which is 60%of the volume of the design domain. The Young’s modulus and Poisson’s ratio for the material in Ω1areE1= 10and ν1= 0.3, respectively. For the material in Ω2, these values are E2= 1andν2= 0.25. The optimization started with the initial level-set configuration that is shown in Figure 23 where the zero-level-set is highlighted with red lines. The initial geometrical model is shown in Figure 24a. This model consists of 3474control points and 2046Bézier patches. The values of the objective function Jand the volume V1corresponding to this topology are J= 107001 and V1= 1025 .4, respectively. The optimization process concluded after 32design cycles, resulting in a final objective function of J= 72256 andV1= 900 .3. This design leads to a 33%reduction in the initial compliance value, while also meeting the volume constraint. Snapshots of the 7th,13th, and the final ( 32nd) design cycles are shown in Figures. 24b-24d, respectively. The convergence graph of the optimization process is displayed in Figure 25 where the objective is scaled according to its value in the initial topology. This graph indicates the stability of the optimization procedure. In the first two iterations, the optimization process prioritizes the fulfillment of the volume constraint, resulting in an increase in the objective function value. Once the volume constraint is satisfied, the reduction of the objective function value begins. The results shown in this example demonstrate the ability of the proposed approach to solve moving interface problems, particularly topology optimization, while maintaining a precise and explicit representation of the interfaces. This paves the way for further exploration of the pro- posed approach in addressing other optimization problems. In particular, those that involve engineering problems that are sensitive to the quality of the boundary representation, such as stress-constrained topology optimization or topology optimization for fluid-structure interac- 31 Figure 22: Problem settings for the topology optimization example. Figure 23: Initial level-set configuration for the topology optimization example. (a) Initial design, J= 107001 ,V1= 1025 .4. (b) Iteration 7,J= 84467 ,V1= 893 .6. (c) Iteration 13,J= 73110 ,V1= 896 .3. (d) Iteration 32,J= 72256 ,V1= 900 .3. Figure 24: Snapshots from the optimization procedure with the corresponding values of the objective, J, and constrained volume, V1. tion problems. Additionally, and particularly beneficial for topology optimization applications, the optimal designs obtained using the proposed approach are compatible with state-of-the-art 32 0 2 5 1 0 1 5 2 0 2 5 3 0 3 2 I t e r a t i o n- 0 . 200 . 20 . 40 . 60 . 6 70 . 811 . 2C o n v e r g e n c e g r a p h J i = J 0 V 1 i = V 1$! 1Figure 25: Convergence graph for the topology optimization example. computer-aided design software due to their spline-based representation. This makes the pro- posed approach able to maintain a consistent representation across geometric, mechanical, and computer-aided design models, in all design cycles eliminating the transitional stages typically associated with adopting different representations for the different models, which can introduce inaccuracies into the optimization procedure. 5.4.Example 4: Stokes equation on two-phase domain In this example, we demonstrate the ability of the proposed approach to address two-phase flow problems by solving the Stokes equation over a domain of two-phase fluids, neglecting the inertia forces. The problem settings of this example are shown in Figure 26. This example is characterized by constantly moving interfaces while the shape and topology of the different fluids evolve based on the solution of the governing equations. To provide a comprehensive understanding of the implementation of this example, we provide a brief overview of the theory and implementation remarks in the following subsection. 5.4.1. Theory and implementation remarks The Stokes equations for two-phase flow are: ∇p−µi∇2u=ρig inΩi, i= 1,2, (20a) ∇ ·u= 0 inΩi, i= 1,2, (20b) where uandpare the variables of the governing equations and represent the velocity and pressure fields, respectively. µiandρiare the dynamic viscosity and the density of the fluids inΩi, respectively; gis the gravitational acceleration field. The boundary conditions on the interface between the phases are: [u] = 0 onΓ, (21a) [σ·nΓ] =−γκnΓ onΓ, (21b) 33 where σis the Cauchy tensor and it is defined by σ=−pI+µ(∇u+ (∇u)T);nΓis the normal vector to the interface Γ;γis the surface tension coefficient; and κis the curvature of Γ. The brackets [□]represent the difference between the physical quantities in the different phases along the interface, e.g. [u] =u1−u2|Γ,ui=u|Ωi, i= 1,2. The boundary conditions on ∂Ωare presented in Figure 26. The term in (21a) ensures a continuous velocity field between the different phases. This condition is achieved by adopting MC0mesh for the velocity field according to Section 3. Con- versely, the pressure field is constructed using MC−1mesh leading to the required C−1continuity between phases. The term in (21b) represents the difference in the Cauchy tensor due to the surface tension force that acts along the interface. To numerically solve the equations in (20), the numerical solutions for the velocities and the pressure fields and their corresponding test functions are constructed using the same bi-cubic spline-based representations following the Bézier extraction procedure discussed in Section 2. To stabilize the pressure field, the Pressure-Stabilizing Petrov-Galerkin approach (PSPG) has been adopted. Denoting the B-spline basis functions by N, the stabilized weak form of (20) can be represented as: Z Ωµ∇N∇uh− ∇NphdΩ =Z ΩNρgdΩ−Z ΓNΓγκnΓdΓ, (22a) Z Ω∇NphdΩ +NpatchX patch =1τPSPG ρZ Ωpatch∇N·[∇ph−µ∇2uh−ρg]dΩpatch= 0,(22b) where uh= (uh x, uh y)andphare the numerical solutions of the velocity and the pressure fields, respectively; Ωpatchis the domain of a single patch in Ω;τPSPGis the stabilization parameter and is computed by τPSPG = 0.1l2 4υ, where lin the mean value of the maximum patch lengths in the two directions and υis the kinematic viscosity, υ=µ/ρ; One of the advantages of the proposed approach for solving two-phase problems is its ability to accurately evaluate physical measures on the exact geometry of the interface. This feature is effectively utilized in the present example to determine the curvature and normal vector of the interface, which are essential for evaluating the surface tension force term in equation (22a). Given the explicit B-spline representation of the interface, these measures can be precisely evaluated as follows: κ=x′y′′−y′x′′ (x′2+y′2)1.5, (23) nΓ=1 (x′2+y′2)0.5" −y′ x′# , (24) where the □′and□′′are the first and second derivatives that can be directly obtained by the derivatives of NΓ. The geometry of the fluids evolves by updating the level-set function based on the velocity field in the perpendicular direction to the interface, uh Γ⊥. Since the interface is explicitly defined, 34 uh Γ⊥can be evaluated at precise locations along the interface as follows: uh Γ⊥=uh Γ·nΓ onΓ. (25) The evolution of the interface is governed by modifying its corresponding control points. De- noting the interface control points by PΓ(PΓ∈PMℓ), their temporal derivatives are computed by: ∂PΓ ∂t=Z ΓNΓuh Γ⊥dΓ. (26) To update the level-set function, whose control parameters are denoted by α, the temporal derivatives∂α ∂tare computed by solving the following system of equations: ∂PΓ ∂α∂α ∂t=∂PΓ ∂t, (27) where the term∂PΓ ∂αis analytically computed following the sensitivity analysis presented in Section 4; and∂α ∂trepresents the unknown variables. This equation ensures finding the level-set function that leads to the desired geometrical changes of the interface according to equations (25) and (26). The system of equations in (27) consists of 2NPΓequations in 2D and Nα variables where NPΓandNαare the number of interface control points in the mechanical model and the number of control parameters αthat support PΓ. The number of equations in (27) does not necessarily equals to the number of variables, leading to either an underdetermined or overdetermined system of equations, depending on NPΓandNα. In the current example, this system of equations is solved using the least squares method. Although the level-set method is commonly used for two-phase flow problems with moving interfaces, it does not naturally guarantee that the amount of fluid in each phase is preserved over time. To ensure that this amount remains equal to its value in the initial configuration, the level-set function is modified at every 5thtime step following Algorithm 4. This algorithm receives as input the areas of phase 2at times t= 0andt=ti, denoted by At0 2andAti 2, re- spectively, along with the level-set controlling parameters at time t=ti, represented by αti and the allowed tolerance, denoted by ϵ. It returns the modified level-set controlling parameters αti modified. The algorithm begins by computing the error between Ati 2andAt0 2(Line 1) and ini- tializing αti modified(Line 2). Then, as the error is larger than a specific threshold ϵ(Line 3) the negative values of αti modifiedare then multiplied by At0 2 Ati 20.5 (Line 4-5). This adjustment leads to increasing or decreasing these values in a manner that minimizes the current error ϵti(Line 6). The tolerance ϵis chosen to be 0.01 based on numerical tests that showed a good balance between the accuracy of the fluid’s area and the efficiency of modifying it. 35 Algorithm4: Thisalgorithmmodifiesthelevel-setfunctiontopreservetheinitial amount of fluids in each phase over time. Data: At0 2,Ati 2,αti,ϵ Result: αti modified 1ϵti=Ati 2−At0 2 At0 2; 2αti modified =αti; 3while absolute (ϵti)> ϵdo 4 idx=indices where α <0; 5 αti modified [idx] =αti modified [idx] At0 2 Ati 20.5 ; 6Evaluate ϵti(αti modified ); 5.4.2. Implementation The procedure presented above is applied to the two-phase flow domain illustrated in Figure 26. The level-set function that is used to generate the initial multiphase domain is shown in Figure 27. This level-set function consists of 18×33control points and 15×30bi-cubic Bézier patches. The corresponding geometrical and mechanical models, along with zoomed-in views of the interface region, are presented in Figure 28. The geometrical model consists of 804 control points and 517 Bézier patches. The mechanical model consists of 8,053 control points and 7,403 Bézier patches. l y = 4 0 : 5 8 l y 0 : 4 2 l y l x = 2+ 10 : 0 6 l y+ 2 . = 1 !n o s l i p , p = 0 n o s l i ps l i p s l i p + 2 d = l x = 3g y = ! 1 xy7 2 = 1 0 ; 2 = 1 0 0 7 1 = 1 0 0 ; 1 = 1 0 0 0! Figure 26: Problem settings for the two-phase flow example. The evolution of the interface is computed over the time interval t∈[0,6], using a uniform time step of 0.025. The velocities and pressure solution fields at t= 0are presented in Figure 29. This figure demonstrates the ability of the presented approach to assign distinct continuity between phases for the different types of solution fields. While using the algorithms is Section 36 Figure 27: The level-set function that is used to generate the initial multiphase domain (a) (b) (c) (d) Figure 28: Two phase flow - initial configuration. (a) Geometrical model. (b) Mechanical model. (c) Zoom-in to the green box in (a). (d) Zoom-in to the green box in (b). 3 the continuity in the perpendicular direction to the interface is C0for the velocity fields and C−1for the pressure field. As stated earlier, the topological evolution of the fluids is governed based on the normal velocity that is evaluated along the precise interface. A demonstration of this velocity is shown in Figure 30. The geometry of the fluids along with contour levels of the corresponding total velocity field at various time points are shown in Figure 31. At time t= 3.25the two domains associated with Ω2merge, demonstrating the ability of the presented approach to deal with topological changes. At time t= 6.0the interface between the different phases flattens, leading to an equilibrium state with very small fluid velocities. 37 (a)uh x(t= 0) (b)uh y(t= 0) (c)ph(t= 0) Figure 29: The Velocities and pressure solution fields at t= 0. (a) The numerical velocity field in the x direction. (b) The numerical velocity field in the y direction. (c) The numerical pressure field with the expected C−1 continuity between the different phases is clearly visible. 38 Figure 30: The velocity in the direction perpendicular to the interface/ (a)t= 0.0 (b)t= 0.5 (c)t= 0.8 (d)t= 1.8 (e)t= 3.2 (f)t= 3.3 (g)t= 3.8 (h)t= 5.5 (i)t= 6.0 Figure 31: The geometry of the fluids and the contour levels of the corresponding total velocity field at various points in time. The color bars in these figures represent the total velocity of the fluids. 39 6. Discussion In most level-set-based applications, the level-set-based topology is typically considered as the final geometrical model. However, due to the implicitness of the level-set-based topology, constructing analysis-suitable mechanical models and evaluating physical measures along the interface often require approximating the geometry. Additionally, since the level-set-based rep- resentation is not directly suitable for analysis, a mesh-based model is reconstructed for analysis purposes, leading to two different representations for the geometrical and mechanical models. A key advantage of the proposed approach is its ability to avoid the undesired implicitness of the level-set-basedtopologybymappingit, usingadifferentiableprocedure, intoanexplicitlydefined spline-based topology that represents the geometrical model. We clarify that the level-set-based topology in this work serves only as an intermediate step in constructing the geometrical model. The mapping procedure in this paper is developed to produce an explicit geometry that can be precisely replicated by the mechanical model. Future implementations might adopt different mapping procedures leading to different geometrical models. However, as long as the mapping procedure is differentiable and does not lead to overlapping patches, it will still be considered valid to replace the mapping procedure used in this paper. The proposed approach relies on refinement techniques to generate analysis-suitable mechan- ical models, specifically to ensure the linear independence of basis functions near extraordinary points. These refinement techniques lead to an increased number of degrees of freedom in re- gions where it is not needed for the accuracy of the numerical solution. Nowadays, there is ongoing research focused on constructing analysis-suitable basis functions that achieve linear independence near extraordinary points without relying on refinement, such as the recent work by Takacs and Toshniwal [57]. Therefore, investigating these new approaches for constructing basis functions in the proposed approach could potentially yield mechanical models with fewer degrees of freedom, while still ensuring suitability for analysis. This paper primarily focuses on solving engineering problems over two-phase domains. Nev- ertheless, it is worth noting that the methodology can be extended to address domains with more than two phases. The representation of such domains using the level-set method is well- established in the literature and is achieved by employing multiple level-set functions simultane- ously to distinguish between different phases [58]. While the untrimming procedure in this paper already encompasses patches that accommodate exactly two phases, this extension should focus on enriching the trimming cases presented in Table 1 with trimmed patches that accommodate 3 and 4 different phases. One of the important aspects that should be examined to assess the effectiveness of the proposed approach is the conditionality of the obtained linear systems. As this problem often arises in approaches involving enhancing the integration near the interfaces, linear systems with high condition numbers may be obtained when Bézier patches with very small sizes are generated near the interface. In this work, the Laplacian smoothing scheme presented in Section 2 is adopted as a remedy for this problem. However, a broader study should be conducted to examine different Laplacian smoothing schemes and their influence on the conditionality of the 40 obtained linear systems. Note that each Bézier patch along the interface layers is defined by 4 layers of control points surrounding the interface curve due to its bi-cubic representation, while the two closest layers to the interface are separated with zero-knots intervals. The conditionality study might examine the effect of choosing different weights for the different layers on the element sizes near the interface and by that on the conditionality of the obtained linear systems. Additionally, choosingdifferentLaplaciansmoothingschemesforthegeometricalandmechanical models should be considered, along with the effect of including the interface control point in the Laplacian smoothing procedure of the geometrical model. One major advantage of the proposed approach is its spline-based interface representation, which is beneficial for both geometric modeling and mechanical simulations. Spline-based repre- sentations are widely used in most CAD applications nowadays, making the proposed approach CAD-compatible and eliminating the need for transition procedures between spline-based CAD models and the different representations used for simulations or other purposes. This eliminates the inaccuracies typically associated with such transition procedures. Furthermore, the explicit feature of the interface allows for a straightforward control on its physical characteristics, such as curvature, perimeter, and slope. 7. Conclusion This paper presents a new approach for solving two-phase engineering problems. This ap- proach is based on untrimming techniques to generate unstructured spline-based representations to define mechanical and geometrical models for geometrically complex two-phase domains. The interface between the phases is precisely and explicitly defined by a set of cubic B-spline curves. This allows for the evaluation of physical measures on points that are precisely located along the interface. The mechanical model precisely replicates the geometrical model using h-refinement techniques. The IGA approach is adopted to solve the engineering problem where all required integrations are performed on the exact physical domain without the need to use approximation techniques. Algorithms that allow for changing the continuity of the solution field between the phases are provided. An analytical sensitivity analysis for generating the geometrical and me- chanical models is developed to allow for solving engineering problems with moving interfaces following the gradient-based approach. The applicability of the proposed approach is demonstrated in four different engineering problems. The results presented in Section 5 reveal the ability of the proposed approach to solve two-phase engineering problems while accurately representing the interface. Specifically, the third and fourth examples demonstrate the proposed approach’s ability to address moving interface problems. In particular, the fourth example illustrates the proposed approach’s ability to solve problems requiring different continuity conditions between phases for different types of unknown solution fields. This paves the way for further studies to examine the advantages of the proposed approach in each two-phase engineering problem separately. Acknowledgments TheauthorwouldliketothankProf. BertJüttlerandassociateProf. Oded Amir for the informative discussions and valuable suggestions during the work on this paper. 41 Funding No funding was received to assist with the preparation of this manuscript. Conflict of interest The author declares that he has no conflict of interest.