## Abstract

In this paper, both singular and hypersingular fundamental solutions of plane Cosserat elasticity are derived and given in a ready-to-use form. The hypersingular fundamental solutions allow to formulate the analogue of Somigliana stress identity, which can be used to obtain the stress and couple-stress fields inside the domain from the boundary values of the displacements, microrotation and stress and couple-stress tractions. Using these newly derived fundamental solutions, the boundary integral equations of both types are formulated and solved by the boundary element method. Simultaneous use of both types of equations (approach known as the dual boundary element method (BEM)) allows problems where parts of the boundary are overlapping, such as crack problems, to be treated and to do this for general geometry and loading conditions. The high accuracy of the boundary element method for both types of equations is demonstrated for a number of benchmark problems, including a Griffith crack problem and a plate with an edge crack. The detailed comparison of the BEM results and the analytical solution for a Griffith crack and an edge crack is given, particularly in terms of stress and couple-stress intensity factors, as well as the crack opening displacements and microrotations on the crack faces and the angular distributions of stresses and couple-stresses around the crack tip.

## 1. Introduction

This paper presents both a complete derivation of the fundamental solutions arising in plane Cosserat elasticity and a demonstration of their use within a boundary element framework based on the simultaneous use of two types of boundary integral equations (BIEs) on the crack faces.

The Cosserat (also known as micropolar) theory of elasticity was first introduced by the Cosserat brothers [1] and further developed by Eringen [2], Nowacki [3] and others as a generalization of the classical elasticity, which takes into account the effects of a material's microstructure by enriching material infinitisimal elements with additional rotational degrees of freedom and introducing material internal length-scale parameters directly into the constitutive equations. This theory is known to represent well a number of natural and engineered materials, e.g. fibre-reinforced composites, metal foams, concrete, synthetic polymers and human bones.

A number of analytical and numerical methods, which have been successfully used to treat boundary value problems in classical elasticity, have been developed for micropolar elasticity as well. While the finite-element method (FEM) remains the most popular tool of numerical analysis, the boundary element method (BEM) is evolving as an efficient alternative, especially in modelling problems with discontinuities.

The mathematical foundation of BEM is the BIE method, which was introduced in [4,5] to study the solvability of boundary value problems in plane micropolar elasticity. It has been shown that the solutions of these boundary value problems can be found in terms of single layer and double layer potentials, i.e. in the form of an integral of the product of the unknown densities and the kernel functions, known as the fundamental solutions. This representation allows one to reduce a boundary value problem to the systems of weakly singular, singular and hyper-singular BIEs, which can be subsequently solved by the BEM.

The BEM for singular integral equations, which are also known as displacement/microrotation boundary integral equations (DBIEs), was developed in [6–8]. To the authors' knowledge, no BEM solutions for the traction boundary integral equations (TBIEs) of micropolar elasticity have been published yet. This might be attributed to the fact that the hyper-singular fundamental solutions are available in the literature only in an implicit form, which requires additional derivations before it can be implemented into a computer code. However, application of the DBIE-based BEM to crack problems results in a singular system matrix, due to the coincident crack surfaces, and therefore it is limited to problems with symmetry [9], which avoid use of one of the coincident boundaries. The most common approach to overcome this difficulty, which has been successfully used in classical elasticity, is known as the dual BEM and it consists in simultaneous use of the both types of the BIEs on the crack surfaces.

For this purpose, in this work, the BEM is developed for the BIEs of both types and its accuracy is demonstrated in various numerical examples, including problems with cracks. Only a limited amount of research in Cosserat fracture is available in the literature and it is limited to the studies of some specific cases (e.g. [10–13]). This work contains the first BEM approach to a general case of a crack problem in Cosserat elasticity. For the sake of comparison with the analytical solution, a classical problem of a Griffith crack is chosen as a numerical example. The obtained data are shown to be in excellent agreement with the exact solutions in terms of stress and couple-stress intensity factors and the crack opening displacements and microrotations along the crack faces.

The BEM results are also compared with the analytical asymptotical solutions and the FEM data obtained in [13] for an edge crack problem in Mode I and Mode II in terms of angular distributions of the stresses and couple-stresses around the crack tip.

The paper is organized as follows. Section 2 contains mathematical foundations of Cosserat elasticity; in §3, the BIEs are formulated; in §4, the BEM is outlined. Numerical examples are presented in §5 and the main results are summarized in §6.

## 2. Mathematical foundations of plane micropolar elasticity

