A variational multiscale immersed meshfree method for heterogeneous materials

We introduce an immersed meshfree formulation for modeling heterogeneous materials with flexible non-body-fitted discretizations, approximations, and quadrature rules. The interfacial compatibility condition is imposed by a volumetric constraint, which avoids a tedious contour integral for complex material geometry. The proposed immersed approach is formulated under a variational multiscale based formulation, termed the variational multiscale immersed method (VMIM). Under this framework, the solution approximation on either the foreground or the background can be decoupled into coarse-scale and fine-scale in the variational equations, where the fine-scale approximation represents a correction to the residual of the coarse-scale equations. The resulting fine-scale solution leads to a residual-based stabilization in the VMIM discrete equations. The employment of reproducing kernel (RK) approximation for the coarse- and fine-scale variables allows arbitrary order of continuity in the approximation, which is particularly advantageous for modeling heterogeneous materials. The effectiveness of VMIM is demonstrated with several numerical examples, showing accuracy, stability, and discretization efficiency of the proposed method.


Introduction
The accurate description of heterogeneous materials is important in a wide range of engineering problems such as composite solids and structures [1], bi-material heat conduction devices [2], bone structures in biological systems [3], among others. These problems are represented by partial differential equations (PDEs) with discontinuous material properties in space, leading to rough solution with strain discontinuities across material interfaces. The finite element method (FEM) with C 0 continuity employs a bodyfitted mesh across the material interface for proper jumps of strains across material interfaces. However, the construction of such body-fitted meshes is time-consuming, and it becomes non-trivial when material heterogeneity involves complex geometry in multi-dimensional space. As a result, methods that can handle unfitted meshes are attractive in heterogeneous material analysis.
Body-unfitted mesh approaches which employ a boundary integral on the material interface for interfacial compatibility (termed "embedded approach" hereafter) have been proposed under different embedded domain methods. Mortar FEM [4] couples unfitted meshes across the interface by enforcing the interfacial constraint with Lagrange multipliers. Unfortunately, such an approach suffers from the violation of the LBB (Ladyzhenskaya-Babuška-Brezzi, or called inf-sup) conditions [5,6] if the function space of Lagrange multipliers is not properly selected. The violation of the LBB condition may lead to instabilities and oscillations in the interfacial traction fields [7]. Nitsche's method by Stenberg [8] and Hansbo [9] was introduced where the Lagrange multipliers are replaced by the tractions on the interface, and a least-square type stabilization is added with a penalty parameter based on material properties. The determination of the penalty parameter is non-trivial, as a local eigenvalue analysis is often required [10]. To address this issue, the extended finite element method (XFEM) stabilized by bubble functions has been proposed [11,12] under the embedded approach, where the displacement field is decomposed into coarse-and fine-scale fields following the variational multiscale framework [13,14], and the fine-scale field is enriched with element bubble function, similar to the residual-free bubble approach [15]. While the bubble-stabilized XFEM [11] and Nitsche's method [9] capture the interfacial constraint and achieve good convergence, they rely on high order quadrature over the contour integral at the material interface, which is tedious for complex geometric interfaces. Further, the employment of the interfacial constraint under the embedded framework often yields ill-conditioning due to the near-zero internal energy at the overlapping domain (also known as the fictitious domain), in which removal of additional degrees of freedom (DOF) is necessary.
The alternative to an interfacial constraint is to employ a volumetric constraint (enforcement of Dirichlet condition on the overlapping domain, termed "immersed approach" hereafter). Such an approach has been employed with various immersed frameworks, especially in the field of fluid-structure interaction. Similar to the interfacial constraint approach, the penalty method [16] and the Lagrange multiplier method [17] were implemented to enforce such constraint. Glowinski et al. [18] developed the distributed Lagrange multiplier and fictitious domain method for particulate flow, where the kinematics of rigid body motion of particles follows background fluid flow through a perturbed Lagrange multiplier. Zhang et al. [19,20] developed the immersed finite element method (IFEM) for fluid-structure interaction problems, where the compatibility condition of solid/fluid velocity is imposed weakly through a reproducing kernel (RK) function. The volumetric constraint enforcement not only avoids ill-conditioning in the embedded method but also avoids tedious boundary integrals along with the material interfaces, although it still exhibits low order accuracy near the interface. Due to the limited literature on the application of volumetric constraints for heterogeneous solid material problems, the accuracy and robustness of the immersed approach with volumetric constraints for heterogeneous solids deserve further investigation.
Meshfree methods, such as the reproducing kernel particle method (RKPM) [21,22], are attractive alternatives to FEM and are constructed based on a set of scattered points without any mesh connectivity, thus avoiding mesh distortion or mesh entanglement issues [21,23,24]. Further, the flexibility in controlling the local smoothness in the approximation is attractive for problems with varying smoothness in the problem domain. The application of the meshfree method for heterogeneous materials has been studied extensively with techniques such as element free Galerkin (EFG) [25][26][27], meshless local Petrov-Galerkin (MLPG) [2], moving least-squares (MLS) [28][29][30], reproducing kernel particle method (RKPM) [31][32][33][34][35], and Smoothed Particle Galerkin (SPG) [36], to name a few. Due to the lack of Kronecker delta properties, it is not natural for the meshfree methods to impose the compatibility condition, and various techniques have been developed to enforce the constraint strongly or weakly. A meshfree approximation with local enrichment across the material interface [26,29,30,35] was constructed to enable the direct imposition of the compatibility condition, but requires additional DOFs employed on the interface for the enrichment. Wang et al. [32] used a discontinuous Galerkin formulation based on an incompatible patch-based RK approximation, such that additional DOFs are not required since the compatibility and equilibrium conditions were imposed weakly in the variational equation. Wu et al. [33,36] introduced a meshfree discretization strategy for composite materials based on an immersed approach and employed the singular kernel method [37] for approximation such that the point-wise continuity across the material interface was naturally built-in. Koester et al. [34] employed the conforming RK shape function constructed on the boundary-fitted local triangulation, creating C 0 approximation functions at material interfaces. Although meshfree methods have demonstrated their applicability in modeling heterogeneous materials, most of the aforementioned works are based on the embedded framework with interfacial constraints and thus have intrinsic difficulties in modeling heterogeneous materials with complex geometries.
With all the above-mentioned issues and challenges, in this work we present an immersed meshfree formulation under the variational multiscale framework for modeling heterogeneous materials using a volumetric constraint. The proposed method provides flexibility in adopting independent discretizations, approximations, and quadrature rules in the foreground and background domains and avoids the issue of ill-conditioning in the immersed discrete system. The method employs the variational multiscale formulism to enhance accuracy and has built-in stabilization to reduce the oscillations common in conventional immersed methods. One of the proposed methods with an approximated finescale solution requires higher order continuity in the background approximation, and methods such as isogeometric analysis (IGA) [38] and RKPM offer such regularity. In this paper, RKPM is employed due to its flexibility in adopting arbitrary smoothness and domain discretization in solving the variational multiscale equations. Interested readers can refer to [39][40][41] for a comprehensive review of RKPM and its applications.
The paper is organized as follows. The problem description of the heterogeneous, linear elasticity problem and the immersed weak form based on a volumetric constraint are introduced in Sect. 2. The variational multiscale immersed method with different approaches for obtaining the fine-scale and coarse-scale solutions are presented in Sect. 3. Domain integration and stabilization techniques for the foreground and background weak forms are also discussed in this section. Finally, benchmark problems are studied in Sect. 4 to examine the performance of the proposed method. The findings of this research are summarized in Sect. 5.

