Enriched virtual elements for plane elasticity with corner singularities

Abstract

We construct a nonconforming virtual element method for the approximation of singular solutions to isotropic linear elasticity problems on polygonal domains. Standard nonconforming virtual element spaces are enriched with suitable singular functions. The enrichment is based on the nonconforming structure of the discrete spaces and not on partition of unity techniques. We prove optimal convergence and assess numerically the theoretical results of the method. The proposed scheme naturally paves the way for an efficient linear elastic fracture solver.

On some corner singularities for the isotropic linear elasticity problem on polygons

Recall splitting (4) and denote the measure of the angle at $\textbf{A}$ by $\omega$. In the following, we review the behaviour of the singular components ${\mathscr {S}}_{\textbf{A}}$ in the neighbourhood of vertex $\textbf{A}$ of $\Omega$. Such singular functions are typically given by the sum of several singular terms. In Sect. A.1, we focus on each single term of such expansion on a couple of benchmark domains; notably, we exhibit possible choices of such terms. Instead, in Sect. A.2, we review the full singular expansion.

1.1 Singular functions at corners

We compute explicitly singular solutions at corners in the kernel of the operator $\textbf{M}$ defined in (1) of the form (5). These solutions represents singular displacement fields, corresponding to self-equilibrated stress states for the 2D material body exhibiting a re-entrant corner, which tend to infinity at the tip of the corner. Notably, the parameter $\alpha$, and the functions $\mathscr {S}_{1}$ and $\mathscr {S}_{2}$ in (5) are computed in a two step procedure:

• impose the constraint

  $$\begin{aligned} \textbf{M}(r^{\alpha } (\mathscr {S}_{1}(\theta ), \mathscr {S}_{2}(\theta ))^T) = \textbf{0}\quad \text {in } \Omega ; \end{aligned}$$

  (32)

• impose the boundary conditions for $\theta =0$ and $\theta =\omega$.

Introduce

$$\begin{aligned} F(\alpha ):= \lambda + 3\mu - \alpha (\lambda + \mu ), \quad G(\alpha ):= \lambda + 3\mu + \alpha (\lambda + \mu ). \end{aligned}$$

As in [30, equation (4.7)], imposing (32) for functions in the form (5) yields, for $\alpha \in \mathbb R^+$ to be fixed later,

$$\begin{aligned} \begin{aligned} {\mathscr {S}}_{\textbf{A}}(\theta )&= c_1(\alpha ) \begin{pmatrix} \cos [(1+\alpha ) \theta ] \\ -\sin [(1+\alpha ) \theta ] \end{pmatrix} \\&\quad + c_2(\alpha ) \begin{pmatrix} \sin [(1+\alpha ) \theta ] \\ \cos [(1+\alpha ) \theta ] \end{pmatrix} \\&\quad + c_3(\alpha ) \begin{pmatrix} F(\alpha ) \cos [(1 - \alpha ) \theta ] \\ -G(\alpha ) \sin [(1 - \alpha ) \theta ] \end{pmatrix} \\&\quad + c_4(\alpha ) \begin{pmatrix} F(\alpha ) \sin [(1 - \alpha ) \theta ] \\ G(\alpha ) \cos [(1 - \alpha ) \theta ] \end{pmatrix}. \end{aligned} \end{aligned}$$

(33)

So far, $\alpha$ denotes an arbitrary singular exponents at corner $\textbf{A}$, which has to be fixed so that the function in (33) fulfils suitable homogeneous boundary conditions. Explicit choices of $\alpha$ are determined solving a $4\times 4$ system with unknowns given by $c_j(\alpha )$, $j=1,\dots , 4$, which is obtained imposing the homogeneous boundary conditions.

Recall that $\omega$ denotes the measure of $\textbf{A}$ and that the angular part of the polar coordinates takes value in $(0,\omega )$. For the sake of exposition, we focus on Dirichlet boundary conditions only, which correspond to clamped edges of the re-entrant corner.

