Hybrid coupling of CG and HDG discretizations based on Nitsche's method

TL;DR

This paper introduces a minimally-intrusive hybrid coupling strategy for CG and HDG discretizations using Nitsche's method, enabling effective simulation of thermal and elastic problems with multiple materials.

Contribution

It presents a novel coupling approach based solely on the HDG hybrid variable, preserving the core structures of CG and HDG methods and facilitating integration into existing libraries.

Findings

Achieves optimal convergence of stress fields.

Demonstrates superconvergence of displacement.

Handles multiple materials with different compressibility.

Abstract

A strategy to couple continuous Galerkin (CG) and hybridizable discontinuous Galerkin (HDG) discretizations based only on the HDG hybrid variable is presented for linear thermal and elastic problems. The hybrid CG-HDG coupling exploits the definition of the numerical flux and the trace of the solution on the mesh faces to impose the transmission conditions between the CG and HDG subdomains. The continuity of the solution is imposed in the CG problem via Nitsche's method, whereas the equilibrium of the flux at the interface is naturally enforced as a Neumann condition in the HDG global problem. The proposed strategy does not affect the core structure of CG and HDG discretizations. In fact, the resulting formulation leads to a minimally-intrusive coupling, suitable to be integrated in existing CG and HDG libraries. Numerical experiments in two and three dimensions show optimal global convergence of the stress and superconvergence of the displacement field, locking-free approximation, as well as the potential to treat structural problems of engineering interest featuring multiple materials with compressible and nearly incompressible behaviors.

Figures

Equations

$$\{\begin{array}{rlrl}-\text{boldmath∇unboldmath}\cdot \text{boldmath∇unboldmath}u& =f& & \text{in Ω,}\\ u& =u_D& & \text{on ΓD,}\\ \text{boldmathnunboldmath}\cdot \text{boldmath∇unboldmath}u& =t& & \text{on ΓN,}\end{array}$$

$$\Omega =\bigcup _{e=1}^{{\text{n}}_{\text{el}}}{\Omega }_e,{\Omega }_i\cap {\Omega }_j=\varnothing \text{ for }i\ne j.$$

$${\mathcal{V}}^h(\Omega )$$

$${\mathcal{V}}_0^h(\Omega )$$

$$(u,v{)}_D:={\int }_{D}uvd\Omega ,⟨u,v{⟩}_S:={\int }_{S}uvd\Gamma .$$

$$(\text{boldmath∇unboldmath}v,\text{boldmath∇unboldmath}{u}^h{)}_{\Omega }=(v,f{)}_{\Omega }+⟨v,t{⟩}_{{\Gamma }_N}$$

$$\Gamma :=\left[\bigcup _{e=1}^{{\text{n}}_{\text{el}}}\partial {\Omega }_e\right]\setminus \partial \Omega .$$

$$\{\begin{array}{rlrl}\text{boldmathqunboldmath}+\text{boldmath∇unboldmath}u& =\text{boldmath0unboldmath}& & \text{in Ωe,e=1,…,nel,}\\ \text{boldmath∇unboldmath}\cdot \text{boldmathqunboldmath}& =f& & \text{in Ωe,e=1,…,nel,}\\ u& =u_D& & \text{on ∂Ωe∩ΓD,}\\ u& =\hat{u}& & \text{on ∂Ωe∖ΓD,}\\ \text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}& =-t& & \text{on ∂Ωe∩ΓN,}\\ [[u\text{boldmathnunboldmath}]]& =\text{boldmath0unboldmath}& & \text{on Γ,}\\ [[\text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}]]& =0& & \text{on Γ,}\end{array}$$

$$[[\odot ]]={\odot }_l+{\odot }_r.$$

$$\{\begin{array}{rlrl}{\text{boldmathqunboldmath}}_e+\text{boldmath∇unboldmath}{u}_e& =\text{boldmath0unboldmath}& & \text{in Ωe,e=1,…,nel,}\\ \text{boldmath∇unboldmath}\cdot {\text{boldmathqunboldmath}}_e& =f& & \text{in Ωe,e=1,…,nel,}\\ u_e& =u_D& & \text{on ∂Ωe∩ΓD,}\\ u_e& =\hat{u}& & \text{on ∂Ωe∖ΓD.}\end{array}$$

$$\{\begin{array}{rlrl}[[u\text{boldmathnunboldmath}]]& =\text{boldmath0unboldmath}& & \text{on Γ,}\\ [[\text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}]]& =0& & \text{on Γ,}\\ \text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}& =-t& & \text{on ΓN,}\end{array}$$

$${\mathcal{W}}^h(\Omega )$$