Strong form
Consider a solid occupying closed subdomains Ω 1 and Ω 2 , where Ω =Ω 1 ∪Ω 2 describes the total domain, Ω 1 ∩Ω 2 = Γ denotes the material interface, and Ω 2 N and Ω 2 D . denote the Neumann and Dirichlet boundary on boundary Ω 2 , respectively, with Ω 2 N ∪ Ω 2 D = Ω 2 . and Ω 2 N ∩ Ω 2 D = � . An illustration of the heterogeneous problem is given in Fig. 1a. To aid the development of the immersed formulation, we intruce a foreground domain Ω + =Ω 1 and a back- The strong form of the problem can be described by the following boundary value problem (BVP): where here u(x) is the displacement vector, (u) is the strain tensor, n 0 is the surface normal vector corresponding to boundary Ω , b is the body force vector, t and g denote the prescribed traction and displacement vectors on Ω N and Ω D , respectively, C is the elastic stiffness, and is the Cauchy stress tensor. The continuity conditions on the interface Γ are as follows: where [[⋅]] the jump operator is defined in (8) and n is the surface normal on the interface.

Remark 2.1.1
Here we classify the boundary value problem with heterogeneous materials as the problem with varying material constants; that is, C(x) is piecewise constant for the problems under consideration as described in Eq. (5). However, the proposed formation is applicable to problems where C(x) is a piecewise higher-order function.