Imposing homogeneous Dirichlet boundary conditions in (33), the resulting $4\times 4$ system reads

$$\begin{aligned} \begin{pmatrix} \cos [(1+\alpha )0] &{} \sin [(1+\alpha )0] &{} F(\alpha ) \cos [(1-\alpha )0] &{} F(\alpha ) \sin [(1-\alpha )0]\\ \cos [(1+\alpha )\omega ] &{} \sin [(1+\alpha )\omega ] &{} F(\alpha ) \cos [(1-\alpha )\omega ] &{} F(\alpha ) \sin [(1-\alpha )\omega ]\\ -\sin [(1+\alpha )0] &{} \cos [(1+\alpha )0] &{} -G(\alpha ) \sin [(1-\alpha )0] &{} G(\alpha ) \cos [(1-\alpha )0]\\ -\sin [(1+\alpha )\omega ] &{} \cos [(1+\alpha )\omega ] &{} -G(\alpha ) \sin [(1-\alpha ) \omega ] &{} G(\alpha ) \cos [(1-\alpha ) \omega ]\\ \end{pmatrix} \begin{pmatrix} c_1(\alpha )\\ c_2(\alpha )\\ c_3(\alpha )\\ c_4(\alpha ) \end{pmatrix}= \begin{pmatrix} 0\\ 0\\ 0\\ 0 \end{pmatrix}, \end{aligned}$$

i.e.,

$$\begin{aligned} \begin{pmatrix} 1 &{} 0 &{} F(\alpha ) &{} 0 \\ \cos [(1+\alpha )\omega ] &{} \sin [(1+\alpha )\omega ] &{} F(\alpha ) \cos [(1-\alpha )\omega ] &{} F(\alpha ) \sin [(1-\alpha )\omega ]\\ 0 &{} 1 &{} 0 &{} G(\alpha ) \\ -\sin [(1+\alpha )\omega ] &{} \cos [(1+\alpha )\omega ] &{} -G(\alpha ) \sin [(1-\alpha ) \omega ] &{} G(\alpha ) \cos [(1-\alpha ) \omega ]\\ \end{pmatrix} \begin{pmatrix} c_1(\alpha )\\ c_2(\alpha )\\ c_3(\alpha )\\ c_4(\alpha ) \end{pmatrix}= \begin{pmatrix} 0\\ 0\\ 0\\ 0 \end{pmatrix}. \end{aligned}$$

(34)

System (34) admits a nontrivial solution if the determinant of the system matrix is zero. Imposing such a determinant to be zero is equivalent to solving the so-called characteristic equation.

In the case of homogeneous Dirichlet, i.e., clamped boundary conditions, the characteristic equation reads as follows, see, e.g., [30, Sect. 4],¹ find $\alpha \in \mathbb {R}_+$ such that

$$\begin{aligned} \sin ^2(\alpha \omega ) - \alpha ^2 \lambda \frac{\lambda +2\mu }{(\lambda +3\mu )^2} \sin ^2 \omega = 0. \end{aligned}$$

(35)

Next, we look for explicit $\alpha$ solutions to (35) on two benchmark domains.

Example A.1

We focus on the re-entrant corner (0, 0) of the L-shaped domain

$$\begin{aligned} \Omega _1:= (-1,1)^2 \setminus [0,1) \times (-1,0], \end{aligned}$$

and fix the following Lamé parameters:

$$\begin{aligned} \lambda = -\frac{7}{2} + 9 \frac{\sqrt{2}}{2}, \quad \mu = \frac{7}{2} - 3 \frac{\sqrt{2}}{2}. \end{aligned}$$

We have $\omega = 3\pi /2$; the first three real smallest solutions to (35) read

$$\begin{aligned} \begin{aligned} \alpha _1&= 0.58917798046021; \quad \alpha _2 = 0.769651906609214; \\ \alpha _3&= 7/6. \end{aligned} \end{aligned}$$

