Plane crack problems within strain gradient elasticity and mixed finite element implementation

Abstract

An alternative approach is proposed and applied to solve boundary value problems within the strain gradient elasticity theory. A mixed variation formulation of the finite element method (FEM) based on the concept of the Galerkin method is used. To construct finite-dimensional subspaces separate approximations of displacements, deformations, stresses, and their gradients are implemented by choosing the different sets of piecewise polynomial basis functions, interrelated by the stability condition of the mixed FEM approximation. This significantly simplifies the pre-requirement for approximating functions to belong to class C¹ and allows one to use the simplest triangular finite elements with a linear approximation of displacements under uniform or near-uniform triangulation conditions. Global unknowns in a discrete problem are nodal displacements, while the strains and stresses and their gradients are treated as local unknowns. The conditions of existence, uniqueness, and continuous dependence of the solution on the problem’s initial data are formulated for discrete equations of mixed FEM. These are solved by a modified iteration procedure, where the global stiffness matrix for classical elasticity problems is treated as a preconditioning matrix with fictitious elastic moduli. This avoids the need to form a global stiffness matrix for the problem of strain gradient elasticity since it is enough to calculate only the residual vector in the current approximation. A set of modeling plane crack problems is solved. The obtained solutions agree with the results available in the relevant literature. Good convergence is achieved by refining the mesh for all scale parameters. All three problems under study exhibit specific qualitative features characterizing strain gradient solutions namely crack stiffness increase with length scale parameter and cusp-like closure effect.

1 Introduction

Gradient elasticity theories, a natural generalization of classical elasticity, allow accounting for the microstructure of actual materials. In particular, it is important in problems where classical elasticity has limitations. These limitations arise when we encounter size-dependent phenomena (e.g., homogenization problems, accounting for surface or interface energies) or discontinuity in the boundary conditions (edges, corners, or concentrated forces). In such cases, to capture the size effect or avoid singularity in solution, it is necessary to consider the material's microstructure. One such generalized continuum theory that has received considerable attention in recent years is the so-called strain gradient elasticity introduced by [1,2,3].

The key point of practical applications of any gradient theory is the specification of nonclassical elastic parameters. It is evident that even using the symmetry properties associated with the reversibility of the deformations and the conjugation of tangential stresses, the theory requires such a huge number of experiments for determining the nonclassical moduli that it will be practically inapplicable. Therefore, the crucial moment of research within any gradient theory shifts to searching for specific models that would minimize the number of nonclassical moduli while preserving the possibility of explaining the most significant size effects, removing the singularities and boundary discontinuities. By introducing certain assumptions for higher-order constitutive relations, one can reduce the number of material constants of the gradient theories to obtain models with one [4,5,6,7,8,9,10,11], two [3, 6, 12], or three additional scale parameters [6, 12], so that the identification procedure is simplified.

Most results have dealt with the simplified strain gradient model [4, 13]. In this model, a sixth-rank tensor ${\mathbb{C}}_{6}$ is taken in the form ${{\mathbb{C}}_{6}=l}^{2}{\mathbb{C}}_{4}\otimes {\mathbf{I}}_{2}$, where ${\mathbb{C}}_{4}$ is the classical stiffness tensor, ${\mathbf{I}}_{2}$ is the unit second-rank tensor, and the characteristic length $l$ is the only nonclassical parameter describing scale effects. Additionally, the model has been extensively used in the last two decades to solve various boundary value problems (see, e.g., [14,15,16,17,18,19,20,21,22,23,24,25]), and it will be employed in this article to solve the boundary problem of strain gradient elasticity theory.

Various authors have adopted this model and its variational counterpart for studying dislocation, disclination, and crack problems [5, 6, 24, 26,27,28,29,30,31,32]. Within the context of Toupin-Mindlin strain gradient elasticity, analytical studies for plane-strain crack problems were carried out by [14, 33]. Numerical studies employing finite and boundary elements (see, e.g., [21, 28, 31, 34,35,36,37] among other) have been developed to address problems within the Toupin-Mindlin theory.

In the present work, an alternative mixed variation formulation of FEM is proposed. This approach is based on the concept of the Galerkin method, and unlike the classical formulation of FEM in the form of a displacement method, displacements, strains, stresses, and their gradients are direct arguments of the variational equations of the finite-dimensional problem, and are not determined by further transformation (differentiation) of the problem solution in displacements. To construct finite-dimensional subspaces based on such variational formulation of the problem, separate (formally independent) approximations of displacements, deformations, stresses, and their gradients are used by choosing the different sets of piecewise polynomial basis functions, interconnected by the stability condition of the mixed approximation [38, 39].

The presentation is organized as follows. In the next section, we introduce a brief outline of the Toupin-Mindlin strain gradient elasticity theory, the simplifications within this theory that are employed in the solution, and discuss the physical meaning and estimation of the length scale parameter. Section 3 contains a variational formulation of the boundary value problem and condition of solvability of the discrete equation of mixed FEM, a description of the solution method, and numerical results. The results are discussed in the context of the similar ones available in the literature. The last section presents conclusions.

2 A brief outline of strain gradient elasticity theory

2.1 Toupin-Mindlin strain gradient elasticity

Within the first strain gradient theory, the strain and strain gradient energy is defined as follows:

$$w\left({ \mathbf{E}}_{2},{\mathbb{E}}_{3} \right)=\frac{1}{2}{ \mathbf{E}}_{2}:\,{\mathbb{C}}_{4}:\,{ \mathbf{E}}_{2}+{ \mathbf{E}}_{2}:\,{\mathbb{C}}_{5}\vdots {\mathbb{E}}_{3}+\frac{1}{2}{\mathbb{E}}_{3}\vdots {\mathbb{C}}_{6}\vdots {\mathbb{E}}_{3},$$

(1)

where ${\mathbb{C}}_{4}, {\mathbb{C}}_{5}, {\mathbb{C}}_{6}$ are stifnesses, while ${\mathbf{E}}_{2}$ and ${\mathbb{E}}_{3}$ are the strains and the second gradient of displacement vector $\mathbf{u}(\mathbf{x})$:

$${\mathbf{E}}_{2}={\text{sym}}\left(\mathbf{u}\otimes \nabla \right),{\mathbb{E}}_{3}=\mathbf{u}\otimes \nabla \otimes \nabla .$$

(2)

The following notations are utilized hereinafter: scalars, vectors, second- and higher-rank tensors are denoted by italic letters (such as $\mathbf{e}$ or $\mathbf{E}$), bold minuscules (such as $\mathbf{e}$), bold majuscules (such as $\mathbf{E}$), and bold blackboard majuscules (such as ${\mathbb{E}}$), respectively. The indices indicate the tensor rank of ${\mathbf{E}}_{2}, {\mathbb{E}}_{3},$ ${\mathbb{C}}_{4}, {\mathbb{C}}_{5}, {\text{and}} {\mathbb{C}}_{6}$. These tensors possess the following symmetries:

$$ \begin{aligned} & E_{ij} = E_{ji} , \\ & E_{ijk} = E_{ikj} , \\ & C_{ijkl} = C_{klij} = C_{jikl} = C_{ijlk} , \\ & C_{ijklm} = C_{jiklm} = C_{ijkml} , \\ & C_{ijklmn} = C_{lmnijk} = C_{ikjlmn} = C_{ijklnm} . \end{aligned} $$

(3)

“$\nabla$” is the three-dimensional nabla operator.

$$\nabla \equiv \frac{\partial }{\partial {x}_{{\varvec{i}}}}{\mathbf{e}}_{i},\left(i=1, 2, 3\right),$$

(4)

with the orthonormal base vectors ${\mathbf{e}}_{1}\boldsymbol{ },\boldsymbol{ }{\mathbf{e}}_{2}\boldsymbol{ },\boldsymbol{ }{\mathbf{e}}_{3}$. Here, a summation over repeated indices is implied. “$\otimes$” denotes the dyadic product, and the displacement gradient is defined as

$$\mathbf{u}\otimes \nabla \equiv \frac{\partial {u}_{{\varvec{i}}}}{\partial {x}_{{\varvec{j}}}}{\mathbf{e}}_{i}\otimes {\mathbf{e}}_{j}={u}_{i,j}{\mathbf{e}}_{i}\otimes {\mathbf{e}}_{j}\boldsymbol{ }.$$

(5)

The dots denote scalar contractions, where the double and triple scalar contractions are defined concerning orthonormal base vectors ${\mathbf{e}}_{i}$ as follows:

$$ \begin{aligned} & {\mathbf{e}}_{i} \otimes {\mathbf{e}}_{j} :{\mathbf{e}}_{k} \otimes {\mathbf{e}}_{l} = \delta_{ik} \delta_{jl} , \\ & {\mathbf{e}}_{i} \otimes {\mathbf{e}}_{j} \otimes {\mathbf{e}}_{m} \vdots {\mathbf{e}}_{k} \otimes {\mathbf{e}}_{l} \otimes {\mathbf{e}}_{n} = \delta_{ik} \delta_{jl} \delta_{mn} .\\ \end{aligned} $$

(6)

The stresses and the double stresses are given by

$$ \begin{aligned} {\mathbf{T}}_{2} & = \frac{\partial w}{{\partial {\mathbf{E}}_{2} }} = {\mathbb{C}}_{4} :\,{\mathbf{E}}_{2} + {\mathbb{C}}_{5} \vdots {\mathbb{E}}_{3} , \\ {\mathbb{T}}_{3} & = \frac{\partial w}{{\partial {\mathbb{E}}_{3} }} = {\mathbb{C}}_{5}^{{\text{T}}} :\,{\mathbf{E}}_{2} + {\mathbb{C}}_{6} \vdots {\mathbb{E}}_{3} . \\ \end{aligned} $$

(7)

2.2 Simplified model

The stiffness tensors ${\mathbb{C}}_{4}$, ${\mathbb{C}}_{5}$, and ${\mathbb{C}}_{6}$in the case of hemitropy (the absence of improper rotations preserves the existence of ${\mathbb{C}}_{5}$,which is the only difference from isotropic strain gradient elasticity) are characterized by eight independent parameters and have the following form (cf. [2, 40, 41]).