Weak form
It can be easily seen that the conventional discretization of weak form for the setting in Fig. 1 requires a geometric conforming domain discretization and numerical integration for subdomains Ω 1 and Ω 2 , which is tedious for problems involving complex geometry or evolving interfaces. In the immersed setting of Fig. 1b the weak form for this problem is first constructed as follows: where (⋅, ⋅) denotes the L 2 inner product, and u + and u − are the displacement fields on the foreground and background domains, respectively. L − (⋅) = − ⋅ (C − ∶ (⋅)), L + (⋅) = − ⋅ C + ∶ (⋅) a r e stress divergence operators acting on the foreground and background subdomains, respectively. The imposition of the continuity jump condition (7) on the interface Γ in Eq. (9) yields the weak form: The function spaces U + , U − , V + , and V − are given in Eq. (11) where n d is the number of spatial dimensions. The weak form associated with Eq. (10) and interfacial condition u − = u + on Γ is usually referred to as the weak form for the "embedded method," which employs the constraint " u − = u + " on the material interface when solving the composite material system. However, the employment of an interfacial constraint requires a contour integral along material interfaces, which is tedious for complex material geometry. To remedy this issue, a volumetric constraint approach can be introduced [20] as given in the next section.

Volumetric constraint
Due to the redundancy of u − in Ω + , a volumetric constraint for the overlapping domain is employed in the immersed method [17]: The employment of a volumetric constraint (12) for Eq. (10) avoids the need for front tracking and contour integrals along material interfaces [4]. With the constraint (12), the internal energy and body forces terms in Eq. (10) can be rewritten as From Eqs. (13) and (14), the new weak form reads: The function spaces U + , V + , U − , and V − here are given as (11) The imposition of the volumetric constraint in Eq. (12) can be achieved by various methods. Here three common approaches are introduced to impose the volumetric constraint: penalty method, Lagrange multiplier method, and Nitsche's method.

Penalty method
The penalty parameter has been proposed in [9] as where nor is the normalized penalty parameter, E is the Young's modulus, and h is the nodal spacing.  u)]] employed in the weak form, the function spaces Ũ + , Ũ − , Ṽ + and Ṽ − are taken as

Lagrange multiplier method
The three methods introduced in Eqs. (17) to (21) are regarded as conventional immersed methods and their performance will be investigated and compared to the proposed VMIM in the numerical examples.

Approximation and discretization
In this work we employ the Reproducing Kernel (RK) approximation [43] for its flexibility in controlling order of continuity, order of consistency, and locality for the foreground and background displacement fields. The RK approximation is constructed based upon a set of scattered points x 1 , x 2 , … x NP in the domain Ω as shown in Fig. 2, where x I is the position vector of node I, and NP is the total number of nodes. The RK approximation of displacement u i is expressed as where x is the spatial coordinates, u iI is the associated nodal coefficient to be determined, and I (x) is the reproducing kernel (RK) shape function of node I expressed as: where the basis vector H x − x I is defined as and M(x) is the so-called moment matrix: The set G x = I| a x − x I ≠ 0 shown in Eqs. (22) and (25) contains the nodal indexes of point x ′ s neighbors, and a x − x I is the kernel function centered at x I with compact support size a I defined as In the above equation, c is the normalized support size, and h I is the nodal spacing associated with nodal point x I . The kernel function controls the order of continuity and locality of the approximation as shown in Fig. 3, where the C 0 tent kernel function is compared with the following C 2 cubic B-spline kernel function: Tent Function ( C 0 ): Cubic B-spline kernel function ( C 2 ): for z I > 1 . By construction, the RK shape functions satisfy the following nth order reproducing conditions: where n is the specified order of completeness, which determines the order of consistency in the solution of PDEs.

Variational multiscale immersed method
To effectively account for the fine-scale features in heterogeneous materials, a variational multiscale immersed method formulated under the variational multiscale (VMS) formulation originally proposed by Hughes et al. [13,14] is employed. Consider the following decomposition of displacement fields u −h (x) and u +h (x) on the background and foreground domains, respectively, as Here superposed "-"and "^" denote coarse and fine scale components, respectively. Here we consider a scale decomposition for the background displacement field, while the foreground displacement field is represented by the coarsescale component only.

Remark 3.1
Note that the scale decomposition in the foreground or background domain can be introduced based on the material composition in heterogeneous materials. The fine-scale field can also be introduced to the foreground domain Ω + or both foreground and background domains.
The following boundary conditions are imposed for the coarse-scale and fine-scale solution: where Ū − and Û − are the function spaces for the coarse-scale and fine-scale, respectively, and U − =Ū − ⊕Û − . Following [13], we assume the Dirichlet and Neumann boundary conditions on Ω − are fully satisfied by the background coarse-scale displacement, this results in the background fine-scale homogeneous Dirichlet boundary conditions for all boundaries.

