Interphases and surfaces play an essential role in heterogeneous materials and structures [1][2][3]. This role has become even more important with recent advances in surface chemistry [4,5] and nanotechnology [6,7] where it is now possible to create nanoscale-sized high-performance thin films and coatings with tailor-made properties. Therefore, there exists a broad interest in advanced modeling and simulation technologies that enable their proactive design. In this context, the multiscale nature of these materials whose characteristic length scales vary across many orders of magnitude constitutes a significant challenge. On the one hand, the direct numerical resolution of thin layers with continuum-based standard finite element or boundary element methods requires prohibitively fine mesh sizes [8]. On the other hand, molecular dynamics type approaches are often prohibitively expensive for sample sizes and at scales of engineering interest [9,10].
Asymptotic models are a classical attempt at circumventing this problem. They replace volumes of thin layers by dimensionally reduced surfaces, along which appropriate jump conditions of the fields are formulated. We note that the term "interphase" refers to the physical volume of the layer and the term "interface" to its surface model. The jump conditions for interface models are typically derived by expanding the fields inside the layer with respect to the (non-dimensionalized) thickness in a Taylor series that is then truncated to first-order accuracy [11][12][13][14][15][16][17][18][19][20][21]. Benveniste and Miloh [22] proposed a set of more elaborate second-order models for 2D isotropic elasticity, identifying seven distinct interface regimes. Based on the work of Gurtin and Murdoch [23][24][25] and Steigmann and Ogden [26,27], more general theories of material surfaces have been developed that incorporate the effect of surface tension [3,28,29]. They can be related to specific regimes of the Benveniste-Miloh classification, if the surface tension is neglected and the properties and the thickness of the thin layer are appropriately chosen [30][31][32][33]. Following [11,12,34], we will use the general term "imperfect" interface that covers both classes of problems.
The solutions of general imperfect interface models involve non-smooth functions with jump conditions at arbitrarily complex surfaces that are expressed in terms of higher-order differential operators. This constitutes a significant challenge for their efficient discretization, as standard C 0continuous finite element methods are unable to represent strong discontinuities and do not allow higher-order derivatives. Discretizing higher-order imperfect interface models with C 0 -continuous basis functions has been shown to work for some test cases [35], but is a variational crime in the sense of Strang and Fix [36] that can potentially lead to uncontrolled errors in the numerical solution fields. Instead, research on the numerical treatment of imperfect interface models has mainly focused on first order models that do not contain higher-order derivatives [37,38]. For instance, spring-and membrane-type models for planar and curved three-dimensional interfaces have been integrated into the extended finite element method [39][40][41]. Javili, Steinmann and collaborators developed a thermodynamically consistent theory for general imperfect interfaces that is open to be discretized with standard C 0 -continuous finite elements [42,43].
In this article, we present a variationally consistent framework for the direct discretization of imperfect interface models of general order. It is based on two essential components: (a) higherorder accurate smooth spline basis functions that have been used as finite element basis functions in the context of isogeometric analysis [44,45]; (b) cut-cell methods [46,47] that circumvent the meshing challenge for complex embedded interfaces; and (c) variational coupling methods [48,49] that tie together unfitted finite element meshes along embedded surfaces. This combination constitutes a natural way to tackle the key challenges stated above. Due to the higher-order smoothness of spline basis functions, derivatives of arbitrary order can be directly evaluated without affecting the consistency of the variational formulation. Unfitted meshes can accommodate complex surface geometries without meshing related obstacles, while variational coupling enables the representation of weak and strong discontinuities at embedded imperfect interfaces.
The main objective of this article is to establish the fundamental components of this concept and to demonstrate its computational efficiency for the discretization of existing higher-order imperfect interface models. In Section 2, we review the Benveniste-Miloh asymptotic model and the classical Steigmann-Ogden material surface model that involve higher-order derivatives up to fourth order. In Sections 3 and 4, we consolidate our discretization strategy based on isogeometric cut-cell finite elements and variational interface formulations in mixed and primal form. In Section 5, we demonstrate with the example of a coated circular inhomogeneity that the proposed discretization framework naturally achieves higher-order accuracy. In Section 6, we summarize the key points of our framework, including aspects that require further elaboration, and outline its potential to open up avenues for future theoretical and computational research.
In this work, we focus on two classes of higher-order imperfect interface models proposed by Benveniste and Miloh [22] and by Gurtin and Murdoch [23,24] and Steigmann and Ogden [26,27]. We briefly discuss the underlying assumptions of each model class, their interrelation, and provide specific interface conditions for three specific models.
We consider the isotropic, linearly elastic plane-strain problem shown in Fig. 1a. This problem involves an inhomogenity embedded in a matrix, separated by a uniform interphase layer, where the Lamé parameters λ and µ are given for the matrix, interphase and inhomogeneity materials, indicated by the subscripts M , Γ and I, respectively. The problem is characterized by the dimensionless length scale ε = h/R, where h is the thickness of the layer and R denotes a typical radius of curvature of the interphase midplane Γ (see in Fig. 1b). In [22], Benveniste and Miloh assumed that ε 1 and that the elastic properties of the matrix and the inhomogeneity are of the same order of magnitude. Using asymptotic analyses of the fields inside the interphase in terms of ε, they transferred the thin interphase problem of Fig. 1b into the analogous interface problem of Fig. 1c, given that appropriate conditions for the fields across the interface Γ are enforced.
The interface conditions depend on the dimensionless length scale ε and the dimensionless interphase material parameters λ Γ and µ Γ which characterize the interphase geometry and the relative elastic response of the interphase, respectively. The latter are defined as:
Figure 1: The inhomogeneity problem, zoom-in of the interphase, and the corresponding asymptotic interface limit.
Benveniste and Miloh then showed that the elastic response of the composite system changes for different interphase thicknesses. Depending on the integer N involved in the following expression, We now focus on the membrane type interface and the inextensible shell type interface that are defined by the following conditions across Γ: Membrane type interface (N = 1):
in which s is the arc length and R(s) is the radius of the curvature of the interface Γ. The subscripts t and n are the tangential and normal directions as shown in Fig. 1c. The constant P is given by:
If we assume a circular inhomogeneity as depicted in Fig. 1a with a constant radius of curvature R and a polar coordinate system {r, θ} attached to the center of the inhomogeneity, the jump conditions (3) for the tractions simplify to:
where
, and the subscript "," indicates differentiation, e.g. u θ,θ = ∂u θ /∂θ.
Inextensible shell type interface (N = 3):
in which
where D is the flexural rigidity. For a circular inhomogeneity with a polar coordinate system, the jump conditions (6) for the tractions simplify to:
supplemented by the following inextensibility condition:
It can be seen from ( 3), ( 5), ( 6), ( 8) and ( 9) that the interface conditions are governed by higherorder differential operators (second-order for N = 1 and fourth-order for N = 3). Therefore, higher-order smoothness is required for the functions that represent the displacements.
Another class of interface models are based on the material surface theories proposed by Gurtin and Murdoch [23,24] and Steigmann and Ogden [26,27]. In both theories, the interface is a material surface that has zero thickness. It is treated as a membrane in the Gurtin-Murdoch (G-M) theory and as a shell in the Steigmann-Ogden (S-O) theory. The material surface is characterized by the surface tension σ 0 , the Lamé parameters µ 0 and λ 0 (G-M and S-O), and the bending parameters ζ 0 and χ 0 (S-O). All interface parameters have dimensions of force per length in contrast to the elastic bulk parameters whose dimensions are force per area.
The equilibrium equations and boundary conditions are presented in component form in [31] and [32,33] for the Gurtin-Murdoch theory and the Steigmann-Ogden theory, respectively. They include the standard Navier equations of linear elasticity for the bulk materials and a set of equations that describe the conditions at the material surface. In both theories, the displacements are continuous across the surface and the tractions undergo jumps that are described by the surface equilibrium conditions.
We again assume a circular inhomogeneity as shown in Fig. 1a with a polar coordinate system {r, θ} attached to its center. The Steigmann-Ogden model prescribes the following set of conditions across the interface between the inhomogeneity and the matrix:
The conditions (11) reduce to those for the Gurtin-Murdoch theory if the bending parameters are neglected, i.e. χ 0 = ζ 0 = 0. For a circular inhomogeneity and for the case of vanishing surface tension, the material surface theories of Gurtin-Murdoch and Steigmann-Ogden can be related to the corresponding regimes by 50]. For example, the conditions of the membrane type interface correspond to the interface conditions of the Gurtin-Murdoch model with σ 0 = 0, if
The conditions for the inextensible shell type interface correspond to the interface conditions of the Steigmann-Ogden model with σ 0 = 0, if: i) the parameters µ 0 , λ 0 satisfy (12), ii) the parameters µ Γ , λ Γ are such that the inextensibility condition (10) is satisfied, and iii) the parameters χ 0 , ζ 0 satisfy the following equation:
In the next step, we consolidate our discretization strategy that is based on isogeometric cut-cell finite elements [51,52]. It uses higher-order continuous splines as basis functions, defined over patches of structured elements that are cut by a boundary or interface. In this section, we outline the essential technical aspects of B-splines and cut-cell quadrature methods, while variational interface formulations will be described in Section 4.
Isogeometric analysis [44,45] is an isoparametric finite element method that uses smooth and higher-order spline basis functions, ubiquitous in computer aided geometric design (CAD) to represent geometric objects, for the approximation of physics-based solutions fields. While its original objective has been a better integration of CAD and finite element analysis, it offers significant additional benefits. Approximations of derivative fields are smooth and their degree can be adjusted to what is required by the primal variational formulation [53][54][55][56]. In this work, we exploit the higher-order smoothness of splines for discretizing imperfect interface models that involve higher-order derivatives. In the following, we outline the construction of smooth B-spline basis functions [57].
A B-spline basis of degree p is formed from a sequence of knots called a knot vector Ξ = {ξ 1 , ξ 2 , . . . , ξ n+p+1 }, where ξ 1 ≤ ξ 2 ≤ . . . ≤ ξ n+p+1 and ξ ∈ R is called a knot. A repeated knot in Ξ is said to have multiplicity k. In this case, the smoothness of the B-spline basis is C p-k at that location. A univariate B-spline basis function N i,p (ξ) is defined using a recurrence relation, starting with piecewise constants (p = 0),
For p > 0, basis functions are defined using the Cox-de Boor recursion formula,
Figure 2 illustrates a B-spline basis of polynomial degree p = 3, where knots at the beginning and the end are repeated to make the basis interpolatory. We observe that the parameter space over which B-spline basis functions are defined can be split into knot spans Ωe that we denote as elements in the following.
B-spline basis functions satisfy a number of beneficial properties: They constitute a partition of unity, that is, ∀ξ, i N i,p (ξ) = 1, the support of each B-spline is compact and contained in the interval [ξ i , ξ i+p+1 ], and each B-spline function is point-wise non-negative over the entire domain. For structured spline meshes, global refinement of the spline basis functions is simple. One can either increase the number of elements (h-refinement) by inserting more knots in the knot vector or increase the polynomial degree (p-refinement) by elevating the order of the basis functions. Figures 3a and 3b illustrate hand p-refinement of the original cubic B-spline patch, where the latter maintains the maximum C p-1 -continuity (also sometimes denoted as k-refinement).
The Cox-de Boor formula (15) can be used for point-wise evaluation of the basis function values. It leads to the construction of a small, lower triangular table of intermediate values. First order derivatives and higher-order derivatives are also computed with the help of this reccursion. In particular, these derivatives are computed as a bi-product of this triangular table computation. Let us note that the evaluation is both efficient (the asymptotic complexity is the same as computing the values) and numerically stable (due to the use of the B-spline recursion).
Multivariate B-splines are a tensor-product generalization of univariate B-splines. We use d s and d p to denote the dimension of the physical and parameter spaces, respectively. Multivariate B-spline basis functions are generated from d p univariate knot vectors
where = 1, . . . , d p , p indicates the polynomial degree along parametric direction , and n is the associated number of basis functions. The resulting univariate B-spline basis functions in each direction can then be denoted by N i ,p , from which multivariate basis functions B i,p (ξ) can be constructed as Moreover, for a given direction w ∈ R d , w = 1 the directional derivative of order k is given in terms of partial derivatives
where k k 1 ,...,k dp are the multinomial coefficients and the bold multiindex which appears is applied coordinate-wise. This is a multivariate polynomial in variables w, which we evaluate with Horner's scheme, for ensuring numerical stability. We note that it is straightforward to map directional derivatives from the parametric to the physical space, since we operate on structured spline meshes with affine mappings.
Cut-cell finite element methods approximate the solution of boundary value problems using non-boundary-fitted meshes [46,47]. Their primary goal is to alleviate meshing related obstacles that often appear for geometrically complex domains. In this work, the flexibility of cut-cell methods with respect to complex geometries is important, since generating boundary-fitted spline discretizations of inhomogeneity problems outside of a CAD environment is not a trivial task. First, cut-cell methods need to be able to evaluate surface and volume integrals in cut elements [58][59][60][61][62], which we will describe in the following. Second, they need to be able to enforce interface conditions at embedded surfaces variationally [48,[63][64][65][66], which will be the target of the next section. A further important challenge outside of the scope of this work revolves around the generation of well-conditioned equation systems despite the presence of small cuts [64,67,68].
In this work, we integrate cut elements via the approach of composed Gauss quadrature, based on a hierarchical decomposition of the original element into quadrature sub-cells [69]. Figure 5 illustrates the concept for two element patches that represent a circular inhomogeneity domain and a surrounding matrix domain. First, we remove all elements and corresponding basis functions from each patch that have no support in the physical domain. We then subdivide each cut element into sub-cells that are constructed in the sense of a quadtree. Adaptive decomposition is repeated until a predefined maximum depth is reached. Each sub-cell is equipped with quadrature points of the standard Gauss rule. We note that in addition to standard Gauss point weights, a geometric weight is required that takes into account the quadtree level of each point. All quadrature points that are located outside the physical domain are discarded.
To clearly distinguish between finite elements and quadrature sub-cells in Fig. 5, finite elements are plotted in black lines, while quadrature sub-cells are plotted in blue lines. Basis functions for the approximation of the solution fields are still defined exclusively on the black mesh. Composed Gauss quadrature is robust, automatic and easy to implement, therefore fitting the purpose of the this work. For quadrature in three-dimensional cut elements, composed quadrature involves a prohibitively large number of quadrature points. In 3D, alternative quadrature schemes should therefore be considered, for instance based on geometrically accurate reparametrization of cut elements [59][60][61][62] or based on the moment fitting scheme [58,70,71].
The numerical integration of surface integrals over unfitted surfaces can be achieved independently of the underlying unfitted spline mesh. One option that maintains full geometric accuracy is to describe the surface via an analytical parametric model, if available, or via a parametric spline surface given by a CAD model. Another option is to introduce a sufficiently fine triangulation of each embedded boundary, where each triangle serves as the basis for generating independent quadrature points based on standard monomial rules. In this work, we focus on circular inhomogeneities whose boundary can be described exactly via a parametric model. Integration points, e.g. based on Gaussian quadrature, can be introduced by splitting the parametric space into suitable quadrature elements. They can then be mapped forward to the physical space using the parametric surface mapping, where each quadrature point can be related to a specific element via its global coordinate. It is important to note that quadrature points must be sufficiently dense such that each cut element contains surface quadrature points, including small cut cases.
In this section, we complement the isogeometric cut-cell finite element technology by variational formulations that enforce the displacement and traction conditions of an imperfect interface model of our choice weakly along the embedded surface. We discuss three different formulations based on a mixed form, the primal form, and a hybridized primal form. In this context, it is important to note that the cut-cell framework discussed above automatically provides a finite element space that is discontinuous across the embedded surface. As a result, both the equilibrium trans- mission conditions (tractions) and the compatibility transmission conditions (displacements) can be accurately represented in a variationally consistent way.
We start with the strong mixed form of the linear elasticity problem [72]. Under the assumption of small strains, it is expressed by the following system of first-order partial differential equations and associated boundary conditions:
with the displacement vector u, the stress tensor σ, the fourth order material stiffness tensor C, and the standard linearized strain tensor ε. The domains are illustrated in Fig. 1c. The functions u D and Γ N are prescribed boundary data, given for the Dirichlet and Neumann parts of the exterior boundary, ∂Ω D and ∂Ω N . At the imperfect interface Γ, the normal vector points from the inhomogeneity domain to the matrix domain, as indicated in Fig. 1c. Both transmission conditions are written in jump form, where the jump operator is defined as
Each of the interface conditions for the interphase models discussed in Section 2 implies a zero displacement jump, see for instance (19). The particular choice of interphase model, such as (3), ( 6) or (11), is therefore specified via the expression for the traction jump condition, h(u).
We now define a space U that contains all vector functions with square integrable gradients, and a space Σ that contains all symmetric second order tensor functions with square integrable divergence. We obtain a mixed weak formulation by multiplying both equations with a suitable test function τ ∈ Σ and integrating over the domain. Due to the major symmetry of C and the symmetry of ε(u), we can write the first equation as:
We use the following integration by parts formula in both equations, where T is some symmetric second order tensor:
This leads to:
where we have made use of the following identity to rewrite the interface terms:
where { {•} } denotes the average operator defined as { {v} } = 1 2 (v M + v I ). In the following, we propose three different approaches for incorporating the transmission and boundary conditions leading to three different finite element formulations.
The weak formulation (22) provides direct access to the displacement and traction jumps on the interface and the tractions and displacements on the domain boundary. A straightforward way to arrive at a finite element formulation is thus the substitution of the transmission conditions from (19) and the discretization of the function spaces U and Σ. We also add a consistent penalty term to the second equation to ensure stability of the final formulation, which leads to:
Our isogeometric cut-cell technology enables us to discretize the different displacement and stress components with spline basis functions. Therefore, all higher-order derivatives are welldefined in a variationally consistent sense due to the higher-order smoothness of spline basis functions. We simplify the implementation by using Voigt notation for the representation of the stress and strain tensors, which also naturally incorporates the symmetry requirement. Furthermore, we use { {u h } } in our evaluation of the traction jump condition h(u). This becomes viable in an isogeometric framework, where u h and { {u h } } have a sufficient level of continuity due to the higher-order smoothness of spline basis functions. We note that we assume sufficient smoothness of the interface Γ itself as well.
To reduce the number of degrees of freedom, we can eliminate the auxiliary variable σ h . To this end, we perform reverse integration by parts to the first line of (24), which results in:
As σ h is eliminated, the associated space Σ h is arbitrary. We can choose this space such that ∀ v ∈ U h ∃ τ ∈ Σ h : ε(v) = τ . By identifying the test function τ = ε(v) and by subtracting (25) from the second line in (24), we obtain:
where we enforce the Dirichlet boundary conditions strongly in the space U .
We replace the remaining { {σ h } } on the interface by { {C : ε(u h )} }. The two quantities are generally not identical due to the discrete nature of (24). Therefore, this replacement is a (variationally consistent) approximation that simplifies (26). As a result, the finite element formulation ( 26) is not equivalent to (24).
The stabilization term that we added to ( 24) and ( 26) is intended to improve the stability characteristics of the discretized variational formulations. According to our numerical experiments, choosing the corresponding free parameter α between 1 and 100 supports accurate solutions and monotone convergence. We would like to highlight, however, that the stabilization term and its free parameter α cannot guarantee the coercivity of (26) in general. This is due to the form of the function h({ {u h } }). Since the term that involves this function is nonsymmetric, it can only be bound from below as a negative term. Since h({ {u h } }) involves averages of u h , it cannot be generally balanced by a penalty-like contraint on the jump of u h (i.e., for zero jump the function h({ {u h } }) would still have an adverse effect on the coercivity estimate). Therefore, finding an effective and generally robust stabilization of ( 24) and (26) is still an open question. In our numerical experiments, however, we were not able to find any unstable behavior for any of the cases we looked at.
In this section, we present three numerical tests that demonstrate the accuracy of the isogeometric cut-cell framework in conjunction with the variational interface formulations for the consistent discretization of higher-order imperfect interface models. From a geometric viewpoint, we restrict ourselves to the classical plane-strain configuration of a coated circular inhomogeneity embedded into an infinite matrix. This problem is of particular interest, as there exists analytical solutions for asymptotic interphase models and material surface theories. Within this configuration, we consider three different interface models (membrane type, inextensible shell type and Gurtin-Murdoch models) and two different variational interface formulations (primal and mixed forms).
Our framework is implemented in the C++ library G+Smofoot_0 , a general purpose, open-source library especially designed for isogeometric analysis. It features tensor-product and (truncated) hierarchical splines of arbitrary degree and dimension as well as a flexible syntax for implementing elaborate variational formulations (see e.g. [73]).
We consider a three-phase composite system that involves a circular inhomogeneity, co-centrical thin interphase layer, and an infinite matrix (the configuration and the following notation follow Fig. 1). The dimensionless shear modulus of the inhomogeneity is set to µ I /µ M = 0.5, and the dimensionless Lamé parameters λ M /µ M and λ I /µ M are 1.5 and 0.5, respectively. The dimensionless elastic properties of the thin interphase layer are µ Γ /µ M = 262.5 and λ Γ /µ M = 393.75. According to the classification of Benveniste and Miloh, this layer can be represented by the zero-thickness membrane type interface (5) with P/(µ M R) = 7.5. The problem is subject to a normalized far-field horizontal stress σ xx /µ M = 2. To obtain a numerical solution, we consider the finite domain [-2R, 2R] × [-2R, 2R], in which the inhomogeneity is centrally located, and impose Dirichlet boundary conditions at the outer boundaries taken from the analytical solution in [31,33]. The composite configuration is discretized with the isogeometric cut-cell framework as illustrated in Fig. 5. We first consider a base mesh with 20 × 20 cubic elements for the complete square domain. The base mesh is then used to extract unfitted meshes for the matrix and the inhomogeneity by removing all elements and associated basis functions that are completely outside the physical domain. We recall that all cut elements carry basis functions that involve independent degrees of freedom for the matrix domain and the inhomogeneity domain. Each cut element is equipped with five levels of quadrature sub-cells. A representative example discretization of the classical circular inhomogeneity problem is given in Fig. 5.
We first illustrate the accuracy of our method by juxtaposing analytical and numerical fields. Figures 6 and7 show the distribution of the dimensionless displacements u x /R and the dimensionless von Mises stresses σ/µ M . The numerical results are obtained with the primal finite element formulation, discretized with 20 × 20 cubic elements. We observe that the displacements are continuous, while the von Mises stress undergoes jumps, as dictated by the imperfect interface model of equation (5). The close agreement between the analytical and numerical solutions is a first indication of the higher-order accuracy of the proposed framework.
Figure 8 plots the traction jumps [[σ rr ]] along the circular interface, enabled by the discontinuous cut-cell solution across the embedded interface. Results were obtained with the mixed and primal finite element formulations discussed in Section 4, discretized with 20 × 20 quadratic and quintic elements. We see that the results of the quadratic discretization exhibits oscillations along the interface at the given mesh size, where no significant difference can be observed between the results of the mixed and primal formulations. When we increase the polynomial degree, the oscillations quickly vanish, indicating a clear advantage of higher-order discretizations for these models.
In the next step, we consider a series of base meshes with 5 × 5, 7 × 7, 10 × 10, 15 × 15, and 25 × 25 higher-order elements and polynomial degrees from p = 2 to p = 5. We quantify the performance of the mixed and primal finite element formulations in terms of the displacement and stress errors in the L 2 norm, when the mesh is refined at a constant polynomial degree. In addition to the error defined over the matrix and inhomogeneity domains, we also compute a corresponding error that is defined only at the interface in terms of the displacements and the stress jumps. Figures 9 through 12 plot the resulting convergence curves for for each polynomial degree. We observe that we consistently achieve optimal higher-order rates of convergence in all cases. The error converges with a comparable rate in the bulk domains and at the interface. From the stress plots at the interface and the convergence curves, we can conclude that the mixed and primal formulations provide equivalent solution accuracy. The primal formulation, however, requires significantly fewer degrees of freedom to achieve the same level of accuracy. We therefore consider only the primal formulation in the following examples.
In all computations, we use a stabilization parameter α = 100. Following Remark 1 at the end of Section 4, however, the question of appropriate stabilization is still unclear at this point. We note that some of the curves indicate superconvergent behavior for the first refinement steps, due to the comparatively large error for the coarsest meshes. One could presume that given a more effective stabilization, the large errors in the coarse discretizations and hence the superconvergent behavior of the convergence curve could also be reduced.
In the next step, we consider a significantly stiffer interphase layer with elastic properties µ Γ /µ M = 2, 625, 000 and λ Γ /µ M = 3, 937, 500, but otherwise assume the same three-phase composite system and the same material parameters for the inhomogeneity and the matrix as in Section 5.1. Following (9), the interphase can be modeled as a zero-thickness inextensible shell type interface, with P/(µ M R) = 75, 000 and D/(µ M R 3 ) = 0.625. The system is subjected to a dimensionless far-field loading σ xx /µ M = 2, for which an analytical solution exists [32,33].
To find a numerical solution, we again consider the square domain [-2R, 2R] × [-2R, 2R], where we find Dirichlet boundary conditions from the analytical solution in [32,33]. We employ the same cut-cell discretization procedure as in Section 5.1, but focus on the primal formulation only. Figures 13 and14 plot the dimensionless horizontal displacement u x /R and the dimensionless von Mises stress σ/µ M , respectively, where the numerical solutions are obtained with a 20 × 20 mesh of quintic elements. Figure 15 plots the corresponding normal traction jump [[σ rr ]] and tangential traction jump [[σ rθ ]] along the circular interface. We observe that the proposed cutcell framework yields accurate solutions in all cases, with the numerical and analytical solutions being indistinguishable from each other.
We again consider a series of base meshes with 5 × 5, 7 × 7, 10 × 10, 15 × 15, and 25 × 25 elements, with a fixed polynomial degree p = 5. Figures 16 and17 plot the convergence of the displacement and stress errors in the L 2 norm, when the mesh is uniformly refined. In Fig. 16, the L 2 errors are computed over the matrix and inhomogeneity domains. In Fig. 17, the L 2 errors are computed for the displacements and stress jumps at the interface only. We observe the same characteristic solution behavior as in Section 5.1. The proposed framework achieves higher-order optimal rates of convergence in the bulk and at the interface. We see superconvergence in the preasymptotic range that is likely due to larger errors for the coarsest meshes.
In the next step, we would like to briefly illustrate that our framework also covers the analysis of interfaces based on material surface theories. To this end, we consider the same composite configuration, but at the nanoscale, such that the Gurtin-Murdoch conditions represent the interface behavior between the matrix and the circular inhomogeneity. The Gurtin-Murdoch conditions follow from the Steigmann-Odgen model (11) with χ 0 = ζ 0 = 0. We note that the following material parameters refer to nano-composite materials made from anodic alumina with a radius of inhomogeneities within 5 and 100 nm. The data related to the surface properties of alumina and the magnitude of the surface tension are taken from [74] and [75], respectively. The elastic properties of the inhomogeneity and the matrix are µ I = 17.35 GPa, λ I = 26.025 GPa and µ M = 34.7 GPa, λ M = 52.05 GPa, respectively. The radius of the circular inhomogeities is R = 1 nm. The elastic properties of the interface are µ 0 = -6.2178 N/m and λ 0 = 3.48912 N/m and the surface tension in the interface is assumed to be σ 0 = 0.72 N/m. The composite system is subjected to a far-field horizontal loading of σ xx = 10 GPa, for which an analytical solution exists [31,33].
For our numerical solution, we consider the square domain [-2R, 2R] × [-2R, 2R], finding Dirichlet boundary conditions from the analytical solution in [31,33]. We employ the same cutcell discretization procedure as in Section 5.1, but focus on the primal formulation only. We discretize the problem with a base mesh of 20 × 20 cubic elements. Figures 18 and19 plot the von Mises stress over the complete computational domain and the corresponding normal and tangential traction jumps [[σ rr ]] and [[σ rθ ]] along the circular interface. We observe that the proposed cut-cell framework yields higher-order accurate solutions with respect to the analytical reference in all cases. Figures 20a and20b plot the convergence of the stress error in the L 2 norm for the matrix and inhomogeneity domains and for the interface, when the mesh is uniformly refined. The results confirm higher-order optimal rates of convergence.
Although in the scope of this work, we do not include geometrically complex cases, we would like to provide a flavor of the geometric flexibility of the cut-cell framework as one of its core advantages over standard boundary-fitted methods. To this end, we image that we need to extend the current problem to include two circular inhomogeneities whose radius can be arbitrarily changed and whose position can arbitrarily moved within the computational domain. In particular, this can involve a topology change, when the inhomogeneities are starting to overlap, transitioning from two distinct to one single inhomogeneity domain. It is straightforward to see that the proposed framework can handle all cases without any further adjustments, whereas standard boundary-fitted cases require a potentially cumbersome mesh generation step for each configurational change. In particular, the cut-cell approach enables the accurate representation of the strong discontinuities in the stress at arbitrary locations within cut elements.
Asymptotic interphase models replace volumes of thin layers by dimensionally reduced surfaces, along which appropriate jump conditions of the fields are formulated. They are closely related to theories of material surfaces that incorporate the effect of surface tension in the analysis of nanocomposites. Both approaches formulate interface conditions with higher-order derivatives that lead to solutions with weak and strong discontinuities at the interface, and are therefore examples of the general class of higher-order "imperfect" interface models. Their numerical treatment constitutes a challenge, as standard C 0 -continuous finite element methods are unable to represent strong discontinuities and do not allow higher-order derivatives.
In this article, we presented first steps towards a variationally consistent framework that enables the consistent and direct discretization of imperfect interface models of general order. Its basic idea is the synergistic combination of three components: (a) smooth spline basis functions to directly evaluate higher-order derivatives of arbitrary order without affecting the variational consistency, (b) cut-cell finite element technology to accommodate complex interface geometries without meshing related obstacles, (c) variational coupling to weakly enforce discontinuities dictated by imperfect interface models. We demonstrated the computational efficiency of this idea for the classical benchmark of a coated circular inhomogeneity. We employed three higher-order imperfect interface models, namely the membrane type and inextensible shell type models of the Benveniste-Miloh classification and the classical Gurtin-Murdoch material surface model. We showed for each test case that the proposed framework achieves higher-order accuracy in terms of solution fields and error convergence, evaluated across the bulk domains and directly at the interface.
The results presented in this article open up several directions for future work. On the one hand, we identified several open questions in the context of the proposed framework. The most important one is the derivation of better stabilization mechanisms that are able to generally ensure the robustness of the proposed framework. Another open question is the treatment of rough interface surfaces that include deviations from the perfect geometry, kinks, or holes with free ends. This is particularly interesting for the understanding of real coatings that always contain material or geometric imperfections. On the other hand, the framework is open to be explored for more complex test cases. This includes of course its extension to three dimensions, but also its integration in computational homogenization procedures, where it enables the computation of effective macroscale properties when higher-order imperfect interface models are used at the microscale [76]. It is also desirable to better understand the modeling capabilities of imperfect interfaces, for instance by comparing their accuracy with interphase models based on embedded membranes or shells [77,78]. It is our hope that in a longer-term perspective the proposed framework will help inspire the development of new higher-order imperfect interface models. performed as part of a joint PhD training program. Sofia Mogilevskaya gratefully acknowledges support provided by the Theodore W. Bennett Chair at the University of Minnesota. D. Schillinger gratefully acknowledges support from the National Science Foundation through the NSF CAREER Award No. 1651577 and the NSF grant CISE-156599. The authors also acknowledge the Minnesota Supercomputing Institute (MSI) of the University of Minnesota for providing computing resources that have contributed to the research results reported within this paper (https://www.msi.umn.edu/).
Geometry plus Simulation Modules, http://github.com/gismo