$$ \begin{aligned} {\mathbb{C}}_{4} & = \left[{\lambda \delta_{ij} \delta_{kl} + \mu \left( {\delta_{ik}\delta_{jl} + \delta_{il} \delta_{kj} } \right)} \right]\user2{}{\mathbf{e}}_{i} \otimes {\mathbf{e}}_{j} \otimes {\mathbf{e}}_{k}\otimes {\mathbf{e}}_{l} , \\ {\mathbb{C}}_{5} & = \left[ {\kappa \left( {\varepsilon_{imk} \delta_{jl} + \varepsilon_{ilk}\delta_{jm} + \varepsilon_{jmk} \delta_{il} + \varepsilon_{jlk}\delta_{im} } \right)} \right]\user2{ }{\mathbf{e}}_{i} \\ & \quad\otimes {\mathbf{e}}_{j} \otimes {\mathbf{e}}_{k} \otimes {\mathbf{e}}_{l}\otimes {\mathbf{e}}_{m} , \\ {\mathbb{C}}_{6} & = [c_{1} \left({\delta_{jk} \delta_{im} \delta_{nl} + \delta_{jk} \delta_{in}\delta_{ml} + \delta_{ij} \delta_{kl} \delta_{mn} + \delta_{jl}\delta_{ik} \delta_{mn} } \right) \\ & \quad + c_{2} \left({\delta_{ij} \delta_{km} \delta_{nl} + \delta_{jm} \delta_{ki}\delta_{nl} + \delta_{ij} \delta_{kn} \delta_{ml} + \delta_{jn}\delta_{ik} \delta_{ml} } \right) \\ & \quad + c_{3} \left({\delta_{jm} \delta_{kl} \delta_{in} + \delta_{jl} \delta_{in}\delta_{km} + \delta_{jn} \delta_{im} \delta_{kl} + \delta_{jl}\delta_{im} \delta_{nk} } \right) \\ & \quad + c_{4} \left({\delta_{jn} \delta_{il} \delta_{km} + \delta_{jm} \delta_{kn}\delta_{il} } \right) + c_{5} \delta_{il} \delta_{jk} \delta_{mn}]\user2{ }{\mathbf{e}}_{i} \otimes {\mathbf{e}}_{j}\\ & \quad \otimes {\mathbf{e}}_{k} \otimes {\mathbf{e}}_{l} \otimes {\mathbf{e}}_{m}\otimes {\mathbf{e}}_{n} , \\ \end{aligned}$$

(8)

where ${\delta }_{ij}$ is the Kronecker symbol, ${\varepsilon }_{imk}$ is the Levi‒Civita permutation symbol, $\lambda$ and $\mu$ are standard Lamé’s coefficients, and $\kappa$ and ${c}_{1}-{c}_{5}$ are the higher order material parameters.

The development of physically simplified models and reducing the number of material parameters are necessary to solve the practical problems within the strain gradient elasticity. One of the possibilities is a simplified model of [4] and [13] adopted here to solve some crack problems. It is assumed that

$$\kappa =0, {c}_{1}=0, {c}_{2}=0, {c}_{3}=\frac{\lambda }{2}{l}^{2} , {c}_{4}={\mu }{l}^{2}, {c}_{5}= 0$$

(9)

where $l$ has the length dimension, and $\lambda$ and $\mu$ are Lamé parameters.

Positive definiteness of the strain and strain gradient energy requires the following inequality constraints for material parameters [3, 42]:

$$3\lambda +2\mu >0, \mu >0, {l}^{2}>0.$$

(10)

In this case, the strain and strain gradient energy can be written as follows [42]:

$$w=\frac{1}{2}{ \mathbf{T}}_{2}: { \mathbf{E}}_{2}+\frac{1}{2}{l}^{2}{\mathbf{T}}_{2}\otimes \nabla\, \vdots \,{ \mathbf{E}}_{2}\otimes \nabla ,$$

(11)

The elastic energy is symmetric concerning the strain and the stress as well as the strain gradient and the stress gradient

$$\frac{\partial w}{\partial {(\mathbf{E}}_{2}\otimes \nabla )}={l}^{2}{\mathbf{T}}_{2}\otimes \nabla , \frac{\partial w}{\partial {(\mathbf{T}}_{2}\otimes \nabla )}={l}^{2}{\mathbf{E}}_{2}\otimes \nabla .$$

(12)

Then, the double stresses are

$${\mathbb{T}}_{3}={l}^{2}{\mathbf{T}}_{2}\otimes \nabla .$$

(13)

The equations of equilibrium take the following form [13]:

$$(1-{l}^{2}{\nabla }^{2}){(\nabla \cdot \mathbf{T}}_{2})=0.$$

(14)

The boundary conditions for prescribed traction and prescribed double traction in the normal direction on the part of the boundary ${{\text{S}}}_{\tau }$:

$$ \begin{aligned} {\mathbf{p}}_{pr} & =\left( {1 - l^{2} \nabla^{2} - l^{2} \nabla \cdot \nabla_{{\varvec{s}}} } \right)({\mathbf{T}}_{2} \cdot {\mathbf{n}})\\&\quad+ l^{2} {\mathbf{T}}_{2} :\left[ {\left( {\nabla_{{\varvec{s}}}\cdot {\mathbf{n}}} \right){\mathbf{n}} \otimes {\mathbf{n}} -{\mathbf{n}} \otimes \nabla_{{\varvec{s}}} } \right], \\{\mathbf{r}}_{n pr} & = l^{2} {\mathbf{T}}_{2} :{\mathbf{n}} \otimes {\mathbf{n}}, \\ \end{aligned} $$

(15)

and the kinematic boundary conditions in terms of the prescribed displacement and their normal derivative on the part of the body surface ${{\text{S}}}_{u}$:

$$ \begin{aligned} & {\mathbf{u}}_{pr} = {\mathbf{u}}, \\ & {\mathbf{Du}}_{pr} = {\mathbf{Du}} \\ \end{aligned} $$

(16)

Here “${\nabla }^{2}$” denotes Laplacian, ${{\text{S}}}_{\tau }\cup {{\text{S}}}_{u}={\text{S}}$, and $\mathbf{D}\mathbf{u}=\left(\mathbf{u}\otimes \nabla \right)\cdot \mathbf{n}= \frac{\partial {u}_{i}}{\partial {x}_{j}}{n}_{j}{e}_{i}$.

2.3 Identification and estimation of the length scale parameter

To clarify the physical meaning and to estimate the length scale parameter in the exploited models within strain gradient elasticity, namely, scaling parameters (parameter $l$ in Eqs. 9–15), the strategy, which originates from continualization (homogenization) of a discrete microstructure (see, e.g., [43,44,45,46,47]), is considered. The key idea is to consider interactions between adjacent particles and links between next-to-adjacent particles. In this case, the associated continuous formulation is expected to be multiscale, with the first scale greater than the distance between the micro-inhomogeneity and the second reflecting the peculiarity of nonlocal interactions compared to local interactions. The latter will play the role of an intermediate scale, for which a typical length scale parameter is still greater than the distance between the micro-inhomogeneity but less than that for the first (i.e., overall) long scale. As a result, it leads to the derivation of strain gradient constitutive relations and potential energy of the strain gradient continuum as a function of strain and strain gradient (e.g., [45, 48]). It might be expected that, along with the overall stiffnesses, the length scale parameters are also dependent on the size of the micro-inhomogeneity and the distances between them.

3 Finite element implementation of plane crack problems

The specific feature of solving the variational equations of the gradient elasticity theory consists of accounting for the first derivatives from the strain tensor components (second derivatives from displacements) [37, 49,50,51,52 among others]. A precondition for the convergence of finite element method solutions (FEM) is the property of approximating functions to be continuously differentiable. In other words, approximating functions must belong to a class of smooth C¹ functions [53, 54]. This leads to significant mathematical and computational difficulties since the dimension of “local” spaces of finite elements increases, and their structure becomes essentially complicated [e.g., 54, 55, among others]. For example, the triangular finite element introduced by Zienkiewicz for the thin plate bending problem satisfies the condition of continuity of the normal derivative on the triangle sides just if the triangle mesh is formed by a system of three equidistant parallel lines [53, 54].

To take into consideration the condition of continuity of the first derivatives from displacements, special finite elements (e.g., Argyris triangle, Bogner-Fox-Schmit rectangle, Zienkiewicz singular triangle, among others) are used in which the nodal unknowns are the displacements, their first and even second derivatives (see, e.g., [37, 49, 54, among others]). Using elements of this type significantly increases the matrix bandwidth and the dimension of the global stiffness matrix, which complicates the computational process.

In the present work, an alternative approach is developed in which a mixed variation formulation of FEM concerning displacement‒strain–stress, and their gradients is used to solve the boundary problem of the gradient elasticity theory.

Here, the classical formulation of FEM in displacements is changed so that displacements, strains, stresses, and their gradients are the direct arguments of the finite-dimensional problem and not are determined from the solution of the problem in displacements. It significantly simplifies the selection of approximating functions since it is not necessary to use high-order finite elements having continuous first derivatives of displacements. The construction of finite-dimensional spaces is based on the separate approximations of the fields of displacements, deformations, stresses and their gradients using a different set of the piecewise polynomial basis functions, interconnected by the stability condition of the mixed approximation [38, 39].

The original premises for mixed FEM formulations may be different. One of them is based on the use of the variational principles of mechanics, according to which the solution of the boundary problem is reduced to finding a stationary value of some variational functional with respect to the corresponding arguments [56,57,58]. The second is based on the methods of the duality theory and on the dual-basic formulation of the problem under assumption of the duality principle [53, 56, 59]. The more general approach is reached from an abstract formulation of the boundary value problem using the results of the theory developed by Brezzi [60]. Within this approach, a mixed formulation of boundary value problem is considered, regardless of whether it is obtained from the problem of the saddle point or not.

A characteristic feature of the mixed FEM formulations of the problems of gradient elasticity theory is that the displacements, deformations, stresses and their gradients are determined simultaneously in the process of solving a finite-dimensional problem, in contrast to the classical FEM, in which they are calculated from the results of the solution in displacements. The mixed variation FEM formulation is based on the concept of the Galerkin method, according to which the equations of the boundary problem are multiplied by arbitrary continuous functions, which can be interpreted as variations in displacements, strains-stresses, and their gradients. Integration of the resulting relations over the region occupied by the body leads to integral identities that include partial derivatives from displacements of the first order only [38, 39]. This is a main feature of the formulated variational equations since the variational formulation in displacements using the Lagrange equations assumes a twofold differentiation of displacements, and the classical formulation contains partial derivatives from displacements up to the fourth order.