Variational multiscale immersed weak form
By introducing the multiscale decomposition in the Eq. (30) into the immersed weak form (19) with Galerkin approximation, the variational multiscale immersed Galerkin equation reads: By the virtue of the linear independence of test functions at different scales, Eq. (35) can be decoupled into the coarsescale equations: and the fine-scale equation: From Eq. (39), it can be seen that the fine-scale solution serves as a correction to the residual of the coarse-scale solution.
By introducing the boundary condition of the fine-scale solution in (32) and (34), and with integration-by-parts operations, the fine-scale Eq. (39) yields Since the domain Ω + is immersed in Ω − , the problem of Eq. (40) can be rewritten as The solution of Eq. (41) can be solved with the aid of Green's function [14], or by an approximation method to be discussed in the following sections.
Similarly, a different form can be employed for the coarse-scale background Eq. (36). Performing integrationby-parts on Eq. (36) and employing the fine-scale boundary conditions in (32) results in the following expression for the coarse-scale solution: From Eq. (42), the fine-scale solution û −h can be directly employed without taking additional differentiation.  (40) can both be used for solving the fine-scale solution. In this study, two approaches are introduced to solve the multiscale solutions. In the first approach, the coarse-scale Eq. (36) and the finescale Eq. (39) are solved by introducing the residual-free bubble method [15] for the fine-scale solution. In the second approach, the coarse-scale Eq. (42) is solved by the standard Galerkin RKPM with stabilized nodal integration, and the fine-scale Eq. (40) is solved by a collocation-type method with an approximation of the local Green's function. Both approaches are introduced in the following sections. (40) and (42) require a second order gradient performed on ū −h , the derivatives on fine-scale terms û −h are avoided, which enables the fine-scale features to be directly introduced in the static condensation process. In addition, since no numerical derivative is performed on û −h , it is more straightforward to employ the resolved fine-scale solution in the coarse-scale equation. The second order gradient L − ū −h required for ū −h can be easily handled by RKPM as shown later.

Residual-free bubble method
One option to obtain the fine-scale solution in Eq. (39) is to employ the residual-free bubble method [15], where the û −h can be approximated by an enrichment function [11] or so-called bubble function [15]. The bubble basis can be computed from the integral of the local Green's function [13]. In the beginning, let the domain Ω − be subdivided into subdomains Ω − L , The fine-scale solution is approximated by bubble basis function N − L in each subdomain Ω − L as [15]: and where S − = I|x I ∈ Ω − is the node set containing all nodes within domain Ω − , û − L is the approximation coefficient, and N − L is the nodal bubble basis function for subdomain Ω − L . In this study we employ the cubic B-spline kernel function in Eq. (28) as the nodal bubble function: The strain vector ε (in Voigt notation) for û − can then be written as where Remark 3.2. 1 The construction of nodal bubble functions can be done with conforming nodal representative domains such as Voronoi cells. Constructing such bubble functions on Voronoi cells with arbitrary geometry can be achieved by, for example, the conforming reproducing kernel [34]. In this work, we relax this condition by introducing cubic B-spline functions on the non-conforming circular nodal domain (45) to construct the bubble functions. The fine-scale homogeneous Dirichlet and Neumann boundary conditions are satisfied by using RK kernel functions as bubble functions with small support size not covering neighboring nodes.
The coarse-scale approximation for coarse-scale variables ū − , ū + and λ are introduced as follows: where N − I , N + I and N I are shape function for ū − , ū + and λ, respectively, B − I and B + I are gradient matrix for N − I and N + I , respectively, S + = I|x I ∈ Ω + is the node set containing all nodes within domain Ω + . Standard C 0 finite element shape functions or smooth RK shape functions given in Sect. 2 can be used as the approximation in Eqs. (48) to (50). More discussion on the choice of shape functions is given in the numerical examples in Sect. 4.
The Galerkin equation of the VMIM reads: find The function spaces and W here are given as By using the approximation in Eqs. (43) to (51) With the static condensation of Eq. (60) into Eqs. (57) to (59), we have:  81) where the term K −1 in G * IJ plays the role of residual-based stabilization. In the conventional augmented Lagrange multiplier formulation, such stabilization is introduced through a penalization.

Remark 3.2.3
The left hand side matrix in Eq. (80) is a symmetric matrix which is expressed entirely in terms of the coarse-scale solution. The fine-scale solution is obtained from Eq. (60) using the coarse scale solution in Eq. (80).