In what follows, a linear homogeneous micropolar elastic solid is considered. The solid occupies an open domain *S*. The plane strain state is described by two in-plane displacements *u*_{1},*u*_{2} and one out-of-plane microrotation *ϕ*_{3}. In addition to stresses *σ*_{11},*σ*_{12},*σ*_{21},*σ*_{22} (note, that *σ*_{12}≠*σ*_{21}), two couple-stresses *m*_{13},*m*_{23} are introduced (figure 1), and the following constitutive equations hold [4].
*δ*_{ij} is the Kronecker delta, *ε*_{ijk} is the alternating symbol and λ,*μ*,*κ*,*γ* are material parameters.

In absence of body forces and body couples, the equilibrium equations [4] are given as
*u*_{3}=*ϕ*_{3} and introduce a vector of generalized displacements ** u**=(

*u*

_{1},

*u*

_{2},

*u*

_{3})

^{T}and generalized tractions

**=(**

*t**t*

_{1},

*t*

_{2},

*t*

_{3})

^{T}, where boundary tractions

*t*

_{1},

*t*

_{2}and a couple traction

*t*

_{3}are defined as

**=(**

*n**n*

_{1},

*n*

_{2})

^{T}is a unit outward normal to ∂

*S*.

The equations of equilibrium (2.2) then can be rewritten in the form
*L*(∂_{x})=*L*(*ξ*_{α}) is given by [5]
*ξ*_{α}=∂/∂*x*_{α} and

Together with *L*(*ξ*_{α}), the boundary stress operator *T*(∂_{x})=*T*(*ξ*_{α}) is considered [5], which is defined by the following equation:
*T*(∂_{x}) is defined according to (2.1) and (2.3) in such a way that
*S*_{u}, and *S*_{t}.

Together with λ,*μ*,*κ*,*γ*, the following (engineering) micropolar material constants are introduced in [16]:
*N*=0 or *l*=0, the micropolar theory reduces to classical elasticity. Case *N*=1 corresponds to well-known so-called ‘couple-stress theory’, which has been studied independently, for example, in [17–20]. In this theory, the couple-stresses are taken into consideration; however, the microrotations are constrained, i.e. defined analogously to macrorotations in classical elasticity.

Two typical approaches to determine Cosserat material constants are experimental methods [16,21,22] and analytical derivation, based on various homogenization schemes for materials with periodic microstructure, which have been proposed, for example, in [23–26].

## 3. Boundary integral equations

In what follows, a source point is denoted by **x**=(*x*_{1},*x*_{2}) and a field point by **y**=(*y*_{1},*y*_{2}).

The matrix of fundamental solutions *D*=*D*(**x**,**y**) of system (2.4) is derived in [5], in the implicit form, using the method of associated matrices, described in [14], in the form
*L**(∂_{x}) is the adjoint of *L*(∂_{x}) (matrix consisting of cofactors of *L*(∂_{x})) and *t*(**x**,**y**) is given as
*a* and *k* are defined by
*K*_{0} is the modified Bessel function of order zero. The full expression of matrix *D* in polar coordinates
*D*, the matrix of singular solutions *P*=*P*(**x**,**y**;*n*_{y}) is introduced by
_{y} implies that in equation (2.6) *ξ*_{α}=∂/∂*y*_{α} and normal *n*_{y} is applied at the point **y**. The full expression of matrix *P* is given, for example, in [27]. It can be verified by direct differentiation that the columns of *D*(**x**,**y**) and *P*(**x**,**y**;*n*_{y}) satisfy (2.4) at all

In order to formulate TBIEs, two more matrices are required. Namely, matrix *H*=*H*(**x**,**y**;*n*_{x}) and *S*=*S*(**x**,**y**;*n*_{x},*n*_{y}) which are obtained by applying stress operator *T*(∂_{x}) to matrices *D*(**x**,**y**) and *P*(**x**,**y**;*n*_{y}):
*D*, *P*, *H* and *S* is presented in [28] with the final expressions in a ready-to-use form in both, symbolic and C/C++ formats.

In order to investigate the behaviour of matrices *D*, *P*, *H* and *S* in the vicinity of **x**=**y**, they are expanded in Taylor series in polar coordinates (3.4). Straightforward derivations show that as *D* are [4]
*P* and *H* with the highest order of singularity being *ρ*^{−1} are listed below:
*μ*^{′}=(1−2*ν*)/2(1−*ν*)−*N*^{2} [4]. Components *P*_{α3}, *P*_{3α}, *H*_{α3} and *H*_{3α} are weakly singular, i.e. of order