Using mixed FEM significantly simplifies the selection of approximating functions while eliminating the request to employ finite elements of the C¹ class. Thus, it is possible to use the simplest triangular finite elements with a linear approximation of displacements, provided that uniform or near-uniform triangulation is utilized. Uniform triangulation in the vicinity of the stress concentrators (crack mouth, concentrated forces, reinforcements, etc.) provides super-convergence in calculating deformations and stresses in mesh nodes [38]. Global unknowns in a discrete problem are nodal displacements, while the strains and stresses at mesh nodes are considered as local unknowns due to the use of special numerical integration formulas of the interpolation type. Specifically, such a formulation of FEM, based on a separate approximation of displacements-strains-stresses, and their gradients is used in this work.

In this regard, it is necessary to mention the article of Amanatidou and Aravas [28], in which three types of equations of the strain gradient elasticity theory are examined in detail, alternative mixed variational formulations of FEM were developed on their basis, and special rectangular finite elements of the mixed type of the class C⁰ were constructed. In these elements, both displacements and their gradients are used as independent unknowns, and their relationship is provided in the “integral sense.” A patch test in a mixed formulation is used to justify the correctness of the constructed finite elements, a necessary condition for the solvability of the discrete problem [54]. At the same time, this condition is not sufficient for stability and convergence of the mixed approximation since it is still necessary to ensure the fulfillment of the stability condition of Brezzi [60], which is a sufficient condition for a unique solution of the saddle point problem in which the solution continuously depends on the input data. Otherwise, discretization based on mixed approximation can be unstable and cause parasitic solution fluctuations.

In the present work, an alternative condition of the solvability and stability of the mixed approximation is formulated, equivalent to the Brezzi condition (inequality (38) presented below). This condition is obtained based on the concretization of the equations of the mixed method for strain gradient elasticity problems. Triangular finite elements of mixed type are constructed, providing the uniqueness of the solution and stability of the mixed approximation.

A modified iteration procedure with a preconditioning matrix is implemented to solve the system of discrete equations of mixed FEM [38]. A global stiffness matrix of the classical FEM is used as a preconditioning matrix for piecewise linear approximation of displacements on triangular elements with fictitious elastic moduli. The elastic moduli in this matrix are modified to provide convergence of the iterative process. It is also possible to use the method of conjugate gradients with a preconditioning matrix, which has a higher convergence rate [61, 62].

In the above methods, there is no need to form a global stiffness matrix of the mixed FEM, to perform each iteration, it is enough to calculate only the residual vector in the current approximation. This significantly simplifies the solution of the equations of the mixed method and reduces requirements regarding computer memory.

An alternative to the mixed FEM approach for solving the strain gradient elasticity theory equations is the hybrid FEM approach that provides the continuity of the first displacement derivatives [54]. One of such finite elements is a hybrid finite element constructed based on the Zienkiewicz triangle [63], which guarantees solution convergence with mesh clustering regardless of the method of partitioning into triangular elements. The results obtained by this approach will be presented in the next publication of the authors.

3.1 Variational formulation of the boundary value problems

Let us consider a body occupying the area $\Omega$ with a regular boundary S. Displacements are defined on the part of the boundary ${{\text{S}}}_{u}$, and surface traction ${f}_{p}$ is defined on the rest part ${{\text{S}}}_{\tau }$ (${{\text{S}}}_{\tau }\cup$ ${{\text{S}}}_{u}={\text{S}}$). In addition, the body is subjected to volume forces ${F}_{p}$. It is assumed that the displacements satisfy boundary conditions on the part of the body surface, excluding its rigid motion.

Let $W$ be a set of admissible vector functions for displacements

$$W=\left(\mathbf{v}|\mathbf{v}=\left({v}_{p}\right),\boldsymbol{ }\boldsymbol{ }1\le p\le n,\boldsymbol{ }{v}_{p}\in {H}^{2}\left(\Omega \right),{ v}_{p}=0\,on\,{{\text{S}}}_{u}\right),$$

(17)

where ${H}^{2}\left(\Omega \right)$—space of functions, square-integrable on $\Omega$ including its first and second derivatives; $n$ denotes either two- or three-dimensions.

Then, the boundary value problem of linear gradient elasticity is formulated as a variational equation w.r.t. displacement $\mathbf{u}\in W$ [52]:

$${\int }_{\Omega }\left({\mathbf{T}}_{2}\left(\mathbf{u}\right): {\mathbf{E}}_{2}\left(\mathbf{v}\right)+{\mathbb{T}}_{3}\left(\mathbf{u}\right)\vdots {\mathbb{E}}_{3}\left(\mathbf{v}\right)\right)d\Omega ={\int }_{\Omega }\mathbf{F}\cdot \mathbf{v}d\Omega +{\int }_{{{\text{S}}}_{u}}\mathbf{f}\cdot \mathbf{v}d{{\text{S}}}_{u}, \forall \mathbf{v}\in W.$$

(18)

An alternative variational formulation of the problem. According to (2), (8), (9), (12), and (13), let us write down the boundary value problem of the gradient elasticity theory by a system of integral identities:

$$ \begin{aligned} & \mathop \int \limits_{{\Omega }} {\mathbf{E}}_{2} :\delta {\mathbf{T}}_{2} d{\Omega } = \mathop \int \limits_{{\Omega }} {\text{sym}}\left( {{\mathbf{u}} \otimes \nabla } \right):\delta {\mathbf{T}}_{2} d{\Omega ,} \\ & \mathop \int \limits_{{\Omega }} {\mathbf{T}}_{2} :\delta {\mathbf{E}}_{2} d{\Omega } = \mathop \int \limits_{{\Omega }} \left( {\lambda {\mathbf{I}}_{2} {\text{tr}}\left( {{\mathbf{E}}_{2} } \right) + 2\mu {\mathbf{E}}_{2} } \right):\,\delta {\mathbf{E}}_{2} d{\Omega }, \\ & \mathop \int \limits_{{\Omega }} {\mathbb{E}}_{3} \vdots \delta {\mathbb{T}}_{3} d{\Omega } = \mathop \int \limits_{{\Omega }} {\mathbf{E}}_{2} \otimes \nabla \vdots \delta {\mathbb{T}}_{3} d{\Omega }, \\ & \mathop \int \limits_{{\Omega }} {\mathbb{T}}_{3} \vdots \delta {\mathbb{E}}_{3} d{\Omega } = l^{2} \mathop \int \limits_{{\Omega }} {\mathbf{T}}_{2} \otimes \nabla\, \vdots\, \delta {\mathbb{E}}_{3} d{\Omega }, \\ & \mathop \int \limits_{{\Omega }} \left( {{\mathbf{T}}_{2} :\,\delta {\mathbf{E}}_{2} + {\mathbb{T}}_{3}\, \vdots \, \delta {\mathbb{E}}_{3} } \right)d{\Omega } = \mathop \int \limits_{{\Omega }} {\mathbf{F}} \cdot \delta {\mathbf{u}}d{\Omega } + \mathop \int \limits_{{{\text{S}}_{u} }} {\mathbf{f}} \cdot \delta {\mathbf{u}}d{\text{S}}_{u} \\ \end{aligned} $$

(19)

Using the variational formulation of the boundary value problem in the form of Eqs. (19) makes it possible to weaken requirements for the smoothness of approximating functions since it is enough to consider just continuous piecewise-polynomial functions to construct an approximate FEM solution of the problem, which essentially simplifies it.

Remark 3.1

The effective tools for analysis of the conditions of solvability, stability and convergence of approximate solutions of boundary problems are functional analysis theory and general theory of linear operator equations in Hilbert spaces. Formulation of variational equations of a mixed approach using linear operators is, apparently, the most convenient for the study of these conditions. It is a reason why the variational FEM formulation of the boundary value problem will be presented in a compact and easy-to-use and to analysis equivalent form utilizing linear operators in Hilbert spaces.

The set of admissible vector functions for displacements will be considered as elements of the functional space

$$U=\left(\mathbf{v}|\mathbf{v}=\left({v}_{p}\right),\boldsymbol{ }\boldsymbol{ }1\le p\le n,\boldsymbol{ }{v}_{p}\in {H}^{1}\left(\Omega \right),{ v}_{p}=0\,on\,{{\text{S}}}_{u} \right),$$

(20)

where ${H}^{1}\left(\Omega \right)$—space of functions, square-integrable on $\Omega$ including its first derivatives. Following Korn's inequality [53], the set $U$ is complete concerning the norm associated with the scalar product

$${\left(\mathbf{u},\mathbf{v}\right)}_{U}={\int }_{\Omega }{\mathbf{E}}_{2}\left(\mathbf{u}\right):{\mathbf{E}}_{2}\left(\mathbf{v}\right)d\Omega , \forall \mathbf{u},\mathbf{v}\in U.$$

(21)

Therefore, $U$ is a Hilbert space with the norm $|| \mathbf{v}{ ||}_{U}={\left(\mathbf{v},\mathbf{v}\right)}_{U}^{1/2}$, $\forall \mathbf{v}\in U$, which is equivalent to the norm of the space ${\left({H}^{1}\left(\Omega \right)\right)}^{n}.$

The range of permissible values for strains ${\mathbf{E}}_{2}$ and stresses ${\mathbf{T}}_{2}$ is defined as the set of symmetric tensor functions of the second rank, all components of which belong to the space ${H}^{1}\left(\Omega \right)$:

$$X{=}\left({\mathbf{G}}_{2}|{\mathbf{G}}_{2}{=}\left({G}_{pq}\right), 1{\le} p,q{\le} n , { G}_{pq}{=}{G}_{qp},{ G}_{pq}\in {H}^{1}\left(\Omega \right)\right).$$

(22)

As the range of permissible values for the strain gradient ${\mathbb{E}}_{3}$ and double stresses ${\mathbb{T}}_{3}$, a set of third-rank tensor functions, the components of which belong to ${L}_{2}\left(\Omega \right)$, will be considered:

$$Z{=}\left({\mathbb{M}}_{3}\boldsymbol{ }|\boldsymbol{ }{\mathbb{M}}_{3}{=}\left({M}_{rpq}\right),\boldsymbol{ }{ M}_{rpq}{=}{M}_{rqp}, 1{\le} r,p,q{\le} n, { M}_{rpq}\in {L}_{2}\left(\Omega \right)\right),$$

(23)

where ${L}_{2}\left(\Omega \right)$ is the space of square-integrable functions on $\Omega$.

Let us introduce Hilbert spaces L and H, consisting of the sets of strains-stresses tensor functions square-integrable on $\Omega$ and their gradients, respectively. The scalar products are defined in spaces L and H as follows:

$${\left({\mathbf{T}}_{2},{ \mathbf{E}}_{2}\right)}_{L}={\int }_{\Omega }{\mathbf{T}}_{2}:{\mathbf{E}}_{2}d\Omega , {\left({\mathbb{T}}_{3},{\mathbb{E}}_{3}\right)}_{H}={\int }_{\Omega }{\mathbb{T}}_{3} \ \vdots \ {\mathbb{E}}_{3}d\Omega .$$

(24)

Then, the boundary value problem can be formulated in the operator form as follows. To find displacements $\mathbf{u}$ strains ${\mathbf{E}}_{2}$ and stresses ${\mathbf{T}}_{2}$ $\left(\mathbf{u},{ \mathbf{E}}_{2},{\mathbf{T}}_{2}\right)\in U\times X\times X$ and stress gradient ${\mathbb{E}}_{3}$ and double stresses ${\mathbb{T}}_{3}$ $\left({\mathbb{E}}_{3},{ {\mathbb{T}}}_{3}\right)\in Z\times Z$ such that

$$ \begin{aligned} & \left( {{\mathbf{E}}_{2} ,\delta {\mathbf{T}}_{2} } \right)_{L} = \left( {{\text{B}}{\mathbf{u}},\delta {\mathbf{T}}_{2} } \right)_{L} ,~~\forall ~ \delta {\mathbf{T}}_{2} \in X, \\ & \left( {{\mathbf{T}}_{2} ,\delta {\mathbf{E}}_{2} } \right)_{L} = \left( {{\text{D}}{\mathbf{E}}_{2} ,\delta {\mathbf{E}}_{2} } \right)_{L} ,~~\forall ~\delta {\mathbf{E}}_{2} \in X, \\ & \left( { {\mathbb{{E}}}_{3} ,\delta {\mathbb{{T}}}_{3} } \right)_{H} = \left( {{\text{N}}{\mathbf{E}}_{2} ,\delta {\mathbb{{T}}}_{3} } \right)_{H} ,~~\forall ~\delta {\mathbb{ {T}}}_{3} \in Z, \\ & \left( { {\mathbb{{T}}}_{3} ,\delta {\mathbb{{E}}}_{3} } \right)_{H} = l^{2} \left( {{\text{N}}{\mathbf{T}}_{2} ,\delta {\mathbb{{E}}}_{3} } \right)_{H} ,~~\forall ~\delta {\mathbb{ {E}}}_{3} \in Z, \\ & \left( {{\mathbf{T}}_{2} ,\delta {\mathbf{E}}_{2} } \right)_{L} + \left( { {\mathbb{ {T}}}_{3} ,\delta {\mathbb{{E}}}_{3} } \right)_{H} = \left\langle {{\mathbf{f}},\delta {\mathbf{u}}} \right\rangle ,~~\forall ~\delta {\mathbf{u}} \in U, \\ \end{aligned} $$

(25)

where ${\text{B}}$ is the linear differential operator for the determination of strains ${\mathbf{E}}_{2}$

$$ {\text{B}}{\mathbf{u}} = {\text{sym}}\left( {{\mathbf{u}} \otimes \nabla } \right). $$

(26)

${\text{D}}$ is the linear self-adjoined positive definite operator associated with a symmetric tensor of material elastic moduli

$${\text{D}}{\mathbf{E}}_{2}={\mathbb{C}}_{4}:{\mathbf{E}}_{2}.$$

(27)

${\text{N}}$ is the linear differential operator for calculating strain and stress gradients

$${\text{N}}{\mathbf{E}}_{2}={ \mathbf{E}}_{2}\otimes \nabla \;\text{or N}{\mathbf{T}}_{2}={ \mathbf{T}}_{2}\otimes \nabla .$$

(28)

$\langle \cdot ,\cdot \rangle$ is the linear form, identified with the work of applied loads to the body on displacements.

Let us set $Y$ to be defined as the value domain of operator ${\text{B}}$

$${\text{B}}\mathbf{u}: U\to Y,\forall \mathbf{u}\in U,$$

(29)

and $Y$ is considered as a linear subspace of the Hilbert space L with a scalar product ${\left(\cdot ,\cdot \right)}_{L}$. Sets X and Z are the definition and value domains of operator ${\text{N}}$

$${\text{N}}{\mathbf{E}}_{2} : X\to Z,\forall {\mathbf{E}}_{2}\in X.$$

(30)

$X$ and $Z$ are linear subspaces of the Hilbert spaces L and $H$, respectively. Given that $\delta {\mathbf{E}}_{2}={\text{B}}\delta \mathbf{u}$, $\delta {\mathbf{T}}_{2}={\text{DB}}\delta \mathbf{u}$, $\delta {\mathbb{E}}_{3}={\text{NB}}\delta \mathbf{u}$, and $\delta {\mathbb{T}}_{3}={l}^{2}{\text{NDB}}\delta \mathbf{u}$ (see (2), (8), (9), (12), (13), (29), and (30)), the variation formulation of the boundary value problem in displacements follows from (25):

$${\left({\text{DB}}\mathbf{u},{\text{B}}\delta \mathbf{u}\right)}_{L}+{l}^{2}{\left({\text{NDB}}\mathbf{u} ,{\text{NB}}\delta \mathbf{u}\right)}_{H}=\langle \mathbf{f},\delta \mathbf{u}\rangle ,\forall \delta \mathbf{u}\in W.$$

(31)

For the continuum problems of gradient elasticity theory, the variational formulation in displacement in (25) and (31) are equivalent under the condition of $\mathbf{u}\in W$. Still, the latter is more flexible for constructing approximate solutions based on mixed FEM schemes.

Finite-dimensional problem formulation. Let ${U}_{h}$, ${X}_{h}$, ${Z}_{h}$ be finite-dimensional approximating spaces satisfying the following conditions:

$${U}_{h}\subset U, {X}_{h}\subset X, {Z}_{h}\subset Z,$$

(32)

where $h$ is the key parameter of the set of finite-dimensional spaces, which goes to zero in the limit. Then, the finite-dimensional problem is formulated similarly to (25). To find $\left({\mathbf{u}}_{h},{{\mathbf{E}}_{2}}_{h},{{\mathbf{T}}_{2}}_{h}\right)\in {U}_{h}\times {X}_{h}\times {X}_{h}$ and $\left({{\mathbb{E}}_{3}}_{h},{{\mathbb{T}}_{3}}_{h}\right)\in {Z}_{h}\times {Z}_{h}$ such that

$$ \begin{aligned} & \left( {{\mathbf{E}}_{{2h}} ,\delta {\mathbf{T}}_{{2h}} } \right)_{L} = \left( {{\text{B}}{\mathbf{u}}_{h} ,\delta {\mathbf{T}}_{{2h}} } \right)_{L} ,~~\forall ~\delta {\mathbf{T}}_{{2h}} \in X_{h} , \\ & \left( {{\mathbf{T}}_{{2h}} ,\delta {\mathbf{E}}_{{2h}} } \right)_{L} = \left( {{\text{D}}{\mathbf{E}}_{{2h}} ,\delta {\mathbf{E}}_{{2h}} } \right)_{L} ,~~\forall ~\delta {\mathbf{E}}_{{2h}} \in X_{h} , \\ & \left( { {\mathbb{{{E}}}}_{{3h}} ,\delta {\mathbb{{{T}}}}_{{3h}} } \right)_{H} = \left( {{\text{N}}{\mathbf{E}}_{{2h}} ,\delta {\mathbb{{{T}}}}_{{3h}} } \right)_{H} ,~~~\delta {\mathbb{{{T}}}}_{{3h}} \in Z_{h} , \\ & \left( { {\mathbb{{{T}}}}_{{3h}} ,\delta {\mathbb{{{E}}}}_{{3h}} } \right)_{H} = l^{2} \left( {{\text{N}}{\mathbf{T}}_{{2h}} ,\delta {\mathbb{{{E}}}}_{{3h}} } \right)_{H} ,~~~\delta {\mathbb{{{E}}}}_{{3h}} \in Z_{h} , \\ & \left( {{\mathbf{T}}_{{2h}} ,\delta {\mathbf{E}}_{{2h}} } \right)_{L} + \left( { {\mathbb{{{T}}}}_{{3h}} ,\delta {\mathbb{{{E}}}}_{{3h}} } \right)_{H} = \left\langle {\mathbf{f}},\delta {\mathbf{u}}_{h}\right\rangle ,~~\delta {\mathbf{u}}_{h} \in U_{h} . \\ \end{aligned} $$

(33)

System (33) is a variational formulation for the boundary value problem of the gradient elasticity theory in the form of displacement–strain-stress-gradients.

Solvability of a discrete problem solution. To analyze the properties of (33), it is necessary to establish a correspondence between subspaces ${X}_{h}\subset X$ and ${Y}_{h}\subset Y$. Note that ${Y}_{h}$ is a value domain of operator B acting on the closed subspace of ${U}_{h}$, i.e., ${Y}_{h}={\text{B}}{U}_{h}$. Therefore, ${Y}_{h}$ is an approximating subspace for space $Y$. Both ${X}_{h}$ and ${Y}_{h}$ are finite-dimensional subspaces of the Hilbert space L, but neither is in general a subset of the other.

Let us introduce an orthogonal projector operator ${{\text{I}}}_{h}$, corresponding to each element from space ${Y}_{h}$ and its orthogonal projection in ${X}_{h}$. Operator ${{\text{I}}}_{h}$ is associated with the scalar product ${\left(\cdot ,\cdot \right)}_{L}$ and defined according to equality

$${\left({\text{B}}{\mathbf{u}}_{h}-{{\text{I}}}_{h}{\text{B}}{\mathbf{u}}_{h},\delta {{\mathbf{T}}_{2}}_{h}\right)}_{L}=0, \forall \delta {{\mathbf{T}}_{2}}_{h}\in {X}_{h}.$$

(34)

Using the orthogonal projector ${{\text{I}}}_{h}: {Y}_{h}\to {X}_{h}$ set of Eqs. (33) can be written as a variational equation in displacements:

$${\left({{\text{DI}}}_{h}{\text{B}}{\mathbf{u}}_{h},{{\text{I}}}_{h}{\text{B}}\delta {\mathbf{u}}_{h}\right)}_{L}+{l}^{2}{\left({{\text{NDI}}}_{h}{\text{B}}{\mathbf{u}}_{h},{{\text{NI}}}_{h}{\text{B}}{\delta \mathbf{u}}_{h}\right)}_{H}=\langle \mathbf{f},{\delta \mathbf{u}}_{h}\rangle , \forall { \delta \mathbf{u}}_{h}\in {U}_{h}.$$