Remark 3.2.4
In the numerical implementation, the kinematically admissible course-scale test and trial functions defined in (61) have been approximated by the employment of Nitsche's method to weakly impose the Dirichlet boundary conditions in (31)- (34). Since the focus of this study is on the imposition of volumetric constraint and the variational multiscale decomposition, we do not include the Ntische's terms in the Galerkin weak form to avoid distraction. The inclusion of Nitsche's terms on the Dirichlet boundary does not affect the derivation of the proposed variational multiscale immersed method.

Approximated fine-scale solution
The second method introduced here follows the concept of stabilized Galerkin methods [14] such as stabilized upwind Petrov-Galerkin (SUPG) [44], where the fine-scale solution can be derived approximately. Due to the strong form nature on the right hand side of Eq. (40), a collocation-type method for Eq. (40) is introduced as follows: where ̂− is the nodal domain average of fine-scale Green functions from Eq. (41), which will be derived later.
The stress divergence operator ⋅ for the coarse-scale solution ū −h is computed as, with C − following the Voigt notation here. Note that (92) requires 2nd order derivatives of the shape functions, which can introduce significant computational costs when using direct derivatives. Here, the implicit gradient technique [45] is employed for N − I,i to avoid a direct derivative, and N − I,ij is then computed by the differentiation of the first-order implicit gradient N − I,i with respect to x j . Following [45], the implicit gradient in N − I,i is computed as: and D i takes the following form and H x − x I , M(x) , and a x − x I are given in Eqs. (24), (25) and (28), respectively. Finally, by introducing the coarse-scale shape functions in Eqs. (48) to (50) into the weak form (86) to (88), and using the arbitrariness of the test functions, the matrix form of VMIM via the approximated fine-scale solution is obtained as The remaining step in this approach is the estimate of the parameter ̂− following [46]: where h is the characteristic nodal distance, E − is the Young's modulus of background matrix media, and c is a constant usually taken as c ∈ [0, 1] . The parameter ̂− is  which yields a Petrov-Galerkin formulation for the given immersed setting. with the operator = L − chosen to be the Galerkin/leastsquares (GLS) operator [48]. The second term on the lefthand side of Eq. (108) can be regarded as a residual-based stabilization.

Remark 3.3.4
From the matrix form given in Eq. (95), method 2 (approximated fine-scale solution) can also be regarded as a stabilized Lagrange multiplier method, similar to the method 1 (bubble method).

Quadrature rule in immersed framework
A stabilized nodal integration scheme is employed for the proposed variational multiscale immersed method, where the stabilized conforming nodal integration (SCNI) [49] with naturally stabilized nodal integration (NSNI) [50] is employed for the coarse-scale equations. For comparison purposes, the same integration formulation is employed for the conventional immersed method. With this integration method, the stiffness matrices in (61) and (62) are expressed as where ̃− I,i x L is the smoothed gradient at the nodal point x L calculated following Chen et al. [49]:  Fig. 4. Same computation is applied to B + I . It was shown in [49] that with the smoothed gradient of shape function in Eq. (114), the first order integration constraint is exactly satisfied. As discussed in [51], in order to maintain linear consistency of the smoothed gradient of a linearly consistent shape function, a simple one-point Gauss integration rule can be used for the contour integral in Eq. (114). In this study, a Voronoi tessellation is employed to generate each nodal integration cell, as shown in Fig. 4. Spurious low energy modes can still be triggered with SCNI, and therefore the naturally stabilized nodal integration (NSNI) [50] is employed. Using this technique engenders additional stabilization terms K −NSNI where ̄− ∇ Ii,j is calculated by the direct differentiation of the first-order implicit gradient ̄− ∇ Ii with respect to x j [50], as shown in Eqs. (93) and (94). Same computation is applied (115) and (116) is the second moment of inertia in each nodal integration domain, calculated from Eq. (118) The matrices K −NSNI IJ and K +NSNI IJ are introduced to maintain coercivity of the system due to reduced nodal integration.
For the approach using the residual free bubble method (Sect. 3.2), SNNI [52] (the non-conforming version of subdomain integrated SCNI [53]) is employed for the domain integration when computing K − IJ in Eq. (69) since the local bubble functions are constructed over non-conforming domains as explained in Sect. 3.2 For all other domain integration terms, direct nodal integration is employed:    where V L is the volume for Lth nodal representative domain.