The corresponding characteristic equation has a finite number of real solutions.

Next, we show how to recover the singular function in (5). To the aim, we only consider the case $\alpha =\alpha _3$. The coefficients $c_j(7/6)$, $j=1,\dots ,4$, in (33) are found by solving the system (34) for the given values of $\alpha$, $\lambda$, and $\mu$. In particular, the rank-3 system is

$$\begin{aligned}{} & {} \begin{pmatrix} 1 &{} 0 &{} 7\left( 1-\frac{\sqrt{2}}{2} \right) &{} 0 \\ -\frac{\sqrt{2}}{2} &{} -\frac{\sqrt{2}}{2} &{} 7\left( 1-\frac{\sqrt{2}}{2} \right) \frac{\sqrt{2}}{2} &{} -7 \left( 1-\frac{\sqrt{2}}{2} \right) \frac{\sqrt{2}}{2}\\ 0 &{} 1 &{} 0 &{} 7 \left( 1 + \frac{\sqrt{2}}{2} \right) \\ \frac{\sqrt{2}}{2} &{} -\frac{\sqrt{2}}{2} &{} 7 \left( 1 + \frac{\sqrt{2}}{2} \right) \frac{\sqrt{2}}{2} &{} 7 \left( 1 + \frac{\sqrt{2}}{2} \right) \frac{\sqrt{2}}{2}\\ \end{pmatrix}\\{} & {} \quad \times \begin{pmatrix} c_1(7/6)\\ c_2(7/6)\\ c_3(7/6)\\ c_4(7/6) \end{pmatrix}= \begin{pmatrix} 0\\ 0\\ 0\\ 0 \end{pmatrix}. \end{aligned}$$

There are infinite solutions to this problem: for any $c_4(7/6)$,

$$\begin{aligned} \begin{aligned}&c_1(7/6) = 7 \left( 1-\frac{\sqrt{2}}{2} \right) (\sqrt{2} + 1) c_4(7/6), \\&c_2(7/6) = -7\left( 1+\frac{\sqrt{2}}{2} \right) c_4(7/6), \\&c_3(7/6) = -(\sqrt{2} + 1) c_4(7/6). \end{aligned} \end{aligned}$$

In particular, the explicit singular solution belongs to$H^{\frac{13}{6} - \varepsilon }(\Omega )$, for all arbitrarily small, positive $\varepsilon$, and is a multiple of

$$\begin{aligned} r^{\frac{7}{6}} {\mathscr {S}}_{\textbf{A}}(\theta )= & {} r^{\frac{7}{6}} \left[ 7 (1 - \sqrt{2} / 2 ) (\sqrt{2} + 1)\begin{pmatrix} \cos ( 13 \theta /6 ) \nonumber \\ -\sin ( 13 \theta /6) \end{pmatrix} \right. \\{} & {} \left. -\,7 (1 + \sqrt{2} / 2 ) \begin{pmatrix} \sin ( 13 \theta /6) \\ \cos ( 13 \theta /6) \end{pmatrix} \right. \nonumber \\{} & {} \left. -\,(\sqrt{2} + 1) \begin{pmatrix} 7\left( 1 - \frac{\sqrt{2}}{2}\right) \cos ( \theta /6) \\ + 7\left( 1 + \frac{\sqrt{2}}{2}\right) \sin ( \theta /6) \end{pmatrix} \right. \nonumber \\{} & {} \left. + \,\begin{pmatrix} - 7\left( 1 - \frac{\sqrt{2}}{2}\right) \sin (\theta /6) \\ 7\left( 1 + \frac{\sqrt{2}}{2}\right) \cos (\theta /6) \end{pmatrix} \right] . \end{aligned}$$

(36)

Example A.2

We focus on the tip of the the slit domain