For matrix *S* as *μ*′′=*E*(1/(1−*ν*^{2})+2*N*^{2}), *μ*′′′=*EN*^{2}/(1+*ν*) and components *S*_{12},*S*_{21} are of order

The derivation of the BIEs in plane Cosserat elasticity is based on the following analogue of Somigliana identity [5]:
**x** tend to the boundary ∂*S*. Removing a small vicinity of a singular point **x** and using the Taylor series expansions of the kernel functions, a straightforward derivation yields
*P*_{12},*P*_{21} are singular and understood as Cauchy principal values. Factor *u*_{1}, *u*_{2} and microrotation *u*_{3}.

In order to derive the TBIE, the first operator *T*(∂_{x}) is applied to equation (3.12) (with arbitrary direction ** n**(

**x**)) to obtain the following analogue of Somigliana stress identity:

**(**

*n***x**) to be the normal to ∂

*S*and performing limiting process analogously to the procedure described above for DBIE, the following TBIE is obtained:

*H*

_{12},

*H*

_{21},

*S*

_{13},

*S*

_{23},

*S*

_{31},

*S*

_{32}are singular and sign

*S*

_{11},

*S*

_{22},

*S*

_{33}are hyper-singular and understood as Hadamard finite part integrals.

## 4. Boundary element method formulation

In this work, the standard BEM discretization with the quadratic Lagrange basis is applied to equations (3.13) and (3.15), which makes use of the following set of shape functions:
*S* is discretized with *N* elements ∂*S*_{n} and the shape functions (4.1) with λ_{0}=1, while the solutions *u*_{i}(**y**) and *t*_{i}(**y**) are approximated with the discontinuous basis with *N* linear algebraic equations for the unknown nodal values of generalized displacements and tractions.

## 5. Numerical results

In this section, four numerical examples are shown. In the first example, the problem of a bending plate is considered, and it is shown that the exact solution, given by the polynomials of second degree, can be reproduced with a good accuracy by BEM on a coarse mesh.

In the second example, the problem of the stress concentration around a circular hole is studied. The stress concentration factors (SCFs) obtained from both types of equations and the stress and couple-stress field on the boundary and inside the domain are compared with the analytical solutions.

In the third example, the classical problem of a Griffith crack is solved by means of the dual BEM with use of both types of BIEs. The BEM solution is compared with the analytical data for stress and couple-stress intensity factors, as well as for the displacements and microrotations at the crack faces.

In the fourth example, the dual BEM formulation is applied to the problem of an edge crack in both loading modes. The BEM results are shown to be in a good agreement with the analytical and FEM solutions available in the literature for the distribution of stresses and couple-stresses around the crack tip, as well as the stress and couple-stress intensity factors.

### (a) A square plate under pure bending

As a first example, a square micropolar plate under pure bending is considered (figure 3). The following boundary conditions are prescribed:
*e*_{i} of every component of the BEM solution *l* in the Bessel functions in the fundamental solutions. Therefore, in order to capture the behaviour of the kernels, the order of Gaussian quadrature was chosen depending on material parameter *l* (increasing for decreasing values of *l*), and it was kept the same for different values of *N* corresponding to the same value of *l*. The results for both types of BIEs are shown in table 1. The obtained relative errors, as defined by equation (5.5), are of order between 10^{−10} and 10^{−4}. It was observed that, for the same values of material parameters and the order of Gaussian quadrature, TBIEs perform slightly better than DBIEs.

In [6], the solution at the corner point was shown to have an accuracy of order 10^{−4} for the boundary discretization with eight to 20 elements, which is comparable with the errors obtained in this work for the singular BIEs with the use of four elements.

### (b) A circular hole in an infinite plate

In the second example, the problem of an infinite plate in tension, weakened by a circular hole of radius *a*, is considered (figure 5). The full analytical solution for all stresses and couple-stresses is given in [33]. In order to demonstrate the performance of the method for the non-straight boundaries as well as non-polynomial boundary conditions, the problem is modelled as a finite quarter-plate of size *L*×*L*, *L*=4*a* with the analytical tractions and couple tractions, given by rational functions, prescribed at *x*=*L* and *y*=*L* (figure 6).