Numerical examples
In this study, numerical examples are solved to examine the performance of the proposed variational multiscale immersed method (VMIM) for modeling heterogeneous materials. The reproducing kernel approximation with a linear basis and cubic B-spline kernel is adopted for the approximation of the coarse-scale foreground and background displacement fields and the Lagrange multipliers as shown in Eqs. (48) to (50). Circular support with a normalized support size equaling 2.0 is employed for all approximations. In VMIM, the fine-scale solution is approximated by the bubble functions and the approximated Green's functions discussed in Sects. 3.2-3.3. The background Dirichlet boundary conditions over Ω − D are imposed by Nitsche's method with a normalized penalty parameter of 100.
The following methods with their abbreviations are given and tested:   Fig. 9 Numerical solution of scalar variable field along x 2 = 0 by the immersed method with penalty, Lagrange multiplier, and Nitsche's method, compared with the exact solution To assess the accuracy of different numerical schemes, the following normalized displacement and energy norms are employed: in which the superscript "exact" denotes the exact solutions. u −h , −h , −h are the numerical solutions for displacement, strain, and stress from the background domain, respectively, and u +h , +h , +h are the numerical solutions for displacement, strain, and stress from the foreground domain,

Patch test
In the first example, the linear patch test is analyzed to verify the accuracy using linear bases in the RK approximation with nodal integration for the derived variational multiscale immersed formulation. An elasticity equation is considered with the exact solution defined as a linear polynomial function: where the problem contains a square foreground inclusion domain immersed in a square background matrix domain as shown in Fig. 6.
For the patch test, the elastic moduli of the foreground and background are set to be C + = C − , with Young's modulus E + = E − = 1 × 10 3 and Poisson's ratio + = v − = 0.3 , exact is the exact strain with exact = [0.1, 0.1, 0.175] T ; g = u exact is enforced on Ω − D , and the body force is set to b = 0. Next, the solution accuracy with discretization refinement under both uniform and non-uniform discretizations is examined. As shown in Fig. 7, non-uniform discretizations are generated by introducing random numbers between ±0.5h into the uniform discretizations, where h is the nodal distance in the uniform discretizations of a unit square. As shown in Tables 1, 2, 3 and 4, the L 2 and energy norm in uniform and nonuniform discretizations show that all tested methods pass the linear patch test, which demonstrates that the quadrature scheme for both the conventional immersed approach and the variational immersed approach all satisfy the integration constraint.

Heterogeneous material diffusion problems
In the second numerical example, the heterogeneous material diffusion problem given in [54] is tested, where the strong form of the heterogeneous material diffusion problem is given below where u is the scalar variable field. 2 = ⋅ is the Laplacian operator, s is the source term, and q and g denote the prescribed flux and essential boundary condition on Ω N and Ω D , respectively. f is the flux vector. Similarly, the source term s and diffusivity k can exhibit discontinuities across different subdomains, but have smooth distribution in each subdomain Ω + and Ω − �Ω + respectively, given as where k + and k − are scalar diffusivities corresponding to different materials in Ω + and Ω − �Ω + respectively, and s + and s − are the source terms corresponding to Ω + and Ω − �Ω + respectively. In this study, k + and k − are taken as material constants. The following conditions apply on the interface Γ: and the divergence operator can also be defined for each material subdomain:   For this example, the diffusivity is set to be k + = 100 and k − = 1 and the problem is subjected to a source term s = s + = s − = −4 . The problem domain and the discretization are illustrated in Fig. 8. Since this is a scalar problem, the VMIM formulation for this heterogeneous diffusion problem is slightly modified from that in Sect. 3 and is summarized in the "Appendix".
The exact solution of the problem reads: where r 0 = 0.5 is the radius of the circular inclusion. The exact solution g = u exact is employed as the essential boundary condition at boundary Ω − D . The performance of the conventional immersed formulation with different kinematic constraint methods is investigated. As shown in Fig. 9, large error is observed in the scalar solution inside the region Ω + for immersed method with all three constraint enforcement methods. Also, it can be seen in Fig. 10 that flux field exhibits strong oscillation. The numerical investigation shows that the conventional immersed approach yields an unstable and inaccurate solution with each constraint enforcement method.
On the other hand, from Fig. 11, by employing the variational multiscale immersed method, the solution becomes much more accurate. Also, the gradient field becomes much Fig. 22 Convergence rate and computational cost of all tested methods. The values in the legends indicates the average convergence rate more stable near the interface, where the oscillation in the flux field is largely reduced. As shown in Fig. 12, only marginal perturbation from the exact solution near the material interface is found. The introduction of fine-scale features not only increases the accuracy of the solution but also surpresses the numerical instability.
The convergence rate and computational costs for all tested methods are given in Fig. 13. From the convergence study, the VMIM analyses show significantly improved convergence as compared to the conventional IM analyses for both L 2 and energy norms. Additionally, the VMIM runs result in similar computational costs compared to the IM runs as shown in Fig. 13b. The computational cost of VMIM with the bubble method (method 1) is higher than that of the approximate fine-scale approach (method 2) due to the domain subdivision required to numerically integrate the fine-scale equation. Finally, the distribution of solution, its gradient and flux field are shown in Fig. 14, where it can be clearly seen that the flux fields from IM + LM are oscillatory near the interface, while the VMIM results are much more stable, although a minor perturbation from the exact solution is found in the flux field. This perturbation can be further reduced by copying the discretization of the foreground domain into the background domain, as shown in Fig. 15. Improvements are seen in the both VMIM and IM methods. Note that boolean operations provide a straigtfoward means to construct such a discretization by removing the background points and adding the foreground points in the overlapping domain as shown in the left column of Fig. 15.