(35)

Let us introduce the linear operator ${\text{Q}}$ using the relation

$${\text{QN}}={\text{N}}\overline{{\text{D}} },$$

(36)

where $\overline{{\text{D}} }$—is a linear operator defined as follows. For a homogeneous and piecewise homogeneous body, operator $\overline{{\text{D}} }$ is identical to operator ${\text{D}}$. For an inhomogeneous body, variable elasticity modules are approximated by their mean values assigned to the centers or nodes of the finite element mesh. Thus, Q is a self-conjugate positive definite operator acting on the space $Z$. By utilizing the operator ${\text{Q}}$, Eq. (35) will have the following form:

$${\left({{\text{DI}}}_{h}{\text{B}}{\mathbf{u}}_{h},{{\text{I}}}_{h}{\text{B}}{\delta \mathbf{u}}_{h}\right)}_{L}+{l}^{2}{\left({{\text{QNI}}}_{h}{\text{B}}{\mathbf{u}}_{h},{{\text{NI}}}_{h}{\text{B}}\delta {\mathbf{u}}_{h}\right)}_{H}=\langle \mathbf{f},{\delta \mathbf{u}}_{h}\rangle , \forall { \delta \mathbf{u}}_{h}\in {U}_{h}.$$

(37)

Solvability condition. It is assumed that there is such a positive constant ${d}_{0}$ independent of $h$, for which the following inequality holds:

$${d}_{0}{|| {\text{B}}{\mathbf{v}}_{h}||}_{L}\le {|| {{\text{I}}}_{h}{\text{B}}{\mathbf{v}}_{h}||}_{L}, 0<{d}_{0}\le 1, \forall {\mathbf{v}}_{h}\in {U}_{h}.$$

(38)

Remark 3.2

The inequality (38) guarantees the unique existence of the solution for Eq. (33) at any parameter ${\text{h}}$, as well as the convergence of mixed approximation in the problems of classical and gradient elasticity theories. Note that inequality (38) is in fact equivalent to the stability condition of the mixed approximation formulated in [60, 64] for the general case of elliptic problems of mathematical physics, in which the solution is reduced to the determination of the saddle point of the alternating functional. Experience in solving problems of mixed approximation indicates that ignoring condition (38) can lead to ill-conditioned discrete problems in which the solutions have an unstable oscillating character.

Then, it is possible to define the scalar product and the norm in ${U}_{h}$as

$$ \begin{aligned} \left( {{\mathbf{u}}_{h},{\mathbf{v}}_{h} } \right)_{{U_{h} }} & = \left( {{\text{I}}_{h}{\text{B}}{\mathbf{u}}_{h} ,{\text{I}}_{h}{\text{B}}{\mathbf{v}}_{h} } \right)_{L}\\&\quad + \left( {{\text{NI}}_{h}{\text{B}}{\mathbf{u}}_{h} ,{\text{NI}}_{h}{\text{B}}{\mathbf{v}}_{h} } \right)_{H} , \forall {\mathbf{u}}_{h},{\mathbf{v}}_{h} \in U_{h} , \\ || {\mathbf{v}}_{h} ||_{{U_{h} }} &= \left( {{\mathbf{v}}_{h} ,{\mathbf{v}}_{h} } \right)_{{U_{h}}}^{1/2} , \forall {\mathbf{v}}_{h} \in U_{h} . \\ \end{aligned}$$

(39)

It will be shown that the solution of the variational Eq. (37) and, therefore, the solution of the equations of the mixed FEM (33) exists and is unique. To this end, consider a symmetric bilinear form ${a}_{h}\left(\cdot ,\cdot \right):{U}_{h}\times {U}_{h}\to {\mathbb{R}}$:

$${{a}_{h}\left({\mathbf{u}}_{h},{\mathbf{v}}_{h}\right)=\left({{\text{DI}}}_{h}{\text{B}}{\mathbf{u}}_{h},{{\text{I}}}_{h}{\text{B}}{\mathbf{v}}_{h}\right)}_{L}+{l}^{2}{\left({{\text{QNI}}}_{h}{\text{B}}{\mathbf{u}}_{h},{{\text{NI}}}_{h}{\text{B}}{\mathbf{v}}_{h}\right)}_{H}, \forall {\mathbf{u}}_{h},{\mathbf{v}}_{h}\in {U}_{h}.$$

(40)

Since the operator D is self-adjoined and positively defined, there are two positive constants ${m}_{0}$ and ${M}_{0}$, such that

$${m}_{0}{|| {{\text{I}}}_{h}{\text{B}}{\mathbf{v}}_{h}||}_{L}^{2}{\le \left({{\text{DI}}}_{h}{\text{B}}{\mathbf{v}}_{h},{{\text{I}}}_{h}{\text{B}}{\mathbf{v}}_{h}\right)}_{L}\le {{M}_{0}|| {{\text{I}}}_{h}{\text{B}}{\mathbf{v}}_{h}||}_{L}^{2}, \forall {\mathbf{v}}_{h}\in {U}_{h},$$

(41)

where

$$ m_{0} = 2 ess.{\text{~}}\mathop {inf}\limits_{{{{~}}\Omega {{~~~}}}} \mu ,M_{0} = ess.{{~}}\mathop {sup}\limits_{{{{~}}\Omega {{~~}}}} \left( {3\lambda + 2\mu } \right), $$

(42)

Moreover, following the definition of operator ${\text{Q}}$,we have

$$\begin{aligned}{{m}_{0}||{{\text{NI}}}_{h}{\text{B}}{\mathbf{v}}_{h}||}_{H}^{2}\le&\,{\left({{\text{QNI}}}_{h}{\text{B}}{\mathbf{v}}_{h},{{\text{NI}}}_{h}{\text{B}}{\mathbf{v}}_{h}\right)}_{H}\\\le&\,{{M}_{0}|| {{\text{NI}}}_{h}{\text{B}}{\mathbf{v}}_{h}||}_{H}^{2},\forall {\mathbf{v}}_{h}\in {U}_{h}.\end{aligned}$$

(43)

Accounting for Eqs. (40)–(43), we obtain

$$\begin{aligned}\left(1+{l}^{2}\right){m}_{0}{||{\mathbf{v}}_{h}||}_{{U}_{h}}^{2}\le&\,{a}_{h}\left({\mathbf{v}}_{h},{\mathbf{v}}_{h}\right)\\ \le&\,\left(1+{l}^{2}\right){M}_{0}{|| {\mathbf{v}}_{h}||}_{{U}_{h}}^{2},\forall {\mathbf{v}}_{h}\in {U}_{h}.\end{aligned}$$

(44)

Thus, the bilinear form defined by (40) possesses continuity and ellipticity properties in space ${U}_{h}\times {U}_{h}$, where the scalar product and norm are introduced by (39). Therefore, the solution of the equations of the mixed FEM (33) taking into account the assumption (36) exists and is unique.

3.2 Solution method

Application of FEM to the solution of the problems of gradient elasticity theory using mixed variation formulation leads to a system of linear algebraic equations of significantly higher order compared to the number of equations in displacements for problems of classical elasticity theory. In addition, the global matrix of the equation system for the mixed method is not positively defined and, in general, is the dense matrix. Such matrix properties practically exclude using direct methods for solving a system of algebraic equations (the Gauss method and its various modifications).

Here, the system of algebraic equations of the mixed method is converted to an equivalent form for nodal displacements. To this end, interpolation quadrature formulas are used, the weighting points of which coincide with the interpolation nodes of the triangular finite element. As a result, the solution of the system of the mixed method is significantly simplified, since there is no need to use computational algorithms for inverting the Gram matrix to find deformations and stresses from the first and second equations of the system (33). Thus, the procedure for determining nodal deformations and stresses is reduced to calculating those using recurrent formulas. In this case, the global unknowns in the discrete problem are the displacements at the mesh nodes since the node values of strains and stresses are calculated using the “averaging” procedure over the triangular elements adjacent to the mesh nodes and the values of strains, stresses and their derivatives at the centers of the triangles are considered as local unknowns.

Using the above procedure for the strains and stresses determination in mesh nodes leads to a positive definiteness of a symmetric global stiffness matrix. The matrix is sparse, but the number of its non-zero elements is significant. In this case, the application of direct methods for solving a system of equations of the mixed FEM remains problematic from the point of view of the computational implementation. In this regard, we have to use the iterative algorithm for solving equations of the mixed FEM.

To solve system (33), consider the iterative process in the form of the method of successive corrections. Let ${\mathbf{u}}_{h}^{0}\in {U}_{h}$ be an arbitrary initial approximation to the ${\mathbf{u}}_{h}$solution. Then, the process of successive linear approximations is constructed as follows:

$$ \begin{aligned} & \left({{\mathbf{E}}_{2h}^{k} ,\delta {\mathbf{T}}_{2h} } \right)_{L} = \left( {{\text{B}}{\mathbf{u}}_{h}^{k} ,\delta {\mathbf{T}}_{2h} } \right)_{L} , \forall \delta {\mathbf{T}}_{2h} \in X_{h} , \\ & \left( {{\mathbf{T}}_{2h}^{k} ,\delta {\mathbf{E}}_{2h} } \right)_{L} = \left( {{\text{D}}{\mathbf{E}}_{2h}^{k} ,\delta {\mathbf{E}}_{2h} } \right)_{L} , \forall \delta {\mathbf{E}}_{2h} \in X_{h} , \\ & \left( {{\mathbb{E}}_{3h}^{k} ,\delta {\mathbb{T}}_{3h} } \right)_{H} = \left( {{\text{N}}{\mathbf{E}}_{2h}^{k} ,\delta {\mathbb{T}}_{3h} } \right)_{H} , \forall \delta {\mathbb{T}}_{3h} \in Z_{h} , \\ & \left( {{\mathbb{T}}_{3h}^{k} ,\delta {\mathbb{E}}_{3h} } \right)_{H} = l^{2} \left( {{\text{N}}{\mathbf{T}}_{2h}^{k} ,\delta {\mathbb{E}}_{3h} } \right)_{H} , \forall \delta {\mathbb{E}}_{3h} \in Z_{h} , \\ & \left( {{\tilde{\text{D}}\text{B}}\Delta {\mathbf{u}}_{h}^{k + 1} ,{\text{B}}\delta {\mathbf{u}}_{h} } \right)_{L} = \left( {{\mathbf{T}}_{2h}^{k} ,\delta {\mathbf{E}}_{2h} } \right)_{L} + \left( {{\mathbb{T}}_{3h}^{k} ,\delta {\mathbb{E}}_{3h} } \right)_{H} \\ & \qquad\qquad\qquad\qquad\qquad\qquad -\left\langle {{\mathbf{f}},\delta {\mathbf{u}}_{h} } \right\rangle, \forall \delta {\mathbf{u}}_{h} \in U_{h} ; \\ & {\mathbf{u}}_{h}^{k + 1} = {\mathbf{u}}_{h}^{k} - \omega \Delta {\mathbf{u}}_{h}^{k + 1} , \\ \end{aligned}$$