$$\begin{aligned} \Omega _2:= (-1,-1)^2 \setminus \{ \{ 0 \} \times [0,1) \}. \end{aligned}$$

In this case,  $\omega = 2 \pi$ and (35) simply reads

$$\begin{aligned} \sin ^2(2\pi \alpha ) = 0. \end{aligned}$$

This equation has positive solutions $\alpha$ given by

$$\begin{aligned} 2\pi \alpha = k \pi \quad \forall k \in \mathbb N \quad \quad \longrightarrow \alpha = \frac{k}{2} \quad \forall k \in \mathbb N. \end{aligned}$$

Interestingly, all the singular values $\alpha$ do not depend on the Lamé coefficients.

For the sake of exposition, we focus on the strongest singularity $\alpha = 1/2$. Compute

$$\begin{aligned} F(1/2) = \frac{1}{2} \lambda + \frac{5}{2} \mu , \quad \quad G(1/2) = \frac{3}{2} \lambda + \frac{7}{2} \mu . \end{aligned}$$

The coefficients $c_j(\alpha )$, $j=1,\dots ,4$, in (33) are found by solving the system (34) for the given $\alpha$, and for all $\lambda$ and $\mu$. In particular, the rank-2 system is

$$\begin{aligned} \begin{pmatrix} 1 &{} 0 &{} F(1/2) &{} 0 \\ -1 &{} 0 &{} -F(1/2) &{} 0 \\ 0 &{} 1 &{} 0 &{} G(1/2) \\ 0 &{} -1 &{} 0 &{} -G(1/2) \\ \end{pmatrix} \begin{pmatrix} c_1(1/2)\\ c_2(1/2)\\ c_3(1/2)\\ c_4(1/2) \end{pmatrix}= \begin{pmatrix} 0\\ 0\\ 0\\ 0 \end{pmatrix}. \end{aligned}$$

(37)

There are two families of solutions to the linear system (34). They correspond to the two following choices of the coefficients in (33):

$$\begin{aligned} c_1(1/2)&= - F(1/2) c_3(1/2) \quad \quad \forall c_3 (1/2) \in \mathbb R, \\ c_2(1/2)&= - G(1/2) c_3(1/2) \quad \quad \forall c_4 (1/2) \in \mathbb R. \end{aligned}$$

Thus, the singular solutions are

$$\begin{aligned}&r^{\frac{1}{2}} {\mathscr {S}}_{\textbf{A}}^1(\theta ) \nonumber \\&\quad = r^{\frac{1}{2}} \left[ -F(1/2) \begin{pmatrix} \cos ( 3/2 \theta ) \\ -\sin ( 3/2 \theta ) \end{pmatrix} + \begin{pmatrix} F(1/2) \sin ( 1/2 \theta ) \\ -G(1/2) \cos ( 1/2 \theta ) \end{pmatrix} \right] \end{aligned}$$

(38)

and

$$\begin{aligned}&r^{\frac{1}{2}} {\mathscr {S}}_{\textbf{A}}^2(\theta ) \nonumber \\&\quad = r^{\frac{1}{2}} \left[ -G(1/2) \begin{pmatrix} \sin ( 3/2 \theta ) \\ \cos ( 3/2 \theta ) \end{pmatrix} + \begin{pmatrix} F(1/2) \sin ( 1/2 \theta ) \\ G(1/2) \cos ( 1/2 \theta ) \end{pmatrix} \right] . \end{aligned}$$

(39)

The functions in (38) and (39) belong to $H^{\frac{3}{2} - \varepsilon }(\Omega )$, for all arbitrarily small, positive $\varepsilon$. They are found in linear elastic fracture mechanics as reported, e.g., in the classical paper [28].

The Lamé parameters do not come into play in the singular exponents $\alpha$, but only in the linear combination appearing in the above singular functions. This is not the case of Example A.1, where also the singular exponent $\alpha$ used to depend on the Lamé parameters.