In tables 2 and 3, the SCFs are presented for two values of the material length *l*:*l*=*a*, *l*=0.1*a*, respectively, *ν*=0.3 and for different values of the coupling number *N*, for both types of BIEs. In all cases, a mesh consisting of 68 elements, uniformly graded towards the edges of the hole, was used. The number of Gauss points per element was chosen for all cases to be 200. As it is seen from tables 2 and 3, the SCF solutions for a circular crack can be reproduced by both types of equations with the relative error within 1%. The relative error is defined as
*N*=1 corresponds to the couple-stress elasticity and it was studied also by means of the BEM in [20].

The stress distribution along the edges of the quarter plate is obtained as a direct output of BEM. In figure 7, distribution of *σ*_{θθ} along *θ*=*π*/2 is shown for *l*=*a*, *ν*=0.3 and various values of *N*. Next, the Somigliana stress formula (3.14) is used to obtain the stress distribution inside the domain. In figures 8–10, the distribution of normalized stresses *σ*_{rθ}, *σ*_{θr} and couple-stress *m*_{rz}, respectively, are shown along the line *θ*=*π*/4 for *l*=*a*, *ν*=0.3 and various values of *N*. As it can be seen from figures 7–10, the results are in a good agreement with the analytical solution from [33]. Analytical solutions are shown in all plots by solid black lines.

### (c) A Griffith crack

In the third example, we consider a straight crack in an infinite plane in a uniform tension (figure 11*a*) (known as a Griffith crack). The solution to this problem is a superposition of two solutions. The first one corresponds to an infinite uncracked plate in tension (figure 11*b*) and the second solution is the one of a crack opened up by a uniform tension applied to the crack faces (figure 11*c*). The solution to the first problem is given in [12] as
*c* for BEM modelling, i.e. the only BEM boundary corresponds to the crack faces

In order to avoid the degenerated system matrix due to the coincident collocation points, we prescribe DBIE on *Γ*_{+} and TBIE on *Γ*_{−}, an approach known as the dual BEM. The obtained results are compared with the analytical solutions obtained in [12]. Note, that the notations for *N* and *l* in [12] differ from the notations used here by the factor *l**, *N**. Therefore, *τ*=*l**/*c* and *M*=*N**/*τ*. The case *M*=1 and *ν*=0.25 is chosen for comparison. In table 4, the results are given in terms of four parameters. The first parameter is a value of the crack opening displacement at the centre of the crack:
*u*_{0} and *τ*=0.001 was used in BEM to approximate the limiting case of classical elasticity. In all cases, a fine mesh consisting of 196 elements, uniformly graded towards the crack tips, was used. The size of the mesh was chosen for the purpose of extracting stress intensity factors from the limiting values of the crack opening displacement as described below.

The stress intensity factor *K*_{t} and the couple-stress intensity factor *K*_{m} are defined as
*K*_{t} is directly expressed from the BEM solution for the crack opening displacement *u*_{2}(*r*,0) according to its asymptotic expansion in the vicinity of the crack tip:
*K*_{t} for various material parameters are compared with those in [12] and an agreement within 1% is shown.

Following the solution procedure in [12], the couple-stress intensity factor *K*_{m} is derived as
*K*_{m} for a Griffith crack is entirely defined by the integral of the crack opening displacement *u*_{2}(*r*,0), i.e. it can be obtained without using the asymptotics of the microrotation and couple-stress fields. This method allows to use coarser discretization of the crack domain to obtain accurate *K*_{m} for all values of material parameters in comparison with the methods, which require fitting solutions near the crack tip.

Values of

In figures 13 and 14, the full solutions for the crack opening displacement and microrotations are plotted for various values of material parameters, and an excellent agreement between the analytical solutions of [12] and the BEM data is seen.