(45)

where ${\Delta \mathbf{u}}_{h}^{k+1}$ is the correction term for displacement approximation ${\mathbf{u}}_{h}^{k}$.

Parameter $\omega >0$ and operator $\widetilde{D}$ are introduced to control the convergence of the iteration process (45). Their choice depends on the properties of the problem, the approximating functions, the finite element length $h$, and the structural (scale) parameter $l$, i.e., on their ratio $l/h$. It is assumed that $\widetilde{D}$ is a linear self-adjoined positively definite operator corresponding to the tensor of fictitious elastic moduli, which are determined by the relations:

$$\widetilde{\lambda }=\gamma \lambda , \widetilde{\mu }=\gamma \mu , \widetilde{E}=\gamma E, \gamma =1+C{\left( \frac{l}{h} \right)}^{2},$$

(46)

where $C$ is a constant independent of the parameter $h$.

Remark 3.3

If condition (38) is satisfied, iterative process (45), (46) converges to the solution of the mixed method Eqs. (33) at a geometric progression rate at any initial approximation of ${\mathbf{u}}_{{\text{h}}}^{0}\in {{\text{U}}}_{{\text{h}}}$.The proof is constructed by reducing the system of Eqs. (45) to one equation for displacements:

$$\begin{aligned} \left( {{ {\tilde{\text {D}}{\rm B}}}{\mathbf{u}}_{h}^{{k + 1}} ,{\text{B}}{\mathbf{v}}_{h} } \right)_{L} & = \left( {{ {\tilde{\text {D}}{\rm B}}}{\mathbf{u}}_{h}^{k} ,{\text{B}}{\mathbf{v}}_{h} } \right)_{L} \\ & \quad - \omega \left\{ {\left( {{\text{DI}}_{h} {\text{B}}{\mathbf{u}}_{h}^{k} ,{\text{I}}_{h} {\text{B}}{\mathbf{v}}_{h} } \right)_{L} + l^{2} \left( {{\text{NDI}}_{h} {\text{B}}{\mathbf{u}}_{h}^{k} ,{\text{NI}}_{h} {\text{B}}{\mathbf{v}}_{h} } \right)_{H} - \left\langle {{\mathbf{f}},{\mathbf{v}}_{h} } \right\rangle } \right\},\;~\forall {\mathbf{v}}_{h} \in U_{h} . \\ \end{aligned}$$

(47)

Then, according to the general results of the analysis of iterative two-layer schemes, we conclude that taking into account the inequality ${{\text{a}}}_{{\text{h}}}\left({\mathbf{v}}_{{\text{h}}},{\mathbf{v}}_{{\text{h}}}\right)>0$ for any ${\mathbf{v}}_{{\text{h}}}\in {{\text{U}}}_{{\text{h}}}$, the value of the parameter $\upomega$ >0 can be chosen in such a way as to ensure the convergence of the iterative process [57]. Convergence of the iterative scheme (45), (46) is considered more detailed in “Appendix”.

3.3 Basic relations of two-dimensional boundary value problems (plane strain)

The solution of the set of modeling two-dimensional boundary value problems of plate tension with central and edge cracks was realized under plane strain based on a mixed FEM scheme. It is assumed that a plate occupying a domain in the (x, y)-plane where the z-axis is normal to this plane. All tractions acts in the plane (x, y) and don’t depend on z. The 2D displacements are following:

$${u}_{x}(x,y)\ne 0,{u}_{y}(x,y)\ne 0,{u}_{z}\left(x,y\right)=0.$$

(48)

In the plane-strain state, the independent components of the stress tensors acts in the plane (x, y).

Constitutive relations for isotropic material have form:

$$ \begin{aligned} \sigma_{xx} & = \left( {\lambda + 2\mu } \right)u_{x,x} + \lambda u_{y,y} , \\ \sigma_{yy} & = \left( {\lambda + 2\mu } \right)u_{y,y} + \lambda u_{x,x} , \\ \sigma_{xy} & = \mu \left( {u_{x,y} + u_{y,x} } \right) , \\ \end{aligned} $$

(49)

Equilibrium equations:

$$ \begin{aligned} \left( {1 - l^{22} } \right)\left( {\sigma_{xx,x} + \sigma_{xy,y} } \right) & = 0, \\ \left( {1 - l^{22} } \right)\left( {\sigma_{yx,x} + \sigma_{yy,y} } \right) & = 0, \\ \end{aligned} $$

(50)

where

$${\nabla }^{2}\left(\dots \right)={\partial }_{x}^{2}\left(\dots \right)+{\partial }_{y}^{2}\left(\dots \right)$$

(51)

and $\lambda$ and $\mu$ are Lamé parameters and $l$ has the length dimension.

3.4 Numerical analysis

To approximate the displacement, piecewise linear functions on triangular elements were used, and for interpolation of the strains and the stresses, a linear combination of piecewise polynomial approximations and bulb functions was employed [38]. With such a choice of approximating functions, the stability and convergence of an approximate solution to the boundary value problem are provided.

To calculate the scalar product ${\left(\cdot ,\cdot \right)}_{L}$ in (45), interpolation quadrature formulas were utilized, in which weighting points coincide with the interpolation nodes of the triangular finite element. As a result, the solution of the system of equations of the mixed method is significantly simplified since there is no need to apply computational algorithms for inverting the Gram matrix to find strains and stresses from the first and second equations of the system (45), respectively, because the matrix becomes diagonal. In addition, linear interpolation of strains and stresses within the triangular element was used to approximate the strain and stress gradient. Then, the distributions of the derivatives of the strains and stresses are described using piecewise constant functions, and therefore, the scalar product ${\left(\cdot ,\cdot \right)}_{H}$ in (45) is calculated accurately. Thus, the determination of the derivatives from strains and stresses is reduced to their calculation by the recurrent formulas at each iteration. According to the accepted approximation, the global unknowns in the discrete problem are the displacements in the mesh nodes since the node values of the strains and the stresses are calculated using the "averaging" procedure over the triangular elements adjacent to the mesh nodes. The values of strains, stresses, and their derivatives in the centers of triangles are considered as local unknowns.

Remark 3.4

It should be noted that the iterative process (45) is similar in a certain sense to the computational procedure for determining continuous stress fields proposed by Loubignac et al. [65]. However, this analogy is completely formal since the construction of an approximate solution of the equations of the gradient elasticity theory is based on the concept of mixed approximation of the boundary value problem.

Values of the iterative parameter in (45) and (46) were determined approximately based on a series of preliminary calculations. The following values were used for the approximating functions:

$$\omega \approx 1\div \mathrm{1,25}, \gamma =1+C\frac{ {l}^{2}}{\Delta }, C\sim 1\div 2,$$

(52)

where $\Delta$ is the area of a triangle element.

Calculations were realized, considering the symmetry of the problem using a uniform and nonuniform triangular mesh (Fig. 1). For a non-uniform triangular mesh, uniform partitioning was used in the vicinity of the crack tip. Symmetry conditions were paid due regard by taking normal displacements and tangent stresses equal to zero on the axes of symmetry.

Fig. 1

Uniform (a) and nonuniform (b) triangular meshes

For verification of the convergence of the FEM solution, calculations were provided for a uniform rougher mesh with element length h = 2 × 10^–3 m, i.e., the crack length was segmented into 20 elements, and for finer meshes with h = 10^–3 m, i.e., the crack length was segmented into 40 sections (Fig. 2) and with h = 5 × 10^–4 m, i.e., the crack length was cut into 80 pieces and nonuniformly partitioned into a triangular mesh in the vicinity of the crack tip: h = 10^–4 m (Fig. 3).

Fig. 2

Partition into a triangular mesh in the vicinity of the crack: $h={10}^{-3} \mathrm{m}$. The crack length is segmented into 40 sections

Fig. 3

Nonhomogeneous partition into a triangular mesh in the vicinity of the crack tip: $h={10}^{-4} {\text{m}}$

The tension of a square plate with a central crack. The square plate with a crack is shown in Fig. 4. The half-length of the plate edges is W = 0.2 m, and the crack length is a = 0.04 m. The plate is under uniaxial tensile stress q = 100 MPa. Due to symmetry, the problem is solved for a quarter of the plate under the assumption of the plane strain state. It is assumed that the plate is made of strain-gradient, isotropic material with the following properties E = 2.10⁵ MPa and ν = 0.3 where E is the Young modulus and ν is Poisson’s ratio. The normalized scale parameter l/a values are taken: 0.02, 0.05, 0.1, 0.2, 0.4, and 0.6.

Fig. 4

The tension of a square plate with a central crack

The boundary conditions

Static boundary conditions for stresses are natural; they are taken into account in a discrete problem using the Lagrange variational equation written in form of the last Eq. (33). Kinematic boundary conditions for displacements are determined by the symmetry conditions of the problem:

$$ \begin{aligned} u_{y} \left( {x,0} \right) & = 0,u_{x,y} \left( {x,0} \right) = 0,u_{y,x} \left( {x,0} \right) = 0, \\ u_{x} \left( {0,y} \right) & = 0,u_{y,x} \left( {0,y} \right) = 0,u_{x,y} \left( {0,y} \right) = 0. \\ \end{aligned} $$

(53)

The problem of a plate with a central crack has been widely studied in the literature (e.g., [6, 14, 31, 37, 49, 50, among others]). This problem is primarily considered to verify the proposed mixed FEM to solve the boundary value problems within the strain gradient elasticity theory.