Besides, there are infinitely many singular exponents $\alpha$. Instead, in Example A.1, there is only a finite number of real exponents $\alpha$.

Remark 1

So far, we mainly focused on the case of Dirichlet and homogeneous Neumann boundary conditions on the two edges abutting vertex $\textbf{A}$. The above way of reasoning extends to all possible combinations of boundary conditions. We refer to [30] for the counterpart of the characteristic equation to solve. The relevant point is that, in the analysis of the method, the actual representation of the singular function is not relevant as long as it satisfies (6).

1.2 General singular expansions at vertices

We review the full singular expansions at a vertex as in [30]; see also Dauge [16], Grisvard [20], Maz’ya and Plamenevskii [27].

Let $\alpha _0$ be such that (32) is satisfied and let $\textbf{e}_0 = \textbf{e}_0 (\alpha ,\varphi )$ be an associated eigenfunction. The set of fields $\{ \textbf{e}_{0,0},\, \textbf{e}_{0,1},\, \dots , \textbf{e}_{0,k} \}$ with $\textbf{e}_{0,0}=\textbf{e}_0$ is called a Jordan chain to $\alpha$ if

$$\begin{aligned} \sum _{q=0}^m \frac{1}{m!} \left( \frac{\partial }{\partial \alpha } \right) ^q \textbf{M}\ \textbf{e}_{0,m-q} (\alpha ,\theta )_{|\alpha =\alpha _0} = 0 \quad \forall m=1,2,\dots , k. \end{aligned}$$

We call the number $k+1$ the length of the Jordan chain.

The main result of this section was proven, e.g., in [30, Theorem 3.1]; see also the references therein.

Theorem A.1

Let $\textbf{u}$ be the solution to (3) on a polygonal domain $\Omega$. Given $\textbf{A}$ one of the vertices of $\Omega$, $\textbf{u}$ admits an asymptotic expansion in polar coordinates at $\textbf{A}$ of the form

$$\begin{aligned} \textbf{u}(r,\theta ) = \eta (r) \left[ \sum _{i=1}^M \sum _{j=0}^{m_i-1} \textbf{c}_{i,j} \textbf{s}_{i,j} (r,\theta ) \right] + \textbf{w}(r,\theta ). \end{aligned}$$

In the above equation:

• $\eta$ denotes a cut-off function, which localises the singular behaviour of the solution;

• M is the number of eigenvalues of $\textbf{M}$ with $\mathbb{R}\mathbb{e} (\alpha ) \in (0,1)$, where multiple eigenvalues are counted multiple times accordingly;

• $m_i$, for all $i=1,\dots ,M$, is the length of the Jordan chain corresponding to the eigenfunction $\textbf{e}{i}$;

• $\textbf{c}_{i,j}$ are suitable coefficients, called stress intensity factors, which are explicitly known [27];

• $\textbf{s}_{i,j}$ are the singular functions

  $$\begin{aligned} \textbf{s}_{i,j} (r,\omega ) = r^{\alpha _i} \sum _{k=0}^j \frac{\ln ^k(r)}{k!} \textbf{e}_{i,j-k} (\alpha _i,\theta ); \end{aligned}$$

  (40)

• $\textbf{w}$ belongs to $H^2$ in a neighbourhood inside $\Omega$ of $\textbf{A}$.

Remark 2

As discussed, e.g., in [14, Sect. 6], if the multiplicity of $\alpha$ is equal to the kernel of the system matrix in (34), then the length of the Jordan chain is one. Consequently, expansion (40) contains no logarithmic factors. This is for instance the case of the L-shaped domain benchmark in Example A.1. In principle, as for the slit domain benchmark in Example A.2, the multiplicity of all the $\alpha$ solving the characteristic equation is 2, which is the dimension of the kernel of the matrix (37). For this reason, we could also compute the second element of the Jordan chain and investigate the behaviour of singular functions containing an extra logarithmic singularity factor.