Displacement
Finally, the fine-scale solutions from the VMIM solves are shown in Fig. 16, where the peak value lies near the material interfaces. From the theory of the variational multiscale method [13], the fine-scale solution plays a role in enhancing the coarse-scale solution. Here, it can be verified that the fine-scale solution is magnified near the interface, which is where the coarse-scale solution exhibits higher errors. This indicates that the fine-scale solution can be used as a posteriori error indicator consistent with the multiscale theory [13].

Circular inclusion in an infinite plate subject to far-field traction
A heterogeneous, linear elasticity problem with the strong form given in Eqs. (1) to (5)  inclusion with radius R in an infinite plate is subject to a far-field traction P as shown in Fig. 17. The computational domain is a truncated finite square domain with length L . The Young's modulus E + = 1 × 10 5 is used for the inclusion and E − = 1 × 10 3 is used for the background matrix media. A Poisson ratio v + = v − = 0.3 is used for both domains. The plane stress condition is assumed and the exact solution [55] in the background matrix domain reads where (r, ) are polar coordinates and 1 , 1 and 1 are expressed as Kernel function used for the background RK approximation  The exact solution is prescribed as essential boundary conditions. The displacement and strain-stress solutions for the immersed methods are shown in Figs. 18 and 19, respectively. It is seen from Fig. 19b that the stress solutions exhibit strong oscillations, similar to the case observed in the heterogeneous material diffusion problem. The VMIM solutions are presented in Figs. 20 and 21 and exhibit much better accuracy. The oscillation in the stress field is largely reduced by the VMIM. The introduction of fine-scale features not only increases the accuracy of the solution but also reduces the numerical instability.
The convergence rates and computational costs for all methods are compared in Fig. 22. From the convergence rate results in Fig. 22a, the VMIMs exhibit significantly improved convergence rates for both L 2 and energy norm compared to the IM methods. Superior computational efficiency (i.e., the accuracy versus total runtime) of VMIM over IM is also shown in Fig. 22b. Finally, the distributions of displacement, strain, and stress fields are shown in Figs. 23 and 24, where it can be clearly seen that the stress solutions from IM + LM are oscillatory, while the results from VMIM solutions are much more stable.
Next, the influence of the foreground and background RK shape functions constructed from kernel functions with different orders of continuity are investigated. In this study, VMIM(App) is employed due to its better accuracy than VMIM(Bubble), as have shown in Fig. 22. Three different kernel functions are adopted: C 0 linear B-spline kernel, C 1 quadratic B-spline kernel, and C 2 cubic B-Spline kernel (the exact expression of the linear and quadratic B-spline kernel can be found from [41]). Different combinations of kernel functions for the background and foreground approximation are tested here.
The stress distributions under different combinations of kernel functions are shown in Figs. 25 and 26. The solution with the C 0 linear B-spline kernel used in the background approximation is worse than all other cases. In contrast, adopting a smoother kernel for the foreground and background approximation gives the most stable result. This is due to the employment of the strong form approach in obtaining the approximated fine-scale solution.
To enhance the accuracy for VMIM, the discretization strategy employed in diffusion problem (left columns of Fig. 15) is also introduced here. As shown in Fig. 27, such modified discretization improves the solution accuracy in all methods, but VMIM again shows better accuracy and stability compared to conventional immersed methods.
The fine-scale solution from VMIMs is shown in Fig. 28  interface. Since the computational domain traction is applied in the x 1 direction, the fine-scale solution is expected to have higher contribution along the x 1 direction.