$$\displaystyle\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{$\mathcal{W}$}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(S$$

$$\text{boldmathnunboldmath}\cdot \text{stackon}[0pt]\text{boldmathqunboldmath}{\text{stretchto}\text{scalerel}\ast [\text{widthof}\text{boldmathqunboldmath}]\bigwedge {0.5ex}_{\text{belowbaseline}[-2pt]}}^h:=\{\begin{array}{ll}\text{boldmathnunboldmath}\cdot {\text{boldmathqunboldmath}}_e^h+\tau (u_e^h-u_D)& \text{on ∂Ωe∩ΓD,}\\ \text{boldmathnunboldmath}\cdot {\text{boldmathqunboldmath}}_e^h+\tau (u_e^h-{\hat{u}}^h)& \text{elsewhere,}\end{array}$$

$$-(\text{boldmathwunboldmath},{\text{boldmathqunboldmath}}_e^h{)}_{{\Omega }_e}+(\text{boldmath∇unboldmath}\cdot \text{boldmathwunboldmath},u_e^h{)}_{{\Omega }_e}$$

$$(v,\text{boldmath∇unboldmath}\cdot {\text{boldmathqunboldmath}}_e^h{)}_{{\Omega }_e}+⟨v,\tau u_e^h{⟩}_{\partial {\Omega }_e}$$

$$\sum_{e=1}^{\texttt{n}_{\texttt{el}}}\Bigl{\{}\langle\hat{v},\text{\boldmath$n$\unboldmath}{\cdot}\text{\boldmath$q$\unboldmath}_{e}^{h}\rangle_{\partial\Omega_{e}\setminus\Gamma_{D}}{+}\langle\hat{v},\tau u_{e}^{h}\rangle_{\partial\Omega_{e}\setminus\Gamma_{D}}{-}\langle\hat{v},\tau\hat{u}^{h}\rangle_{\partial\Omega_{e}\setminus\Gamma_{D}}\Bigr{\}}={-}\sum_{e=1}^{\texttt{n}_{\texttt{el}}}\langle\hat{v},t\rangle_{\partial\Omega_{e}\cap\Gamma_{N}},$$

$$\{\begin{array}{rlrl}{u}_{\text{CG}}& ={\hat{u}}_{\text{HDG}}& & \text{on ΓI,}\\ -{\text{boldmathnunboldmath}}_I\cdot \text{boldmath∇unboldmath}{u}_{\text{CG}}& ={\text{boldmathnunboldmath}}_I\cdot {\text{boldmathqunboldmath}}_{\text{HDG}}& & \text{on ΓI,}\end{array}$$

$$-{\text{boldmathnunboldmath}}_I\cdot \text{stackon}[0pt]\text{boldmath∇unboldmath}{u}_{\text{CG}}\text{stretchto}\text{scalerel}\ast [\text{widthof}\text{boldmath∇unboldmath}{u}_{\text{CG}}]\bigwedge {0.5ex}_{\text{belowbaseline}[-2pt]}:=-{\text{boldmathnunboldmath}}_I\cdot \text{boldmath∇unboldmath}{u}_{\text{CG}}+\frac{\gamma }{h}(u_{\text{CG}}-{\hat{u}}_{\text{HDG}}),$$

$$\begin{array}{rl}(\text{boldmath∇unboldmath}v,\text{boldmath∇unboldmath}{u}_{\text{CG}}{)}_{{\Omega }_{\text{CG}}}-⟨v,\text{boldmathnunboldmath}\cdot \text{boldmath∇unboldmath}& u_{\text{CG}}{⟩}_{{\Gamma }_I}-⟨\text{boldmathnunboldmath}\cdot \text{boldmath∇unboldmath}v,u_{\text{CG}}{⟩}_{{\Gamma }_I}+⟨v,\frac{\gamma }{h}{u}_{\text{CG}}{⟩}_{{\Gamma }_I}\\ & =(v,f{)}_{{\Omega }_{\text{CG}}}+⟨v,t{⟩}_{{\Gamma }_N}-⟨\text{boldmathnunboldmath}\cdot \text{boldmath∇unboldmath}v,{\hat{u}}_{\text{HDG}}{⟩}_{{\Gamma }_I}+⟨v,\frac{\gamma }{h}{\hat{u}}_{\text{HDG}}{⟩}_{{\Gamma }_I},\end{array}$$

$$\displaystyle\begin{aligned} \sum_{e=1}^{\texttt{n}_{\texttt{el}}}\Bigl{\{}\langle\hat{v},\text{\boldmath$n$\unboldmath}{\cdot}{\text{\boldmath$q$\unboldmath}}_{\texttt{HDG}}\rangle_{\partial\Omega_{e}\setminus\Gamma_{D}}{+}&\langle\hat{v},\tau{u}_{\texttt{HDG}}\rangle_{\partial\Omega_{e}\setminus\Gamma_{D}}{-}\langle\hat{v},\tau{\hat{u}}_{\texttt{HDG}}\rangle_{\partial\Omega_{e}\setminus\Gamma_{D}}{+}\langle\hat{v},\tfrac{\gamma}{h}{\hat{u}}_{\texttt{HDG}}\rangle_{\partial\Omega_{e}\cap\Gamma_{I}}\Bigr{\}}\\ &{=}\sum_{e=1}^{\texttt{n}_{\texttt{el}}}\Bigl{\{}{-}\langle\hat{v},t\rangle_{\partial\Omega_{e}\cap\Gamma_{N}}{-}\langle\hat{v},\text{\boldmath$n$\unboldmath}{\cdot}\text{\boldmath$\nabla$\unboldmath}{u}_{\texttt{CG}}\rangle_{\partial\Omega_{e}\cap\Gamma_{I}}{+}\langle\hat{v},\tfrac{\gamma}{h}{u}_{\texttt{CG}}\rangle_{\partial\Omega_{e}\cap\Gamma_{I}}\Bigr{\}},\end{aligned}$$

$$\left[\begin{array}{cc}{\mathbf{K}}_{\text{CG}}& {\mathbf{K}}_{\text{I}}\\ {\mathbf{K}}_{\text{I}}^T& {\mathbf{K}}_{\text{HDG}}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_{\text{CG}}\\ {\hat{\mathbf{u}}}_{\text{HDG}}\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}_{\text{CG}}\\ {\mathbf{f}}_{\text{HDG}}\end{array}\right],$$

$$\{\begin{array}{rlrl}-\text{boldmath∇unboldmath}\cdot \text{boldmathσunboldmath}& =\text{boldmathfunboldmath}& & \text{in Ω,}\\ \text{boldmathσunboldmath}& ={\text{boldmathσunboldmath}}^T,& & \text{in Ω,}\\ \text{boldmathuunboldmath}& ={\text{boldmathuunboldmath}}_D& & \text{on ΓD,}\\ \text{boldmathnunboldmath}\cdot \text{boldmathσunboldmath}& =\text{boldmathtunboldmath}& & \text{on ΓN,}\end{array}$$

$$\text{boldmathσunboldmath}=\frac{E}{1+\nu }\text{boldmathεunboldmath}(\text{boldmathuunboldmath})+\frac{E\nu }{(1+\nu )(1-2\nu )}\mathrm{tr}(\text{boldmathεunboldmath}(\text{boldmathuunboldmath})){\mathbf{I}}_{{\text{n}}_{\text{sd}}},$$

$$\text{\boldmath$\sigma$\unboldmath}_{\texttt{V}}:=\begin{cases}\bigl{[}\sigma_{xx},\;\sigma_{yy},\;\tau_{xy}\bigr{]}^{T}&\text{in 2D,}\\ \bigl{[}\sigma_{xx},\;\sigma_{yy},\;\sigma_{zz},\;\tau_{xy},\;\tau_{xz},\;\tau_{yz}\bigr{]}^{T}&\text{in 3D,}\end{cases}$$

$$\text{\boldmath$\varepsilon$\unboldmath}_{\texttt{V}}:=\begin{cases}\bigl{[}\varepsilon_{xx},\;\varepsilon_{yy},\;\gamma_{xy}\bigr{]}^{T}&\text{in 2D,}\\ \bigl{[}\varepsilon_{xx},\;\varepsilon_{yy},\;\varepsilon_{zz},\;\gamma_{xy},\;\gamma_{xz},\;\gamma_{yz}\bigr{]}^{T}&\text{in 3D.}\end{cases}$$

$${\text{boldmath∇unboldmath}}_{\text{S}}:=\{\begin{array}{ll}{\left[\begin{array}{ccc}\frac{\partial }{\partial x}& 0& \frac{\partial }{\partial y}\\ 0& \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]}^T& \text{in 2D,}\\ {\left[\begin{array}{cccccc}\frac{\partial }{\partial x}& 0& 0& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}& 0\\ 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial x}& 0& \frac{\partial }{\partial z}\\ 0& 0& \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial x}& \frac{\partial }{\partial y}\end{array}\right]}^T& \text{in 3D,}\end{array}$$

$$\mathbf{D}:=\{\begin{array}{ll}\lambda \left[\begin{array}{ccc}1+(1-\vartheta )\nu & \nu & 0\\ \nu & 1+(1-\vartheta )\nu & 0\\ 0& 0& \frac{1-\vartheta \nu }{2}\end{array}\right]& \text{in 2D,}\\ \lambda \left[\begin{array}{cccc}1-\nu & \nu & \nu & \\ \nu & 1-\nu & \nu & {\text{boldmath0unboldmath}}_{{\text{n}}_{\text{sd}}}\\ \nu & \nu & 1-\nu & \\ & {\text{boldmath0unboldmath}}_{{\text{n}}_{\text{sd}}}& & \frac{1-2\nu }{2}{\mathbf{I}}_{{\text{n}}_{\text{sd}}}\end{array}\right]& \text{in 3D,}\end{array}$$

$$\lambda :=\{\begin{array}{ll}\frac{E}{(1+\nu )(1-\vartheta \nu )}& \text{in 2D,}\\ \frac{E}{(1+\nu )(1-2\nu )}& \text{in 3D,}\end{array}$$

$$\mathbf{N}:=\{\begin{array}{ll}{\left[\begin{array}{ccc}{n}_x& 0& n_y\\ 0& n_y& n_x\end{array}\right]}^T& \text{in 2D,}\\ {\left[\begin{array}{cccccc}{n}_x& 0& 0& n_y& n_z& 0\\ 0& n_y& 0& n_x& 0& n_z\\ 0& 0& n_z& 0& n_x& n_y\end{array}\right]}^T& \text{in 3D.}\end{array}$$

$$\{\begin{array}{rlrl}-{\text{boldmath∇unboldmath}}_{\text{S}}^T{\text{boldmathσunboldmath}}_{\text{V}}& =\text{boldmathfunboldmath}& & \text{in Ω,}\\ {\text{boldmathσunboldmath}}_{\text{V}}& =\mathbf{D}{\text{boldmathεunboldmath}}_{\text{V}}& & \text{in Ω,}\\ \text{boldmathuunboldmath}& ={\text{boldmathuunboldmath}}_D& & \text{on ΓD,}\\ {\mathbf{N}}^T{\text{boldmathσunboldmath}}_{\text{V}}& =\text{boldmathtunboldmath}& & \text{on ΓN.}\end{array}$$

$$\{\begin{array}{rlrl}-{\text{boldmath∇unboldmath}}_{\text{S}}^T\mathbf{D}{\text{boldmath∇unboldmath}}_{\text{S}}{\text{boldmathuunboldmath}}_{\text{CG}}& =\text{boldmathfunboldmath}& & \text{in ΩCG,}\\ {\text{boldmathuunboldmath}}_{\text{CG}}& ={\text{boldmathuunboldmath}}_D& & \text{on ΓD,}\\ {\mathbf{N}}^T\mathbf{D}{\text{boldmath∇unboldmath}}_{\text{S}}{\text{boldmathuunboldmath}}_{\text{CG}}& =\text{boldmathtunboldmath}& & \text{on ΓN,}\end{array}$$

Full text

Keywords: Hybridizable discontinuous Galerkin, Coupling with finite element, Nitsche’s method, Locking-free, Superconvergence

1 Introduction

The finite element (FE) method has been extensively used in simulation-based engineering over the last $50$ years. Strategies to enhance FE computations by coupling them with other discretization methods have received a lot of attention both in the mathematical and engineering community. Examples of successful couplings include, e.g., FE and meshless methods [1, 2, 3, 4], FE and finite volume methods [5, 6, 7], as well as advanced discretization techniques as mixed and discontinuous Galerkin (DG) methods [8, 9].

The idea of coupling distinct numerical techniques in different regions of the computational domain may also be interpreted in the framework of domain decomposition (DD) methods. In this context, and for nonverlapping strategies, a key aspect is represented by the so-called transmission conditions to be imposed between two different regions of the domain. Such conditions enforce continuity of the solution and of the flux/stress across the interface and have been treated in the literature either via Lagrange multipliers [10] or using Nitsche’s method [11]. On the one hand, couplings involving FE discretizations and Lagrange multipliers has lead to the well-known Mortar element method [12, 13, 14, 15, 16, 17, 18, 19, 20]. On the other hand, DD methods based on DG approximations have been studied using both Lagrange multipliers [21, 22, 23] and Nitsche’s method [24, 25].

Among existing numerical methods, the coupling of CG and DG in different regions of the domain has been extensively studied in the literature to exploit the advantages of both approaches. This is of special interest in the context of multiphysics and multimaterial problems in which different regions of the computational domain feature distinct physical properties, for which specific discretizations need to be devised. On the one hand, CG leads to computationally efficient discretizations with a limited number of degrees of freedom [26]. On the other hand, DG provides a flexbile paradigm to handle meshes with hanging nodes and construct nonuniform polynomial degree and high-order approximations [27, 28, 29, 30]. Moreover, DG methods have proved to be very efficient in stabilizing convection terms in conservation laws [31, 32].

The idea of coupling CG and DG, specifically the local discontinuous Galerkin (LDG) [33], via an appropriate definition of the numerical flux was first proposed in [34] to handle nonmatching grids. The advantage of a CG-LDG coupling in presence of convection-dominated problems was later discussed in [35, 36, 37].

More recently, hybrid discretization techniques, especially HDG [38, 39, 40, 41, 42, 43] and the hybrid high-order (HHO) method [44, 45, 46, 47], have gained a lot of attention owing to their reduced computational costs with respect to classical DG approaches like LDG. In the context of HDG, coupling with the boundary element method has been discussed in [48, 49], whereas a first attempt to couple HDG and CG discretizations has been proposed in [50]. This approach requires the introduction of an appropriate projection operator to enforce the transmission condition in the HDG local problem. This results in a coupling of local and global degrees of freedom of the HDG problem with the ones of the CG discretization, making the implementation of this strategy in existing CG and HDG libraries intrusive.

This work proposes a CG-HDG coupling in which solely the HDG hybrid variable is exploited in the transmission conditions. This approach requires two ingredients. On the one hand, an appropriate definition of the trace of the numerical normal flux at the interface between the CG and HDG subdomains. On the other hand, the weak imposition of Dirichlet boundary conditions in the CG formulation via Nitsche’s method. A formulation of the latter in an optimization framework is described in [51]. The resulting hybrid CG-HDG coupling does not affect the structure of the core CG and HDG matrices, thus leading to a minimally-intrusive implementation of this technique in existing CG and HDG libraries.

This approach is of special interest in the context of linear elastic problems since displacement-based formulations fail to provide locking-free approximations in nearly incompressible materials, when low-order CG discretizations are utilized [52]. To remedy this issue, mixed [53, 54, 55, 56] and equilibrium formulations [57], as well as discretizations based on the nonconforming Crouzeix-Raviart element [58] and DG techniques [59, 60, 61] have been proposed. The above mentioned HDG method also provides locking-free approximations [62, 63, 64, 65, 66] while preserving the advantages of a DG discretization with hybridization. It relies on a mixed hybrid formulation [67] with polynomial approximations discontinuous element-by-element. More precisely, the HDG formulation for linear elasticity proposed in [68, 69] is utilized. Exploiting the well-known Voigt notation for second-order symmetric tensors [70], this approach provides optimal convergence of the stress tensor and superconvergence of the displacement field using equal order approximation for all the variables, even for low-order polynomial functions. Hence, a CG approximation in the compressible region of the domain is coupled with an HDG solver for the nearly incompressible one in order to study multimaterial problems of engineering interest.

The remaining of this paper is organized as follows. First, a linear thermal problem is considered to introduce the CG and HDG discretizations, as well as the proposed hybrid coupling (Section 2). In Section 3, the coupled CG-HDG discretization is presented for a linear elastic problem exploiting the HDG-Voigt formulation. Section 4 is devoted to the numerical validation of the hybrid coupling in two dimensions. More precisely, optimal orders of convergence, a sensitivity analysis to the parameters of the method and robustness in the incompressible limit are verified. Special emphasis is put in the development of a strategy based on nonuniform polynomial degree approximation to obtain optimal global convergence of the stress of order $k{+}1$ and superconvergence of the displacement field of order $k{+}2$. In Section 5, two and three-dimensional elastic problems involving composite materials with space-dependent mechanical properties are analyzed and Section 6 summarizes the results of this work.

2 CG-HDG coupling for thermal problems

The hybrid coupling of CG and HDG discretizations is presented for a thermal problem described by a second-order scalar elliptic partial differential equation (PDE). Consider an open bounded domain $\Omega\in\mathbb{R}^{\texttt{n}_{\texttt{sd}}}$ with boundary $\partial\Omega{=}\Gamma_{D}\cup\Gamma_{N}$ such that $\Gamma_{D}\cap\Gamma_{N}{=}\emptyset$ and $\texttt{n}_{\texttt{sd}}$ being the number of spatial dimensions. The strong form of the Poisson equation is

[TABLE]

where $u$ represents the unknown temperature, $f$ is a user-prescribed source term and $u_{D}$ and $t$ denote the Dirichlet and Neumann boundary data, respectively.

In the following subsections, the discrete forms of the CG and HDG approximations are first recalled separately and the coupling strategy based on Nitsche’s method is thus presented.

2.1 CG approximation

Assume that $\Omega$ is partitioned in $\texttt{n}_{\texttt{el}}$ disjoint subdomains $\Omega_{e}$ such that

[TABLE]

The following discrete functional spaces are introduced

[TABLE]

where $\mathcal{P}^{k}(\Omega_{e})$ is the space of polynomial functions of complete degree at most $k{\geq}1$ on each mesh element $\Omega_{e},\,e{=}1,\ldots,\texttt{n}_{\texttt{el}}$. Moreover, $(\cdot,\cdot)_{D}$ and $\langle\cdot,\cdot\rangle_{S}$ are introduced to denote the usual $\mathcal{L}_{2}$ inner products on a generic subdomain $D$ and on a generic surface $S$, namely

[TABLE]

The discrete form of the CG approximation of Equation (1) is thus obtained by integrating by parts the corresponding weak form and it reads: find $u^{h}\in\mathcal{V}^{h}(\Omega)$ such that

[TABLE]

for all $v\in\mathcal{V}^{h}_{0}(\Omega)$.

2.2 HDG approximation

Consider the partition of the domain introduced in Equation (2) and define the internal interface as

[TABLE]

The HDG formulation of the Poisson problem under analysis is obtained by rewriting Equation (1) element-by-element on $\Omega_{e},\,e{=}1,\ldots,\texttt{n}_{\texttt{el}}$ as a system of first-order PDEs, via the introduction of the mixed variable \text{\boldmathq\unboldmath}{=}{-}\text{\boldmath\nabla\unboldmath}u and the hybrid variable $\hat{u}$ representing the trace of the solution on $\partial\Omega_{e}\setminus\Gamma_{D}$, namely

[TABLE]

where the last two equations are transmission conditions enforcing the continuity of the solution and of the normal flux across the interface $\Gamma$. The jump operator $\llbracket\cdot\rrbracket$ is defined according to [71] as the sum of the values from the elements $\Omega_{l}$ and $\Omega_{r}$, respectively on the left and on the right of a given interface, namely

[TABLE]

The solution of the HDG problem is performed in two stages [38, 39, 40, 41, 42, 43]. First, a set of $\texttt{n}_{\texttt{el}}$ local problems is defined to determine (u_{e},\text{\boldmathq\unboldmath}_{e}) as functions of the unknown hybrid variable $\hat{u}$ in each element $\Omega_{e},\,e{=}1,\ldots,\texttt{n}_{\texttt{el}}$, that is

[TABLE]

Second, the trace of the solution on $\Gamma\cup\Gamma_{N}$ is computed by solving a global problem given by the above introduced transmission conditions

[TABLE]

where the first condition is automatically fulfilled owing to the Dirichlet boundary condition $u_{e}{=}\hat{u}$ on $\partial\Omega_{e}\setminus\Gamma_{D}$ imposed in the local problem and to the uniqueness of the hybrid variable on each mesh edge (respectively, face in 3D).

Following the rationale in Section 2.1, the discrete functional spaces

[TABLE]

are introduced for the HDG approximation and $\mathcal{P}^{k}(\Omega_{e})$ and $\mathcal{P}^{k}(\Gamma_{i})$ here denote the spaces of polynomial functions of complete degree at most $k{\geq}1$ in $\Omega_{e}$ and on $\Gamma_{i}$, respectively.

Introduce the definition of the trace of the numerical normal flux

[TABLE]

where $\tau$ is a stabilization parameter critical to ensure stability and convergence of the HDG method, as extensively studied in a series of publications [38, 39, 40, 41, 42, 43] by Cockburn and coworkers.

The discrete HDG local problems are thus obtained integrating by parts the weak form of Equation (6) and exploiting the definition in Equation (8i). Following [72] and integrating the first equation by parts once and the second one twice, the resulting problem is: for $e{=}1,\ldots,\texttt{n}_{\texttt{el}}$, given $u_{D}$ on $\Gamma_{D}$ and $\hat{u}^{h}$ on $\Gamma\cup\Gamma_{N}$, find (u_{e}^{h},\text{\boldmathq\unboldmath}_{e}^{h})\in\mathcal{W}^{h}(\Omega_{e}){\times}\left[\mathcal{W}^{h}(\Omega_{e})\right]^{\texttt{n}_{\texttt{sd}}} such that

[TABLE]

for all (v,\text{\boldmathw\unboldmath})\in\mathcal{W}^{h}(\Omega_{e}){\times}\left[\mathcal{W}^{h}(\Omega_{e})\right]^{\texttt{n}_{\texttt{sd}}}.

Similarly, the discrete form of the HDG global problem (7) is: find \hat{u}^{h}\in\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{\mathcal{W}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(\Gamma\cup\Gamma_{N}) such that

[TABLE]

for all \hat{v}\in\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{\mathcal{W}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(\Gamma\cup\Gamma_{N}).

For the sake of readability and except in case of ambiguity, the superscript h associated with the numerical counterpart of the unknowns in the continuous spaces and the subscript e related to the HDG elemental approximations will be henceforth omitted.

2.3 Hybrid coupling based on Nitsche’s method

Consider a splitting of the domain $\Omega$ under analysis in two nonoverlapping subdomains ${\Omega}_{\texttt{CG}}$ and ${\Omega}_{\texttt{HDG}}$ such that $\Omega{=}{\Omega}_{\texttt{CG}}\cup{\Omega}_{\texttt{HDG}}$. The interface is defined as $\Gamma_{I}{:=}{\overline{\Omega}}_{\texttt{CG}}\cap{\overline{\Omega}}_{\texttt{HDG}}$. In the following, the subscripts ${}_{\texttt{CG}}$ and ${}_{\texttt{HDG}}$ are employed to identify the quantities associated with the CG and HDG discretizations, respectively.

The degrees of freedom of the coupled problem are displayed in Figure 1: for the CG subdomain, the unknown $u$ is depicted in red, whereas the degrees of freedom of HDG local and global problems are represented by blue circles and blue squares, respectively.

The coupled CG-HDG strategy is obtained starting from the discretizations described in Section 2.1 and 2.2, with two additional conditions enforcing continuity of the solution and of the normal flux across the interface $\Gamma_{I}$, namely

[TABLE]

where ${u}_{\texttt{CG}}$ is solution of the CG problem in Equation (4), ${\hat{u}}_{\texttt{HDG}}$ is obtained from the HDG global problem in Equation (8k), the numerical flux \text{\boldmathn\unboldmath}_{I}{\cdot}{\text{\boldmathq\unboldmath}}_{\texttt{HDG}} stems from the definition in Equation (8i) and \text{\boldmathn\unboldmath}_{I} is the outer normal to the domain ${\Omega}_{\texttt{HDG}}$. It is worth noticing that the condition imposing the equilibrium of the flux in Equation (8l) already accounts for the uniqueness of the unit normal vector on the interface $\Gamma_{I}$. Henceforth, and unless in case of ambiguity, on the interface $\Gamma_{I}$ the normal $n$ is assumed to be the outer direction to the subdomain under analysis, that is \text{\boldmathn\unboldmath}{=}\text{\boldmathn\unboldmath}_{I} for ${\Omega}_{\texttt{HDG}}$ and \text{\boldmathn\unboldmath}{=}{-}\text{\boldmathn\unboldmath}_{I} for ${\Omega}_{\texttt{CG}}$.

From a practical point of view, on the one hand, the HDG global unknown ${\hat{u}}_{\texttt{HDG}}$ is now defined also on the interface $\Gamma_{I}$ and the continuity of the solution is imposed as a Dirichlet boundary condition in the CG problem. On the other hand, the equilibrium of the flux is enforced via a Neumann boundary condition in the HDG global problem. Hence, the proposed technique does not affect neither the HDG local problem nor the core routines of the CG solver, leading to a minimally-intrusive coupling strategy.

Recall that ${\hat{u}}_{\texttt{HDG}}$ is solution of Equation (8k) on ${\Omega}_{\texttt{HDG}}$ and is thus known solely at the nodes of the HDG discretization. In this context, the Dirichlet condition ${u}_{\texttt{CG}}{=}{\hat{u}}_{\texttt{HDG}}$ is enforced in the CG problem via the well-known Nitsche’s method for the weak imposition of essential boundary conditions [11]. Moreover, this choice provides a flexible framework for the coupled discretization, allowing for nonconforming meshes at the interface $\Gamma_{I}$, as well as nonuniform polynomial approximations in ${\Omega}_{\texttt{CG}}$ and ${\Omega}_{\texttt{HDG}}$.

Following the rationale utilized for HDG in Equation (8i), the trace of the CG numerical normal flux on the interface $\Gamma_{I}$ is defined as

[TABLE]

where $h$ denotes the characteristic element size of the mesh discretization on $\Gamma_{I}$ and $\gamma$ is a sufficiently large positive parameter, commonly used to enforce coercivity of the discrete bilinear form in CG discretizations with Nitsche’s imposition of essential boundary conditions [2].

Exploiting the definition of the CG numerical normal flux in Equation (8m), the discrete form of the coupled CG-HDG solver given by Equations (4)-(8k)-(8l) is: find ({u}_{\texttt{CG}},{\hat{u}}_{\texttt{HDG}})\in\mathcal{V}^{h}({\Omega}_{\texttt{CG}}){\times}\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{\mathcal{W}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(\Gamma\cup\Gamma_{N}\cup\Gamma_{I}) such that

[TABLE]

for all (v,\hat{v})\in\mathcal{V}^{h}_{0}({\Omega}_{\texttt{CG}}){\times}\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{\mathcal{W}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(\Gamma\cup\Gamma_{N}\cup\Gamma_{I}). Note that the third and fourth terms on the left-hand side of Equation (8na) respectively enforce the symmetry and the coercivity of the resulting CG bilinear form, being $\gamma$ the above mentioned Nitsche’s parameter [2]. An alternative approach which does not require introducing the penalty term to guarantee the stability of the resulting numerical method is presented in [73] and relies on a nonsymmetric formulation of the discrete equation. It is worth noticing that in Equation (8na) $n$ refers to the outer normal to the domain ${\Omega}_{\texttt{CG}}$, whereas in Equation (8nb) $n$ is the outer normal to the element $\Omega_{e}\subset{\Omega}_{\texttt{HDG}}$.

The linear system arising from the discretization of the coupled problem in Equation (8n) is

[TABLE]

where the matrices ${\mathbf{K}}_{\texttt{CG}}$ and ${\mathbf{K}}_{\texttt{HDG}}$ are symmetric and feature the usual structure of the matrices of the CG and HDG global problem, respectively, whereas $\mathbf{K}_{\texttt{I}}$ is responsible for the coupling at the interface and stems from the last two terms on the right-hand side of Equations (8n).

It is worth noticing that despite the structure of the block matrix in Equation (8o), which is similar to the one presented in [50], its construction is extremely different. The strategies described in [50] couple CG and HDG discretizations at the elemental level. More precisely, they rely on the introduction of an appropriate projection operator to define the trace of the numerical flux and to impose the Dirichlet boundary condition on $\Gamma_{I}$ in the HDG local problems. This leads to the intrusive modification of either the block matrix of the elemental problem or the one of the global problem. On the contrary, the proposed coupling strategy solely relies on the hybrid variable ${\hat{u}}_{\texttt{HDG}}$ both to impose the Dirichlet boundary condition in the CG problem, see Equation (8na), and the Neumann one involving the numerical flux in the HDG global problem, see Equation (8nb). Moreover, the HDG local problems are not affected by the coupling and the implementation of the resulting strategy is thus minimally-intrusive with respect to existing CG and HDG solvers.

3 CG-HDG coupling for linear elastic problems

In this section, the coupling strategy is presented for the linear elasticity equation. The goal is to exploit in a minimally-intrusive way both the computational efficiency of CG approximations in presence of compressible materials and the accuracy and robustness of HDG discretizations in the nearly incompressible limit.

3.1 Governing equations

First, recall the governing equations of a linear elastic material

[TABLE]

where $u$ is the unknown displacement field and $\sigma$ is the Cauchy stress tensor. Equation (8p) enforces the equilibrium of linear and angular momentum of an elastic structure $\Omega$ subject to a volume force $f$, a tension $t$ on the boundary $\Gamma_{N}$ and an imposed displacement \text{\boldmathu\unboldmath}_{D} on $\Gamma_{D}$. For a linear elastic material, Hooke’s law describes the relationship between the Cauchy stress tensor $\sigma$ and the linearized strain rate tensor \text{\boldmath\varepsilon\unboldmath}(\text{\boldmathu\unboldmath}){:=}(\text{\boldmath\nabla\unboldmath}\text{\boldmathu\unboldmath}{+}\text{\boldmath\nabla\unboldmath}\text{\boldmathu\unboldmath}^{T})/2, namely

[TABLE]

where $\mathbf{I}_{\texttt{n}_{\texttt{sd}}}$ is the $\texttt{n}_{\texttt{sd}}{\times}\texttt{n}_{\texttt{sd}}$ identity matrix, $\operatorname{tr}(\cdot)$ is the trace operator and $(E,\nu)$ are the Young’s modulus and the Poisson’s ratio describing the mechanical properties of the material under analysis. Henceforth, the material is assumed to be homogeneous and isotropic inside the subdomains ${\Omega}_{\texttt{CG}}$ and ${\Omega}_{\texttt{HDG}}$. Thus, the above mentioned material coefficients depend neither on the spatial coordinate nor on the direction of the main strains.

Following [70], Equation (8p) is rewritten exploiting the Voigt notation for second-order symmetric tensors. The rationale of this approach is to enforce the symmetry of $\sigma$ and \text{\boldmath\varepsilon\unboldmath}(\text{\boldmathu\unboldmath}) pointwise by storing solely the $\texttt{m}_{\texttt{sd}}{:=}\texttt{n}_{\texttt{sd}}(\texttt{n}_{\texttt{sd}}{+}1)/2$ nonredundant components of the tensor, namely

[TABLE]

and

[TABLE]

The strain rate tensor is thus rewritten as \text{\boldmath\varepsilon\unboldmath}_{\texttt{V}}{=}\text{\boldmath\nabla\unboldmath}_{\!\texttt{S}}\text{\boldmathu\unboldmath}, where the symmetric gradient operator is expressed in matrix form by defining \text{\boldmath\nabla\unboldmath}_{\!\texttt{S}}\in\mathbb{R}^{\texttt{m}_{\texttt{sd}}{\times}\texttt{n}_{\texttt{sd}}} such that

[TABLE]

Moreover, the constitutive law is expressed in matrix form as \text{\boldmath\sigma\unboldmath}_{\texttt{V}}{=}\mathbf{D}\text{\boldmath\varepsilon\unboldmath}_{\texttt{V}}, where the fourth-order elasticity tensor linking $\sigma$ and \text{\boldmath\varepsilon\unboldmath}(\text{\boldmathu\unboldmath}), see Equation (8q), is rewritten by means of the $\texttt{m}_{\texttt{sd}}{\times}\texttt{m}_{\texttt{sd}}$ matrix

[TABLE]

where

[TABLE]

and the parameter $\vartheta$ denotes either a plane stress model ($\vartheta{=}1$) or a plane strain model ($\vartheta=2$) in 2D.

In a similar fashion, the Neumann boundary condition is formulated as \mathbf{N}^{T}\text{\boldmath\sigma\unboldmath}_{\texttt{V}}{=}\text{\boldmatht\unboldmath} by introducing the $\texttt{m}_{\texttt{sd}}{\times}\texttt{n}_{\texttt{sd}}$ matrix describing the outer unit normal vector to the boundary

[TABLE]

Hence, the linear elastic problem in Equation (8p)-(8q) is rewritten using Voigt notation as

[TABLE]

3.2 Hybrid coupling based on Nitsche’s method

Following the rationale discussed in Section 2, the strong form of the coupled CG-HDG problem is introduced. More precisely, in the CG subdomain the classical formulation discussed in [70] holds

[TABLE]

whereas in each element $\Omega_{e},\,e{=}1,\ldots,\texttt{n}_{\texttt{el}}$ of the HDG subdomain, the recently proposed HDG-Voigt formulation [68, 74] is considered

[TABLE]

The two methods are thus coupled by a set of conditions equivalent to the ones introduced in Equation (8l), namely

[TABLE]

where the first condition enforces the continuity of the displacement field and the second one the equilibrium of the normal traction at the interface, $\mathbf{N}_{I}$ being the outer normal to the domain ${\Omega}_{\texttt{HDG}}$. Henceforth, and unless in case of ambiguity, on the interface $\Gamma_{I}$ the normal $\mathbf{N}$ is assumed to be the outer direction to the subdomain under analysis, that is $\mathbf{N}{=}\mathbf{N}_{I}$ for ${\Omega}_{\texttt{HDG}}$ and $\mathbf{N}{=}{-}\mathbf{N}_{I}$ for ${\Omega}_{\texttt{CG}}$.

As previously discussed, the discrete form of the coupled problem is obtained by the weak form of the CG problem with Dirichlet boundary conditions imposed via Nitsche’s method on $\Gamma_{I}$ and by the HDG global problem, that is: find ({\text{\boldmathu\unboldmath}}_{\texttt{CG}},{\widehat{\text{\boldmathu\unboldmath}}}_{\texttt{HDG}})\in[\mathcal{V}^{h}({\Omega}_{\texttt{CG}})]^{\texttt{n}_{\texttt{sd}}}{\times} [\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{\mathcal{W}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(\Gamma\cup\Gamma_{N}\cup\Gamma_{I})]^{\texttt{n}_{\texttt{sd}}} such that

[TABLE]

for all (\text{\boldmathv\unboldmath},\widehat{\text{\boldmathv\unboldmath}})\in[\mathcal{V}^{h}_{0}({\Omega}_{\texttt{CG}})]^{\texttt{n}_{\texttt{sd}}}{\times}[\stackon[0pt]{\mathcal{W}}{\stretchto{\scalerel*[\widthof{\mathcal{W}}]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{0.5ex}_{\smash{\belowbaseline[-2pt]{\scriptstyle}}}}^{h}(\Gamma\cup\Gamma_{N}\cup\Gamma_{I})]^{\texttt{n}_{\texttt{sd}}}, where the definition of the trace of the numerical normal flux used for CG on the interface $\Gamma_{I}$ is

[TABLE]

and for HDG on all internal and boundary faces is

[TABLE]

As previously remarked, the proposed coupling only affects the HDG discretization in the global problem, whereas the local element-by-element problems present the usual structure of an HDG approximation. More precisely, the choice of the HDG formulation under analysis relies on its capability to achieve optimal convergence of the stress tensor and superconvergence of the displacement field using equal order polynomial approximation for all the variables. This is due to the definition of a pointwise symmetric mixed variable, namely the strain rate tensor, via Voigt notation [68, 74]. The resulting HDG local problems are: for $e{=}1,\dots,\texttt{n}_{\texttt{el}}$ find ({\text{\boldmathu\unboldmath}}_{\texttt{HDG}},{\text{\boldmathL\unboldmath}}_{\texttt{HDG}})\in[\mathcal{W}^{h}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}}{\times}[\mathcal{W}^{h}(\Omega_{e})]^{\texttt{m}_{\texttt{sd}}} such that

[TABLE]

for all (\text{\boldmathv\unboldmath},\text{\boldmathw\unboldmath})\in[\mathcal{W}^{h}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}}{\times}[\mathcal{W}^{h}(\Omega_{e})]^{\texttt{m}_{\texttt{sd}}}.

Moreover, given the discrete functional space

[TABLE]

where $\mathcal{P}^{k+1}(\Omega_{e})$ is the space of polynomial functions of complete degree at most $k{+}1$, a superconvergent approximation {\text{\boldmathu\unboldmath}}_{\texttt{HDG}}^{\star}\in[\mathcal{W}^{h}_{\star}(\Omega)]^{\texttt{n}_{\texttt{sd}}} of the displacement field is computed for each element $\Omega_{e},\,e{=}1,\dots,\texttt{n}_{\texttt{el}}$ by solving the postprocessed problem

[TABLE]

with the constraint

[TABLE]

to remove the underdetermination associated with rigid body translations and

[TABLE]

to treat rigid body rotations, where the curl operator \text{\boldmath\nabla\unboldmath}_{\!\texttt{W}}\in\mathbb{R}^{\texttt{n}_{\texttt{rr}}{\times}\texttt{n}_{\texttt{sd}}} and the tangential direction to the boundary $\mathbf{T}\in\mathbb{R}^{\texttt{n}_{\texttt{rr}}{\times}\texttt{n}_{\texttt{sd}}}$, with $\texttt{n}_{\texttt{rr}}{:=}\texttt{n}_{\texttt{sd}}(\texttt{n}_{\texttt{sd}}{-}1)/2$ being the number of rigid rotational body modes, are written in matrix form as

[TABLE]

and

[TABLE]

4 Numerical studies

In this section, several numerical examples are presented to test the optimal convergence properties of the proposed hybrid CG-HDG coupling. Special emphasis is devoted to highlight the advantages of this methodology in terms of accuracy, robustness and minimal intrusiveness of its implementation in existing CG and HDG libraries.

4.1 Two-dimensional thermal problem

The first example considers the problem in Equation (1), in two dimensions, to assess optimal convergence and robustness of the coupling approach to the involved parameters. The computational domain $\Omega{=}[-1,1]{\times}[-1,1]$ is decomposed in two nonoverlapping subdomains, namely ${\Omega}_{\texttt{CG}}{=}[0,1]{\times}[-1,1]$ and ${\Omega}_{\texttt{HDG}}{=}\Omega\setminus{\Omega}_{\texttt{CG}}$. The interface is thus identified by $\Gamma_{I}{=}\left\{(x,y)\in\mathbb{R}^{2}:\right.$ $\left.x=0\right\}$. The domain $\Omega$ is discretized using uniform meshes of triangular elements constructed by means of the mesh generator EZ4U [75, 76].The first three levels of mesh refinement are presented in Figure 2. Elements in red (resp., blue) belong to the CG (resp., HDG) subdomain, whereas the interface $\Gamma_{I}$ is drawn in black.

The source term is selected so that the analytical solution is

[TABLE]

and Dirichlet boundary conditions, corresponding to the analytical solution, are imposed on $\Gamma_{D}{=}\partial\Omega$.

The solution of the thermal problem computed on the fifth level of refinement of the mesh with polynomial approximation of degree $k{=}5$ in both the CG and HDG subdomains is displayed in Figure 3.

4.1.1 Sensitivity to Nitsche’s parameter

The influence of Nitsche’s parameter on the accuracy of the proposed hybrid CG-HDG coupling is investigated. The HDG stabilization parameter $\tau$ is considered constant in the subdomain ${\Omega}_{\texttt{HDG}}$ and equal to $10$. It is worth recalling that a stabilization parameter $\tau{=}C\kappa/\ell$, where $\kappa$ is the thermal conductivity, equal to $1$ for the problem under analysis, $\ell$ the characteristic length of the problem and $C$ a positive constant scaling factor, guarantees stability and optimal convergence of the HDG discretization for the Poisson equation [39].

Figure 4 displays the evolution of the error of the primal variable $u$ on the whole domain $\Omega$ as a function of Nitsche’s parameter $\gamma$ using the third level of mesh refinement and for different degrees of the polynomial approximation. For low values of Nitsche’s parameter, oscillations appear in the error, whereas stability is achieved choosing a sufficiently large $\gamma$. It is worth noticing that the minimum value of $\gamma$ guaranteeing stability of the numerical method varies with the degree of the polynomial approximation. Henceforth, $\gamma{=}10^{2}$ is considered for the following numerical experiments.

4.1.2 Optimal convergence of the CG-HDG coupling

The convergence of the error of the temperature measured in the $\mathcal{L}_{2}$ norm as a function of the characteristic element size $h$ is presented in Figure 5 for polynomial degree of approximation $k{=}1,\ldots,5$. Optimal convergence of order $k{+}1$ is achieved both in the CG subdomain (Fig. 5, left) and in the HDG one (Fig. 5, right) by the proposed hybrid CG-HDG coupling.

Comparable results in terms of accuracy are obtained using the coupling strategy discussed in [50] and the convergence studies are omitted for the sake of brevity. It is worth recalling that the main advantage of the proposed hybrid coupling is represented by its minimally-intrusive nature which makes this approach extremely easy to implement in existing CG and HDG libraries.

4.2 Two-dimensional elastic problem

The second example considers the linear elastic problem in Equation (8v) with compressible and nearly incompressible materials. The domain $\Omega{=}[-1,1]{\times}[-1,1]$ is decomposed in two nonoverlapping subdomains, namely ${\Omega}_{\texttt{CG}}{=}[-1,0]{\times}[-1,0]\cup[0,1]{\times}[0,1]$ and ${\Omega}_{\texttt{HDG}}{=}\Omega\setminus{\Omega}_{\texttt{CG}}$. The interface is thus identified by $\Gamma_{I}{=}\left\{(x,y)\in\mathbb{R}^{2}:\right.$ $\left.x=0\text{ or }y=0\right\}$. The computational meshes are constructed as for the above described thermal problem. The first three levels of mesh refinement are presented in Figure 6: the subdomain ${\Omega}_{\texttt{CG}}$ in displayed in red, ${\Omega}_{\texttt{HDG}}$ in blue, whereas the interface $\Gamma_{I}$ is drawn in black. It is worth noticing that for the problem under analysis both the CG and the HDG subdomains are not simply connected

The source term is selected so that the analytical solution of the problem is \text{\boldmathu\unboldmath}{=}(u_{x},u_{y}) such that

[TABLE]

and Dirichlet boundary conditions, corresponding to the analytical solution, are imposed on $\Gamma_{D}{=}\partial\Omega$. This numerical experiment is inspired by the work in [77] and considers inhomogeneous mechanical properties in the two subdomains, namely a compressible and stiff material with Young’s modulus $E{=}250$ and Poisson’s ratio $\nu{=}0.3$ in the red domain ${\Omega}_{\texttt{CG}}$ and a nearly incompressible and soft one such that $E{=}25$ and $\nu{=}0.49999$ in the blue domain ${\Omega}_{\texttt{HDG}}$.

4.2.1 Locking-free approximation of the CG-HDG coupling

Four strategies are considered for the discretization of the previously introduced elastic problem:

CG discretization in the whole domain $\Omega$; 2. 2.

HDG discretization in the whole domain $\Omega$; 3. 3.

hybrid CG-HDG discretization with polynomial approximation of degree $k$ both in ${\Omega}_{\texttt{CG}}$ and ${\Omega}_{\texttt{HDG}}$; 4. 4.

hybrid CG-HDG discretization with polynomial approximation of degree $k{+}1$ in ${\Omega}_{\texttt{CG}}$ and $k$ with local HDG postprocess in ${\Omega}_{\texttt{HDG}}$.

It is worth noticing that the presence of nearly incompressible regions makes a CG discretization unfeasible in the whole domain, whereas employing HDG everywhere leads to a higher computational cost than CG in the subdomain ${\Omega}_{\texttt{CG}}$. Coupling the two approaches allows to devise a robust and flexible numerical scheme in $\Omega$, especially benefitting from nonuniform polynomial degrees of approximation.

For the following simulations, the HDG stabilization parameter $\tau$ is considered constant and equal to $2.5{\times}10^{2}$. It is worth recalling that a stabilization parameter $\tau{=}CE/\ell$, where $E$ is the Young’s modulus, spanning from $25$ to $250$ for the problem under analysis, $\ell$ the characteristic length of the problem and $C$ a positive constant scaling factor guarantees stability and optimal convergence of the HDG discretization for the linear elasticity equation [69]. The Nitsche’s parameter is set to $\gamma{=}2.5{\times}10^{3}$ for the coupling based on uniform degree of approximation, whereas $\gamma{=}2.5{\times}10^{4}$ is considered for the case of nonuniform approximations.

The convergence of the error of the displacement field measured in the $\mathcal{L}_{2}$ norm as a function of the characteristic element size $h$ is presented in Figure 7, for polynomial degree of approximation $k{=}1,\ldots,3$.

Due to the presence of a nearly incompressible material, the CG approximation suffers from classical locking phenomena [52] which prevent convergence for $k{=}1$, whereas suboptimal rates of order $k$ are obtained for higher-degree of polynomial approximation (Fig. 7a). On the contrary, both HDG (Fig. 7b) and the hybrid ${\text{CG}}_{k}$-${\text{HDG}}_{k}$ coupling (Fig. 7c) present optimal convergence of order $k{+}1$ for the displacement field without locking effects [62] when polynomial of degree $k$ are utilized in both ${\Omega}_{\texttt{CG}}$ and ${\Omega}_{\texttt{HDG}}$. Moreover, considering a nonuniform degree of approximation with polynomial functions of degree $k{+}1$ in ${\Omega}_{\texttt{CG}}$ and $k$ in ${\Omega}_{\texttt{HDG}}$ and exploiting the HDG postprocess strategy in Equation (8ad), a displacement field superconverging with order $k{+}2$ in the whole domain $\Omega$ is obtained (Fig. 7d).

Figure 8 displays the convergence of the error of the stress field measured in the $\mathcal{L}_{2}$ norm as a function of the characteristic element size $h$ for polynomial degree of approximation $k{=}1,\ldots,3$.

The CG approximation of the stress field being computed as a postprocess of the displacement field leads to unreliable results due to the locking effects (Fig. 8a). On the contrary, mixed formulations provide a direct approximation of the stress tensor [67]. More precisely, HDG exploits the definition of the pointwise symmetric mixed variable {\text{\boldmathL\unboldmath}}_{\texttt{HDG}} in [68, 74] to obtain optimal convergence of the stress with order $k{+}1$ (Fig. 8b).

To reduce the computational cost of the approximation of the linear elastic problem in the whole domain $\Omega$, a CG discretization is considered in the compressible subdomain ${\Omega}_{\texttt{CG}}$ and an HDG one in the nearly incompressible region ${\Omega}_{\texttt{HDG}}$. The approximation of the stress computed using the ${\text{CG}}_{k}$-${\text{HDG}}_{k}$ coupling in Figure 8c presents suboptimal convergence of order $k$ in $\Omega$ since $\sigma$ is computed as a postprocess of the primal solution {\text{\boldmathu\unboldmath}}_{\texttt{CG}} in the subdomain ${\Omega}_{\texttt{CG}}$. Optimal convergence of the stress field in the whole domain $\Omega$ is recovered by constructing the hybrid ${\text{CG}}_{k{+}1}$-${\text{HDG}}_{k}$ coupling based on a nonuniform degree of approximation with polynomial functions of degree $k{+}1$ in ${\Omega}_{\texttt{CG}}$ and $k$ in ${\Omega}_{\texttt{HDG}}$ (Fig. 8d).

5 Application to multimaterial elastostatics problems

In this section, the proposed hybrid CG-HDG coupling is applied to linear elastostastics problems of engineering interest. More precisely, problems with domains featuring multiple materials and inhomogeneous mechanical properties are investigated.

5.1 Bimaterial Cook’s membrane problem

The classical bending-dominated problem of Cook’s membrane [78] is revisited to take into account two materials with different mechanical properties as proposed in [77]. The tapered plate $\Omega$ defined by the points $(0,0)$, $(48,44)$, $(48,60)$, $(0,44)$ is decomposed in two nonoverlapping subdomains such that ${\Omega}_{\texttt{CG}}$ is the convex region identified by the points $(12,20.25)$, $(36,38.75)$, $(36,50.25)$, $(12,38.75)$ and ${\Omega}_{\texttt{HDG}}{=}\Omega\setminus{\Omega}_{\texttt{CG}}$. The first three levels of mesh refinement are displayed in Figure 9, where the subdomain ${\Omega}_{\texttt{CG}}$ is represented in red, ${\Omega}_{\texttt{HDG}}$ in blue and the interface $\Gamma_{I}$ in black.

The subdomain ${\Omega}_{\texttt{CG}}$ features a compressible and stiff material with Young’s mudulus $E{=}250$ and Poisson’s ratio $\nu{=}0.35$, whereas the properties of the nearly incompressible region ${\Omega}_{\texttt{HDG}}$ are $E{=}80$ and $\nu$ spanning among the values $\left\{0.49,0.4999,0.499999\right\}$. The plate is clamped on the left and is subject to a shear load uniformly distributed along the positive $y$-direction with total load $F{=}100$ on the right. The remaining boundaries are free surfaces on which a homogenous Neumann condition is imposed.

The bimaterial Cook’s membrane problem is solved using the four strategies discussed in Section 4.2.1. The HDG stabilization parameter $\tau$ is considered constant an equal to $10$, whereas the Nitsche’s parameter $\gamma$ is set to $10^{4}$. The HDG approximation being locking-free [62] is considered as reference solution for the problem under analysis.

Figure 10 compares the approximation of the displacement and the stress fields in the nearly incompressible case of $\nu{=}0.4999$, for the different discretization strategies described above, using polynomial of degree $k{=}1$ on the fifth level of mesh refinement. It is straightforward to observe that CG suffers from locking phenomena providing an unreliable approximation of both displacement and stress fields. The two coupling strategies provide comparable results for the displacement field, whereas a slightly more accurate approximation of the stress field is achieved by the ${\text{CG}}_{k{+}1}$-${\text{HDG}}_{k}$ coupling in the compressible region. This is due to the increased accuracy of the CG approximation in the subdomain ${\Omega}_{\texttt{CG}}$, already observed for this approach in Section 4.2.1.

The evolution of the vertical displacement of the top corner of the right end of the plate is displayed in Figure 11 as a function of the number of mesh elements in the domain $\Omega$, for different values of the Poisson’s ratio.

The CG approximation with polynomial functions of degree $k{=}1$ leads to a significant underestimation of the displacement of the tip and this error increases as the Poisson’s ratio tends to $0.5$ due to locking effects (Fig. 11c). A good approximation of the displacement of the tip is recovered using both coupling strategies. As highlighted above, the ${\text{CG}}_{k{+}1}$-${\text{HDG}}_{k}$ approach exploits the additional accuracy of the CG discretization in ${\Omega}_{\texttt{CG}}$ to construct a more reliable approximation, even when coarser meshes are utilized. Moreover, it is worth noticing that the CG-HDG couplings are locking-free and no loss of accuracy is experienced in the incompressible limit (Fig. 11c).

5.2 Three-dimensional laminated composite beam

The last example considers a three-dimensional laminated composite beam to show the capability of the proposed coupling to treat problems of interest for engineering applications. The beam $\Omega{=}[-1,1]{\times}[-1,1]{\times}[0,10]$ consists of four layers with alternating mechanical properties, namely a compressible isotropic material with Young’s modulus $E=10$ and Poisson’s ratio $\nu=0.3$ in ${\Omega}_{\texttt{CG}}{=}[-1,1]{\times}[-0.5,0]{\times}[0,10]\cup[-1,1]{\times}[0.5,1]{\times}[0,10]$ and a nearly incompressible one with $E=1$ and $\nu=0.49999$ in ${\Omega}_{\texttt{HDG}}{=}\Omega\setminus{\Omega}_{\texttt{CG}}$. The interfaces between the two materials are identified by the surfaces $\Gamma_{I}{=}\left\{(x,y,z)\in\mathbb{R}^{3}:\right.$ $\left.y{=}{-}0.5,y{=}0\text{ or }y{=}0.5\right\}$.

The cantilever beam is clamped on the left surface and a load $F{=}1.3{\times}10^{-3}$ is uniformly distributed on the top surface of the beam, along the negative $y$-direction. The remaining surfaces are free boundaries on which a homogenous Neumann condition is imposed.

A reference solution is computed using HDG on a mesh featuring $414,720$ tetrahedral elements with a polynomial approximation of degree $k{=}1$, for a total of $7,568,640$ global unknowns representing the displacement on the faces. The HDG stabilization parameter $\tau$ is considered constant an equal to 10, whereas the Nitsche’s parameter $\gamma$ is set to $10^{2}$.

Contrary to the previous strategy based on nonuniform polynomial approximations, in this example the degree is set to $k{=}1$ in all $\Omega$ and the mesh in ${\Omega}_{\texttt{CG}}$ is refined. The goal is to show that an improved accuracy of the discrete solution is achieved by increasing the number of degrees of freedom in the subdomain ${\Omega}_{\texttt{CG}}$, whereas accurate results are obtained by the corresponsing HDG discretization on the coarse grid. More precisely, two mesh configurations are considered. The coarse mesh counts $7,680$ tetrahedra both in ${\Omega}_{\texttt{CG}}$ and ${\Omega}_{\texttt{HDG}}$, whereas the fine discretization respectively features $30,720$ and $7,680$ elements in the two subdomains, as displayed in Figure 12. The subdomain ${\Omega}_{\texttt{CG}}$ is represented in red, ${\Omega}_{\texttt{HDG}}$ in blue and the interfaces $\Gamma_{I}$ in black.

Figure 13 displays the comparison of the displacement field computed using the proposed ${\text{CG}}_{k}$-${\text{HDG}}_{k}$ coupling strategy with a polynomial approximation of degree $k{=}1$ on the meshes presented in Figure 12. The results are qualitatively comparable for both approaches. Nonetheless, a small discrepancy is observed in the quantitative evaluation of the vertical displacement of the axis at $x{=}0$ and $y{=}0$ along the beam length (Fig. 14). More precisely, the ${\text{CG}}_{k}$-${\text{HDG}}_{k}$ coupling strategy applied to the configuration with a coarse mesh in the subdomain ${\Omega}_{\texttt{CG}}$ introduces a relative error for the vertical displacement and the bending moment at midspan of $8\%$, whereas using the fine mesh in ${\Omega}_{\texttt{CG}}$ they drop to $4\%$ and $2\%$, respectively. For both quantities of interest, the HDG approximation using the coarse mesh of ${\Omega}_{\texttt{CG}}$ provides a relative error smaller than $1\%$. On the contrary, the CG approximation severely underestimates the vertical displacement of the axis.

6 Concluding remarks

A novel strategy to couple CG and HDG discretizations based on Nitsche’s method has been presented in the context of linear thermal and elastic problems. Special emphasis has been devoted to the case of linear elasticity in which both the computational efficiency of CG in the simulation of compressible materials and the robustness of HDG in the incompressible limit are exploited.

The proposed coupling imposes continuity of the solution and of the trace of the normal numerical flux across the interface between the CG and HDG subdomains. The coupled degrees of freedom are only the ones associated with the HDG hybrid variable on the interface. The remaining terms in the formulation, namely the CG discrete matrix and the HDG local and global ones are not modified, leading to a minimally-intrusive implementation in existing CG and HDG libraries.

Numerical examples have been used to verify the optimal convergence of the proposed methodology in the global domain, as well as in the CG and HDG subdomains, separately. Moreover, a strategy based on nonuniform polynomial degree approximation with a CG discretization of degree $k{+}1$ and an HDG one of degree $k$ has been proved to lead to a global convergence of the flux/stress of order $k{+}1$ and a superconvergence of the solution with order $k{+}2$ by exploiting the inexpensive element-by-element HDG postprocess strategy proposed in [68, 74].

Acknowledgements

This work is partially supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie actions (Grant No. 675919), the Spanish Ministry of Economy and Competitiveness (Grant No. DPI2017-85139-C2-2-R) and the Generalitat de Catalunya (Grant No. 2017-SGR-1278). Andrea La Spina is supported by the European Education, Audiovisual and Culture Executive Agency (EACEA) under the Erasmus Mundus Joint Doctorate Simulation in Engineering and Entrepreneurship Development (SEED), FPA 2013-0043.