Research in fracture mechanics of Cosserat materials (e.g. [10,13,34,35]) indicates that both displacements and microrotations in the vicinity of a crack tip have asymptotic expansions of order *ρ* is the distance to the crack tip. The asymptotic expansion of a microrotation in standard system of polar coordinates (*ρ*,*θ*), associated with the crack tip, is given in [13] as
*N*=0.9 and varying values of parameter *l*/*c* together with the plot of the first term of equation (5.16). Figure 15 illustrates the asymptotic property of microrotations, according to which the asymptotic range of equation (5.16) is
*l*. This behaviour indicates that in the extended FEM and BEM, the adequate choice of the enrichment zone must significantly depend on the material parameters. However, the application of enriched Cosserat BEM and FEM need detailed study and remain the subject of future work.

### (d) A plate with an edge crack

In this example, we consider the problem of a plate with an edge crack (figure 16) and compare our results with the FEM data obtained in [13]. For the sake of comparison, we consider the same dimensions of the plate and the loading conditions, as in [13], i.e. the width of the plate *W*=11 mm, height 2*H*=20 mm and the length of the crack *c*=1 mm. The centre of the coordinate system is placed at the crack tip and the polar coordinates (*r*,*θ*) are introduced. However, in order to avoid the rigid body motions, we impose a slightly different Dirichlet condition than the one in [13], i.e. since in the collocation BEM, the boundary condition cannot be imposed at the crack tip, we fix the point (10,0) on the boundary of the plate, i.e.
*σ*_{0}=100 MPa, i.e
*σ*_{0}=100 MPa, i.e.
*E*=100 GPa and *ν*=0.3.

In figures 17–20, we demonstrate the angular distributions of the stresses and couple-stresses (for an equivalent problem of a crack with the traction-free faces) defined as
*f*_{θθ}(0)=*g*_{θz}(0)=1 in Mode I and *f*_{θr}(0)=*g*_{rz}(0)=1 in Mode II. The remaining material parameters in figures 17–20 were taken *l*=0.025495 mm, *N*=0.849837, which correspond to parameters *γ*=100 MPa and *α*/*E*=1 in [13]. The equation (5.21) were evaluated at *r*=10^{−5} mm. In figures 17–20, the excellent agreement of the BEM data with the analytical solutions, derived in [13] is shown for both loading modes. Note, that due to the zero boundary condition imposed in [13] for *σ*_{rθ} instead of *σ*_{θr}, our BEM solution for *σ*_{rθ} corresponds to *σ*_{θr} in [13], and *σ*_{θr} corresponds to *σ*_{rθ} in [13]. Note, that according to [13], couple-stresses *m*_{rz}, *m*_{θz} do not have a *m*_{rz},*m*_{θz}, which, according to [13], for fixed *r*, have the analytical form:
*N*=0.849837 and the four values of the parameter *l*=0.025495, 0.254951, 0.806226,2.54951 corresponding to *γ*=10^{2}, 10^{4}, 10^{5}, 10^{6} in [13]. The results for both modes are given in table 6. The data in table 6 were obtained by fitting the displacements the crack faces according to
*μ*+*α*)(1−*ν*) instead of 2(*μ*−*α*)(1−*ν*) in eqn (72). The corrected constant is consistent with our BEM data.

For all material parameters, the same mesh with 193 elements on each crack face, gradually refined towards the crack tip was employed for the BEM analysis, and the data in the vicinity of *r*=10^{−3}*l* mm were used for fitting. The stress intensity factors in table 6 were normalized as follows:
*k*_{I}=208.1 MPa mm^{1/2}, *k*_{II}=199.3 MPa mm^{1/2} are the values corresponding to the Mode I and Mode II SIF solutions for an edge-cracked specimen in the case of the classical elasticity.

The couple-stress intensity factors were obtained by fitting the microrotations on the crack faces according to (5.16). The normalized couple-stress intensity factors are given by

## 6. Conclusion

This paper presented the derivation of the fundamental solutions for plane Cosserat elasticity and their application within a boundary element framework. Both singular and hypersingular BIEs of plane Cosserat elasticity were formulated using these fundamental solutions and solved by the BEM. The dual BEM, developed in this work, was shown to be an accurate numerical tool which can be used for analysis of problems with singularities, such as certain crack problems without limitations on the geometry or the type of the loading conditions. The excellent accuracy of the method was demonstrated on the classical example of a Griffith crack. The approach can be further used to model crack propagation in micropolar materials and can be extended to the three-dimensional case. The accuracy of the method for crack problems can be further improved, for a given number of degrees of freedom, by incorporating crack tip enrichments into the approximation space which can be derived from the asymptotic behaviour of the displacements and microrotations in the vicinity of a crack tip, also studied within this paper. This approach is known as the eXtendend BEM and is the subject of further study.

## Data accessibility

The analytical solutions, derived in this work, are available at https://sourceforge.net/projects/cosseratfundamentalsolutions/.

## Author' contributions

E.A. derived the fundamental solutions, implemented them into the boundary element method and obtained the numerical results. S.B. participated in the design of the study cases, analysis of the numerical results, structuring and proof-reading the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

The work presented in this paper was partially supported by Fondecyt Chile grant no. 11130259 entitled ‘Boundary element modelling of crack propagation in micropolar materials’. The time of the second author was funded by the European Research Council Starting Independent Research Grant (ERC StG) no. 279578 of the Framework Programme 7.

- Received April 1, 2015.
- Accepted June 5, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.