Modeling of heterogeneous microstructure
Heterogeneous microstructure usually involves complex geometry which leads to a low-quality mesh for FEM, unless a very fine mesh is used, as shown in Fig. 29. The proposed immersed meshfree framework can more effectively handle the discretization for such complicated geometries while providing accurate and stable numerical solutions. In this numerical example, a microstructure with prescribed horizontal and vertical displacements to generate tensile and shear deformation, as shown in Fig. 30, is considered. The geometry of the microstructure is generated by randomly distributing 35 circular inclusions inside a square domain with unit length. The radii of the inclusions are also randomly generated from 0.05 to 0.15 units as shown in Fig. 30. The Young's modulus and Poisson ratio for the The plane stress condition is assumed. The numerical discretization for the proposed immersed method is shown in Fig. 31, where 1,108 nodes are used for the inclusions, and 1681 nodes generated from a 41 by 41 rectilinear grid are used for the background domain. The finite element discretization with a conforming mesh is also shown in Fig. 31. In this test, only VMIM with an approximated fine-scale solution is employed due to its better efficiency and accuracy as shown in the previous examples. A body-fitted FEM result with a highly refined mesh shown in Fig. 31 is employed as the reference solution. Note that when generating the FEM results in this example, a mesh size of 0.005 (corresponding to 44,258 nodes in the FEM solution) was required to avoid a low-quality mesh. In comparison, the proposed VMIM formulation only required 2,789 nodes.
The results of the microstructure under tensile deformation can be found in Figs. 32 and 33, where it can be easily seen that the VMIM results show consistent displacement and stress distribution in the microstructure when compared to reference FEM solution. Also, VMIM produces non-oscillatory strain and stress solutions, and the results are in good agreement with the reference conforming FEM

Conforming FEM
solutions even at a discretization that is 16 times coarser than the conformal mesh.
The results of the microstructure under shear deformation can be found in Figs. 34 and 35. Similar to the tensile test, the VMIM produces smooth and comparable accuracy in comparison with the reference FEM solution obtained with a much finer discretization. Oscillation is observed in the stress solution in the conventional immersed method. For both the tensile and shear tests the strain concentration for nearby inclusions is smoother in the case of VMIM. This may be due to the smooth approximation employed in RKPM, and the coarser discretization used in the VMIM compared to that used in FEM.
Finally, we present numerical results generated by copying the foreground discretization into the background discretization as before. The discretization is shown in Fig. 36 and the numerical results are shown in Fig. 37. The stress and strain field between inclusions are slightly better captured in this case compared to the solution of the other discretization in Fig. 35. However, it is observed that the concave geometry of the foreground inclusions causes local support overlap issues in the region with strong concavity. This issue can be enhanced by the conforming RK method [34] for the foreground inclusion domain.

Summary
In this study, an immersed formulation based on the variational multiscale framework has been proposed for modeling heterogeneous materials. The weak formulation of the problem has been formulated under an immersed setting with independent foreground and background discretizations, and volumetric kinematics constraints have been imposed in the   (2) approximated fine-scale solution. In the first method, the fine-scale equation was obtained by the introduction of nodal bubble functions. The substitution of the fine-scale solution into the coarse-scale equations naturally yields a stabilized Galerkin formulation in an immersed setting. In the second method, due to the stress divergence operator involved in the proposed VMIM with approximated fine-scale field, reproducing kernel approximations with higher-order continuity have been introduced for the solution approximation. Domain integration of the foreground and background weak forms has been accomplished via the SCNI nodal integration technique with NSNI stabilization. The effectiveness of the proposed variational multiscale immersed meshfree framework has been validated by solving several numerical examples. VMIM showed enhanced accuracy and stability compared to the conventional immersed methods with comparable efficiency. The proposed method with an approximated fine-scale solution (method 2) exhibited better efficiency than the residual-free bubble method (method 1). Furthermore, the proposed VMIM approach has been effectively applied to material microstructures with complex geometry. The accuracy of the proposed method is shown to be in good agreement with the solution obtained from FEM with a highly refined body-fitted mesh. In conclusion, the proposed VMIM has been shown to be an effective method for modeling heterogeneous materials when compared to FEM based methods, which require tedious body-fitted mesh generation. Further, VMIM demonstrated superior performance to conventional immersed approaches. The extension of this research to fluid-structure interaction problems is currently in progress.

Appendix: VMIM for heterogeneous material diffusion problems
The immersed weak form for the given heterogeneous material diffusion problem reads: find u + , u − ∈ U + × U − , u − = u + on Ω + , s u c h t h a t ∀ u + , u − ∈ V + × V − : The function space U + , U − , V + , and V − are given as (142)    (178) − = − 1 12 adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.