Figure 5 presents a comparison of the results for the crack-opening displacement along the crack edge using the uniform partitions of the crack length: dashed lines are segmented into 20 elements; solid lines into 40 sections; dotted lines into 80 pieces and different values of the normalized scale parameter l/a = 0.1; 0.2; 0.3; 0.4; 0.6. Good convergence of the approximate FEM solutions is observed by refining the mesh even for high (unrealistic) values of the scale parameter $l/a$. This can be considered as a verification of the mixed FEM proposed for solving such a problem within the strain gradient elasticity theory.

Fig. 5

Variation in crack-opening displacement for uniform partitions of crack length (dashed lines correspond to partitions into 20 sections; solid lines correspond to partitions into 40 sections; dotted lines correspond to partitions into 80 sections) along the crack edge for different values of the normalized parameter $l/a$

Numerical results for the normalized crack-opening displacements $u/u_{cl}$ ($u_{cl}$ is the solution within classical linear elasticity) as a function of normalized distance from the crack tip $x/a$ for the uniform partition of the crack length into 80 sections and nonuniform segmentation of the crack edge with element length in the vicinity of the crack tip h = ${10}^{-3} {\text{m}}$ are shown in Figs. 6 and 7, respectively. The dashed line corresponds to the solution within classical linear elasticity (scale parameter $l/a=0$), and solid lines—to the solution within strain gradient elasticity with scale parameter l/a = 0.1; 0.2; 0.3; 0.4; 0.6. The diagrams of both figures are visually identical, which confirms the good convergence of the approximate solution obtained by the proposed mixed FEM. Moreover, curves corresponding to the solutions within classical linear elasticity and strain gradient elasticity with scale parameter $l/a=0.2;$ $0.6$ are in good agreement with the results obtained by FEM in [37] and with the results of the boundary element method presented in [31] (Fig. 9, parameter l/a is equivalent to $2g$ in [31], which means that $g$ = 0.1; 0.3). This is an additional verification of mixed FEM calculations presented here, especially because the results of [37] were determined using the Tri18 element belonging to a class of C¹ functions.

Fig. 6

Normalized crack-opening displacement as a function of normalized distance from the crack tip for uniform partitioning of the crack length into 80 elements and different values of the normalized parameter $l/a$

Fig. 7

Normalized crack-opening displacement as a function of normalized distance from the crack tip for nonuniform partition of crack length and different normalized parameter values $l/a$. Element length in the vicinity of the crack tip h = 10^–3 m

Remark 3.5

It is expected that the length scale parameters are dependent on the size of the micro-inhomogeneity and the distances between them or micro-structure of material. At any case, the length scale parameters are significantly less than the macro scale. Nevertheless, we present results for unrealistic values of the scale parameter l/a first of all to demonstrate a good convergence of the approximate FEM solutions and secondly to compare with results available in the literature. This can be considered as a validation of the mixed FEM proposed for solving such a problem within the strain gradient elasticity theory.

Figures 5, 6, 7 indicate that crack stiffness increases with the scale parameter, which is in agreement with other results available in the literature (e.g., [14, 31, 37, 50, among others]). As expected, the behavior of the crack profile is different for classical and strain-gradient solutions. The crack faces close more smoothly in the case of strain-gradient elasticity compared to the classical elasticity results. Cusp-like closure of the crack in strain-gradient solutions occurs because the asymptotic solution for crack-opening displacement near the crack tip ($r\to 0, r$ is the radial distance from the crack tip) is governed by dominant term $r^{3/2}$, which has higher order in comparison with classical elasticity ($r^{1/2}$). The cusp-like closure effect has been indicated in experiments by [66] and discussed in many papers (e.g., [14, 31, 33, 37, 50, 67, among others]).

For comparison, normalized crack-opening displacements $\mu u/(aq)$ calculated on the basis of mixed FEM (for the same parameters as in [14] and on the basis of hypersingular integral equations [14] are depicted in Fig. 8. The solid line corresponds to mixed gradient FEM solution for scale parameter $a/l=5$; small circles indicate results [14]. The results obtained by mixed FEM and those determined on the basis of hypersingular integral equations in Gourgiotis and Georgiadis [14] are visible identical, as shown in Fig. 8, in spite of rather fundamental differences between the methods involved and of unrealistic high length scale parameter ($a/l=5$) what significantly complicates mixed FEM calculations (convergence).

Fig. 8

Normalized crack-opening displacement $\mu u/(aq)$ as a function of normalized distance from the crack tip $x/a$ for scale parameter $a/l=5$. Solid line corresponds to mixed gradient FEM solution, small circles indicate results Gourgiotis and Georgiadis (2009)

Normalized normal strains $\mu {\varepsilon }_{yy}/q$ as a function of normalized distance from the crack tip $x/a$ are shown in Fig. 9. The solid line corresponds to results obtained by mixed gradient FEM solution for scale parameter $a/l=50$, small circles indicate results [14]. It is observed a good agreement between the both solutions. It is also seen that the normal strains take the finite values at the crack-tip $x/a=1$, while the corresponding strains in classical elasticity (dashed line) exhibit a square root singularity. In Fig. 10 the distribution of normalized shear strain $\mu {\varepsilon }_{yx}/q$ as a function of normalized distance from the crack tip $x/a$ for scale parameter $a/l=20$ is displayed. Contrary to the classical elasticity case, the shear strains are not zero at the crack-faces.

Fig. 9

Normalized normal strain $\mu {\varepsilon }_{yy}/q$ as a function of normalized distance from the crack tip $x/a$. Solid line corresponds to mixed gradient FEM solution for scale parameter $a/l=50$, small circles indicate results Gourgiotis and Georgiadis [14], dashed line—classical elasticity solution

Fig. 10

Normalized shear strain $\mu {\varepsilon }_{yx}/q$ as a function of normalized distance from the crack tip $x/a$ for scale parameter $a/l=20$

The tension of a rectangular plate with an edge crack. The problem for a rectangular plate with a symmetrically placed edge crack of length a = 0.4 m is shown in Fig. 11. The half-length of the vertical edges of the plate and the width are equal to W = 0.2 m. The plate is under uniaxial tensile stress q = 100 MPa. The material properties and the scale parameters are the same as in the above problem. Due to symmetry, the problem is solved for half of the plate under the assumption of the plane strain state.

Fig. 11

The tension of a rectangular plate with an edge crack

The boundary conditions

Kinematic boundary conditions for displacements are:

$$ \begin{aligned} u_{y} \left( {x,0} \right) & = 0,u_{x,y} \left( {x,0} \right) = 0,u_{y,x} \left( {x,0} \right) = 0, \\ u_{x} \left( {W,0} \right) & = 0. \\ \end{aligned} $$

(54)

Figures 12 and 13 demonstrate the crack-opening displacements as a function of distance from the plate edge with the crack for the uniform partitioning of the crack face into 40 sections and nonuniform segmentation of the crack edge with element length in the vicinity of the crack tip h = 10^–3 m, respectively. The dashed line indicates the solution within classical linear elasticity (scale parameter $l=0$) and solid lines—the solution within strain gradient elasticity with normalized scale parameter l/a = 0.1; 0.2; 0.3; 0.4; 0.6. Diagrams of both figures are visually identical, confirming good convergence of the approximate solutions obtained by the considered mixed FEM.

Fig. 12

Crack-opening displacements as a function of distance from the plate edge with the crack for uniform partitioning of the crack length into 40 sections and different values of the normalized parameter $l/a$

Fig. 13

Crack-opening displacements as a function of distance from the plate edge with crack for nonuniform partition of the crack length and different normalized parameter values $l/a$. Element length in the vicinity of the crack tip h = 10^–3 m

For an additional verification of the obtained solutions, Fig. 14 presents a comparison of the crack-opening displacement along the crack edge for the uniform partitions of the crack length: dashed lines are segmented into 20 sections; solid lines are segmented into 40 sections; dotted lines are segmented into 80 pieces and for different values of the normalized scale parameter l/a = 0.1; 0.2; 0.3; 0.4; 0.6. Good convergence of the mixed FEM solutions is observed by refining the mesh for all values of parameter $l/a$.

Fig. 14

Variation in crack-opening displacement for uniform partitions of crack length (dashed lines correspond to partitions into 20 sections; solid lines correspond to partitions into 40 sections; dotted lines correspond to partitions into 80 sections) along the crack edge for different values of the normalized parameter $l/a$

The same tendencies as for the problem with the central crack are observed for the plate with an edge crack (Figs. 12, 13, 14): increasing the crack stiffness with the scale parameter value, a cusp-like closure effect. The difference is just in the value of the crack-opening displacements. As expected, the opening displacements of the edge crack are higher than those of the central crack for all values of scale parameters.

The tension of a square plate with two symmetrically placed edge cracks. The problem for a square plate with two symmetrically placed edge cracks of length a = 0.04 m is shown in Fig. 15. The half-length of the plate edges is equal to W = 0.2 m. The plate is under uniaxial tensile stress q = 100 MPa. The material properties and the scale parameters are the same as in the above two problems. Due to symmetry, the problem is solved for the quarter of the plate under the assumption of the plane strain state.

Fig. 15

The tension of a square plate with two symmetrically placed edge cracks

The boundary conditions

Kinematic boundary conditions for displacements are:

$$ \begin{aligned} u_{y} \left( {x,0} \right) & = 0,u_{x,y} \left( {x,0} \right) = 0,u_{y,x} \left( {x,0} \right) = 0, \\ u_{x} \left( {0,y} \right) & = 0,u_{y,x} \left( {0,y} \right) = 0,u_{x,y} \left( {0,y} \right) = 0. \\ \end{aligned} $$

(55)

Figures 6, 7, 8, 9 and 16 indicate that increasing the crack stiffness with scale parameter value and the cusp-like closure effect are common to all three problems of the plate with the central, edge, and two symmetrical edge cracks. The difference is just in the value of the crack-opening displacements. For the plate with two symmetrical edge cracks, the opening displacements are higher in comparison with those of the central crack and lower in comparison with those of the single-edge crack for all values of scale parameters.

Fig. 16

Crack-opening displacements of the symmetric edge cracks as a function of distance from the plate edge with the crack for uniform partitioning of the crack length into 40 elements and different values of the normalized parameter $l/a$

Figure 17 demonstrates good convergence of the FEM solutions by refining the mesh for all values of scale parameters $l/a$ and is presented for verification of the FEM solution. A comparison of the crack-opening displacement of the symmetric edge cracks along the crack edge for the uniform partitions of the crack length is shown: dashed lines—segmented into 20 sections, solid lines—into 40 sections, and dotted lines—the nonuniform partitions with element length in the vicinity of the crack tip h = ${10}^{-3} {\text{m}}$ and for different values of the normalized scale parameter l/a = 0.1; 0.2; 0.3; 0.4; 0.6.

Fig. 17

Variation in the crack-opening displacement of the symmetric edge cracks for uniform and nonuniform partitions of crack length (dashed lines represent the uniform partition into 20 elements; solid lines represent the uniform partitioning into 40 sections; dotted lines represent the nonuniform partition with element length in the vicinity of the crack tip h = ${10}^{-3} {\text{m}}$) along the crack edge for different values of the normalized parameter $l/a$

4 Conclusions

An alternative approach has been applied to solve the boundary value problem within the gradient elasticity theory in that a mixed variation formulation of FEM is proposed. This approach is based on the concept of the Galerkin method, and unlike the classical formulation of FEM in the form of a displacement method, displacements, strains, stresses, and their gradients are direct arguments of the variational equations of the finite-dimensional problem, and are not determined by further transformation (differentiation) of the problem solution in displacements. To construct finite-dimensional subspaces based on such variational formulation of the problem, separate (formally independent) approximations of displacements, deformations, stresses, and their gradients are used by choosing the different sets of piecewise polynomial basis functions, interconnected by the stability condition of the mixed approximation. Mixed FEM significantly simplifies the pre-requirement for approximating functions to belong to class C¹. Thus, using the simplest triangular finite elements with a linear approximation of displacements is possible under the condition that uniform or near-uniform triangulation is realized. Uniform triangulation in the vicinity of the stress concentrators provides super-convergence in calculating strains and stresses in mesh nodes. Global unknowns in a discrete problem are nodal displacements, while the strains, stresses, and their gradients are considered as local unknowns.

The conditions of existence, uniqueness, and continuous dependence of the solution on the initial data of the problem have been formulated for discrete equations of mixed FEM. The system of equations of mixed FEM has been solved by a modified iteration procedure with a preconditioning matrix. The global stiffness matrix for classical elasticity problems is considered as a preconditioning matrix with fictitious elastic moduli. Then, there is no need to form a global stiffness matrix for the problem of strain gradient elasticity since it is enough to calculate only the residual vector in the current approximation. This reduces the memory resource requirements of the computer.

A set of modeling plane crack problems has been solved, namely, the tension of the plate with central, single-edge, and two (symmetrically placed) edge cracks. The first problem has been considered in the context of similar results available in the literature. It indicates good agreement with those obtained by FEM in Papanicolopulos and Zervos [37], by the boundary element method in [31] and and by an analytical/numerical technique based on hypersingular integral equations in Gourgiotis and Georgiadis [14]. This verification of mixed FEM calculations is especially important because the results of [37] were obtained using the Tri18 element belonging to a class of C¹ functions. Furthermore, good agreement between the results of mixed FEM and those of Gourgiotis and Georgiadis [14], despite of rather fundamental differences between the methods involved can be interpreted as an additional validation of the proposed approach. To the best author’s knowledge, edge crack problems are considered within strain gradient elasticity theory for the first time. Good convergence has been observed by refining the mesh (finite element length) for all scale parameters, even for unrealistically high ones. Moreover, all three problems demonstrate specific qualitative features characterizing stain gradient solutions, increasing crack stiffness with length scale parameter and cusp-like closure effect, and agree with similar theoretical developments in the literature. Therefore, it should be noted that the mixed FEM based on the accepted approximation for displacements, strains, stresses, and their gradients on triangular elements can be successfully applied to solving two-dimensional boundary value problems within the strain gradient elasticity theory.

Appendix: Convergence of the iterative scheme

Let us write down the variational Eq. (37) in the operator's form:

$${\text{A}}{\mathbf{u}}_{h}= {\mathbf{f}}_{h},{\mathbf{f}}_{h}\in {U}_{h}^{*},$$

(56)

where ${U}_{h}^{*}$—space conjugated to space ${U}_{h}$; ${\text{A}}$—a linear symmetric coercive operator acting from ${U}_{h}$ to ${U}_{h}^{*}$ and defined by equation:

$$ {\text{A}}: {\mathbf{u}}_{h} \in U_{h} \to \left( {{\text{DI}}_{h} {\text{B}}{\mathbf{u}}_{h} ,{\text{I}}_{h} {\text{B}}{\mathbf{v}}_{h} } \right)_{L} + l^{2} \left( {{\text{QNI}}_{h} {\text{B}}{\mathbf{u}}_{h} ,{\text{NI}}_{h} {\text{B}}{\mathbf{v}}_{h} } \right)_{H} , \forall {\mathbf{v}}_{h} \in U_{h} . $$

(57)

To solve Eq. (56) let’s consider iterative process:

$${\mathbf{u}}_{h}^{k+1}={\mathbf{u}}_{h}^{k}-\omega {{\text{K}}}^{-1}\left({\text{A}}{\mathbf{u}}_{h}^{k}-{\mathbf{f}}_{h}\right),$$

(58)

where $\omega$ is the parameter introduced to control process convergence; ${\text{K}}$ is a linear symmetric coercive operator acting from ${U}_{h}$ to ${U}_{h}^{*}$, corresponding to the equations of the classical elasticity theory in displacement:

$$K:{ \mathbf{u}}_{h}\in {U}_{h}\to {\left({\text{DB}}{\mathbf{u}}_{h},{\text{B}}{\mathbf{v}}_{h}\right)}_{L}, \forall { \mathbf{v}}_{h}\in {U}_{h}.$$

(59)

According to the general results of the convergence analysis of iterative methods, process of successive approximations (58) converges to the solution of Eq. (56) if the parameter $\omega$ satisfies the condition [61]:

$$0<\omega <2/\gamma ,$$

(60)

where $\gamma$ is the positive constant, defined from inequality:

$$\left({\text{A}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right)\le \gamma \left({\text{K}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right), \forall { \mathbf{v}}_{h}\in {U}_{h}.$$

(61)

Accounting for the properties of the orthogonal projector ${{\text{I}}}_{h}$, the Cauchy-Bunyakovsky inequality and the “inverse inequalities” we find:

$$ \begin{aligned} \left( {{\text{A}}{\mathbf{v}}_{h} ,{\mathbf{v}}_{h} } \right) & \le \left( {{\text{DB}}{\mathbf{v}}_{h} ,{\text{B}}{\mathbf{v}}_{h} } \right)_{L} + C\left( {l/h } \right)^{2} \left( {{\text{DI}}_{h} {\text{B}}{\mathbf{v}}_{h} ,{\text{I}}_{h} {\text{B}}{\mathbf{v}}_{h} } \right)_{L} \\ & \le \gamma \left( {{\text{DB}}{\mathbf{v}}_{h} ,{\text{B}}{\mathbf{v}}_{h} } \right)_{L} = \gamma \left( {{\text{K}}{\mathbf{v}}_{h} ,{\mathbf{v}}_{h} } \right), \\ \end{aligned} $$

(62)

with

$$\gamma =1+C{\left(l/h \right)}^{2}.$$

(63)

Here $C$ is a positive constant; h is the step of the finite element mesh.

Given Eq. (63), it follows that for dense partitions we have $l/h\gg 1$ and, therefore, $\gamma \gg 1$, and as a result $\omega \ll 1$. In other words, by the mesh refinement, the parameter $\omega$ decreases, and its value can be significantly less than unity, i.e., the choice of the parameter $\omega$ depends significantly on the mesh size h. The use of the parameter $\omega \ll 1$ using the parameter $\omega$ may lead to instability and misconvergence of the computational process.

Let us modify the iterative process (58) in such a way as to provide stable value of the parameter $\omega$, which slightly depends on the mesh size $h$. To this end, the operator ${\text{K}}$ is transformed using fictitious elastic material moduli $\widetilde{\lambda },\widetilde{\mu }$, which are defined within the finite element by the following relations:

$$\widetilde{\lambda }=\gamma \lambda , \widetilde{\mu }=\gamma \mu .$$

(64)

The iterative process (47) can then be written in equivalent operator form

$${\mathbf{u}}_{h}^{k+1}={\mathbf{u}}_{h}^{k}-\omega {\widetilde{{\text{K}}}}^{-1}\left({\text{A}}{\mathbf{u}}_{h}^{k}-{\mathbf{f}}_{h}\right),$$

(65)

where $\widetilde{{\text{K}}}$ is a linear operator defined analogically as the operator ${\text{K}}$, but for fictitious elastic material moduli $\widetilde{\lambda },\widetilde{\mu }$.

To estimate the upper values of the parameter $\omega$ inequality (61) is used:

$$\left({\text{A}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right)\le \left(\widetilde{{\text{K}}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right), \forall { \mathbf{v}}_{h}\in {U}_{h}.$$

(66)

On the basis on inequality (66), it is possible to obtain that the upper estimation of the permissible values of the parameter $\omega$ does not depend on the mesh step h, and the iterative process (65) convergences under the following condition:

$$\omega \in \left(0, 2\right).$$

(67)

The optimal value of the parameter ${\omega }_{opt}$ providing the maximal convergence rate of the iteration process (65) is defined by the formula

$${\omega }_{opt}=2/\left(1+\alpha \right),$$

(68)

where $\alpha$ is a positive constant, defined by inequality:

$$\left({\text{A}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right)\ge \alpha \left(\widetilde{{\text{K}}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right), \forall { \mathbf{v}}_{h}\in {U}_{h}.$$

(69)

Utilizing stability condition Eq. (38) have:

$$\left({\text{A}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right)\ge {d}_{0}^{2}\left({\text{K}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right), \forall { \mathbf{v}}_{h}\in {U}_{h},$$

(70)

and taking into account the inequality

$$\left(\widetilde{{\text{K}}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right)\le \gamma \left({\text{K}}{\mathbf{v}}_{h},{\mathbf{v}}_{h}\right), \forall { \mathbf{v}}_{h}\in {U}_{h},$$

(71)

constant $\alpha$ can be determined:

$$\alpha ={d}_{0}^{2}/\gamma .$$

(72)

From Eqs. (63), (68), (72), it follows that for fine partitions the value of the parameter ${\omega }_{opt}$ non-essentially depends on the mesh size h. The positive constant C in Eq. (63) is choosing to fulfill condition ${\omega }_{opt}\approx 1$. Introducing the operator $\widetilde{{\text{K}}}$ and choice of the parameter $\omega$ with accounting for the above considerations ensures the stability and convergence of the iterative process (65).