Simulation of cracks in a Cosserat medium using the extended finite element method

Abstract

Cosserat (micropolar, asymmetric) elasticity can better predict the mechanical behavior of the materials with a characteristic length scale than the classical theory of elasticity. However, the area of fracture modelling in a Cosserat medium is not widely presented in the literature. The simulation of cracks in a Cosserat medium using the eXtended Finite Element Method (XFEM) is presented in this paper. The proper crack tip enrichment of the translational and microrotational fields is important for the robustness and efficiency of the XFEM/Cosserat model. Using the example of an edge crack of the Mode I in this paper, we have shown that the values of the J-integral for the Cosserat solid are higher than those for the equivalent classical elastic solid. The difference becomes significant for the cases of materials with strong micropolar properties. The dependence of the fracture behavior of the crack of different sizes on the Cosserat elastic constants was found to be significant. Therefore, it is recommended that the special tip enrichment is considered for the microrotational field, and a careful analysis of the material parameters is performed in order to model fracture in a Cosserat medium.

Keywords

Cosserat elasticity fracture mechanics XFEM

1. Introduction

The simulation of the fracture behavior of materials with intrinsic length scales (such as grains, particles, etc.) remains a challenge to the computational mechanics community. However, the scale of the microstructure contributes significantly to the overall behavior and the strength of materials such as polymers, graphite, composites, cellular solid, bones, etc. The Cosserat elasticity theory, being a higher order gradient theory, effectively represents the scale dependent behavior of materials with significant microstructure [1, 2]. The Cosserat theory of elasticity was first proposed by the Cosserat brothers [3]. Later, works by Shaefer [4, 5], Eringen [6] and Nowacki [7] gave further development of the Cosserat theory, and the new name for the theory was introduced: the theory of micropolar or asymmetric elasticity. The Cosserat elasticity theory considers a medium as a continuous collection of particles that behave like rigid bodies having three independent microrotations at each point, in addition to three classical translational degrees of freedom. The isotropic Cosserat continuum has six elastic constants, in contrast to just two constants of the classical elastic continuum. The couple stress tensor is introduced in addition to the conventional stress tensor. This leads to an additional set of equilibrium equations that account for the moments. The essential feature of a micropolar continua is that the stress tensor is not necessarily symmetric, so that the balance of momentum equation should be modified accordingly. The presence of the rotational degrees of freedom allows us to explicitly model the mechanical behavior of a heterogeneous material with microstructure, while considering it as a continuum (see, for example, [8 –19] ). For example, the Cosserat theory is able to model discrete media, where interblock or interlayer slip occurs under loading, such as discontinuous rock masses and layered materials [20, 21]. In both cases, this type of deformation is independent of the macroscopic deformation of the system. Considering the bone structure, the mechanical properties of the bone are better described by the Cosserat theory than the classical theory of elasticity [22, 23]. The size effect in torsion tests of bone specimens was predicted by the Cosserat theory, according to Lakes [24] and Yang and Lakes [25].

The present paper focuses on a way to combine the eXtended Finite Element Method (XFEM), which is widely used in the area of the fracture mechanics, with the Cosserat elastic model. The XFEM was first introduced for the linear elastic fracture mechanics problems by Belytschko and Black [26] and Moes et al. [27]. Then, it was effectively applied for the modelling of different kinds of discontinuities in elastic solids such as various types of cracks [28 –30], interfaces [31], inclusions [32], dislocations [33 –35], and so on.

The XFEM is a numerical method that allows one to model the geometry of the discontinuity independently of the mesh. The geometry of the crack is convenient to define using the level set method [36]. The main idea of the XFEM is to introduce special enrichment functions to the displacement approximation, depending on the type of discontinuity [27]. For the crack problem, the elements containing the crack tip and some elements in the vicinity of the crack tip are enriched with the branch functions that span the asymptotic field. The Heaviside step function is used for the elements cut by the crack, but not enriched with the branch functions. This technique allows the modelling of the crack propagation without the need to re-mesh, which is an advantage over the Finite Element Method (FEM), where the crack surface must conform to the mesh. In the XFEM/Cosserat model presented here, both the translation and microrotation fields in the vicinity of the crack tip are enriched by the functions spanning the asymptotic fields [37] and by the Heaviside function away from the crack tip.

In the area of the Linear Elastic Fracture Mechanics (LEFM) for the Cosserat materials, very few numerical models have been presented, in contrast to experimental and analytical work. Amongst the first investigations in the field, the finite element analysis was performed to study the stress concentrations around a blunt crack in the Cosserat solid by Nakamura and Lakes [38]. The stress concentration factors were shown to be smaller than those of the classical problem. Simultaneously, the experimental studies of the cracks in the bones, considered as the Cosserat solids, were performed [22, 39]. Although, it falls out of the LEFM scope, it is worth mentioning that the wide implementation of the FEM/Cosserat model lies in the field of geomechanics, such as modelling of shear band formations and strain localizations. In Khoei and Karimi [40], a combined XFEM/Cosserat continuum model to simulate the strain localization in elasto-plastic solids was proposed. In Khoei et al. [41], the non-linear behavior of Cosserat materials, particularly plastic deformation, was modelled by representing material interfaces independent of the finite element mesh. Due to complexity of fracture modelling in Cosserat materials and the lack of information on the fracture parameters, the numerical simulations of the cracks in the Cosserat elastic medium are poorly covered up to date. In the present article, the XFEM model of the Mode I and II edge cracks in the Cosserat elastic solid are presented and discussed. The calculation of the energy release rate is performed. The obtained results could give a new rise to the modelling of the crack propagation problems in the Cosserat elastic solids.

This article is organized as follows. In Section 2, the concepts of the linear elastic fracture mechanic applicable to the crack modelling in the Cosserat materials are presented. In Section 3, the problem formulation and the governing equations are given. In Section 4, the extended finite element method formulation and the discrete equations are stated, and the numerical integration scheme is described. In Section 5, several numerical examples are presented. First two examples introduce the verification of the XFEM model and the convergence study for the Mode I and II edge cracks in the Cosserat medium, respectively. In the third example, the calculation of the energy release rate for the Mode I edge crack is presented and verified, and the comparison of the J-integral for the Cosserat and classical elasticity problems is performed for varying coupling levels. The conclusions are given in Section 6.

2. Modelling of fracture in Cosserat media

Crack tip behavior in classical linear elastic fracture mechanics is described by three crack opening modes, modes I–III [42]. Due to the presence of the additional degrees of freedom (microrotation) in Cosserat theory, crack tip behavior is described by three crack opening models plus three higher modes, bending modes IV–VI [43, 44]. These bending crack modes describe discontinuities in the microrotations and the moment stress concentrations at the crack tip. Fracture theory for Cosserat continua is still in its infancy. The J-integral (energy release rate) for a crack in a Cosserat media has six components, compared to three for classical elasticity [45]. While in classical linear elastic fracture mechanics there is a well-established relationship between the energy release rate and the stress intensity factors, $J = K_{I}^{2} / E^{'} + K_{II}^{2} / E^{'}$ , no such relationship exists for a Cosserat material. Furthermore, it is not yet clear what role modes IV–VI play in crack propagation, especially in mixed-mode fracture. To the authors’ knowledge, there is not a single empirical crack propagation law, backed by experimental data, for a Cosserat material based upon stress intensity factors and/or energy release rates. It is also not clear how the crack propagation direction should be chosen. In classical fracture mechanics, crack propagation direction is often chosen to be in the direction of the maximum circumferential stress which can be related in two-dimensions to the stress intensity factors; in a Cosserat continuum, it has not yet been determined what role the coupled stress (if any) should play in determining the crack propagation direction. Due to the limits on current knowledge about the relationships governing crack propagation, we will restrict ourselves to the simulation of static (non-propagating) cracks.

Diegele et al. [37] derived the form of the two-dimensional asymptotic fields for the displacements near a crack tip in a Cosserat medium:

\begin{matrix} (\begin{matrix} u_{r}^{tip} \\ u_{θ}^{tip} \end{matrix}) = \sqrt{\frac{r}{2 π}} \frac{{\tilde{K}}_{I}}{μ [μ + α (3 - 2 ν)]} \\ [\begin{matrix} \cos \frac{θ}{2} [μ (1 - 2 ν) - α + [μ + α (7 - 6 ν)] \sin^{2} \frac{θ}{2}] \\ - \sin \frac{θ}{2} [2 (μ + α) (1 - ν) - [μ + α (7 - 6 ν)] \cos^{2} \frac{θ}{2}] \end{matrix}] + \\ [\begin{matrix} O (r^{3 / 2}) \\ O (r^{3 / 2}) \end{matrix}], \end{matrix}

(1)

ϕ^{tip} = \sqrt{\frac{2 r}{π}} \frac{L_{I}}{γ + δ} \sin \frac{θ}{2} + O (r^{3 / 2}),

(2)

where $u_{r}^{tip}$ and $u_{θ}^{tip}$ are displacements at the crack tip and $ϕ^{tip}$ is the microrotation. Also, ${\tilde{K}}_{I}$ and $L_{I}$ are the Mode I stress intensity factors associated with the stress and couple stress, respectively. The stress intensity factors for the micropolar theory of elasticity ${\tilde{K}}_{I}$ and ${\tilde{K}}_{II}$ are distinct from the stress intensity factors of the classical theory of elasticity $K_{I}$ and $K_{II} .$ The deviation of ${\tilde{K}}_{I}$ from $K_{I}$ varies from 14% to 26% depending on the $α$ and $γ$ values [37].

The degree of coupling of the displacements and rotations is governed by elastic constant $α$ . For $α = 0$ and $α \to \infty$ the problem reduces to classical linear elasticity and to couple stress linear elasticity, respectively. Typical values for $α$ for real materials varies from $α = 10^{- 1}$ to $α = 10^{5}$ [16, 37].

It is noteworthy that ${\tilde{K}}_{I}$ and $L_{I}$ are not independent, and that the relationship between these quantities depends upon the elastic constants. This is consistent with earlier investigations where the effect of the coupled stresses on the crack tip stress concentration, energy release rate and mechanical behavior were discussed [42, 46].

2.1. Functions spanning the asymptotic solution

Careful manipulation of the asymptotic displacement solutions for an elastic crack in a Cosserat continuum (1) shows that the crack tip displacements are asymptotically spanned by the same four (branch) functions which span the Williams solutions of classical linear elastic fracture mechanics, i.e.

{F_{m}^{u}} = {\sqrt{r} \sin \frac{θ}{2}, \sqrt{r} \cos \frac{θ}{2}, \sqrt{r} \sin \frac{θ}{2} \sqrt{r} \sin θ, \sqrt{r} \cos \frac{θ}{2} \sqrt{r} \sin θ},

(3)

where $(r, θ)$ are the local polar coordinates centred at the crack tip.

From direction observation of the asymptotic microrotation fields for an elastic crack in a Cosserat continuum (2), it is clear that the first order term is spanned by:

{F^{ϕ}} = {\sqrt{r} \sin \frac{θ}{2}} .

(4)

The above sets of functions will be the basis for XFEM enrichment of the crack tip translation and microrotation fields.

2.2. The J-integral for a Cosserat continuum

The Cosserat equivalent to the J-integral has been derived in numerous ways [43 –45, 47]. Following the derivation presented by Lubarda and Markenscoff [45], the energy release rate (in the absence of the body forces and body moments) for plane strain micropolar elasticity is:

J_{i} = \int_{Γ} P_{ji} n_{j} d Γ,

(5)

where i and j are the Einstein indices; $Γ$ is the contour about the crack tip, and $n$ is the normal to the contour $Γ$ . $P$ is the energy momentum tensor [48]:

P_{ij} = W δ_{u j} - σ_{ik} u_{k,, j} - m_{i 3} ϕ_{3, j}

(6)

where $W$ is elastic energy per unit volume; $σ_{ik}$ and $m_{i 3}$ are the components of the stress and couple stress tensors, respectively; $u_{k,, j}$ and $ϕ_{3, j}$ are derivatives of the translations and rotations with respect to $x_{j}$ correspondingly. The elastic energy is defined as:

W = \frac{1}{2} σ_{ij} ε_{ij} + \frac{1}{2} m_{ij} ϰ_{ij},

(7)

where $ε_{ij}$ and $ϰ_{ij}$ are the components of the strain and curvature tensors, defined in the following section.

3. Problem formulation

Consider a linearly elastic domain $Ω$ bounded by $Γ$ which contains a traction free crack $Γ_{cr}$ , as shown in Figure 1. Boundary $Γ$ consists of the sets $Γ_{u},$ $Γ_{ϕ},$ $Γ_{t}$ and $Γ_{m},$ where displacements $u$ , microrotations $ϕ$ , surface forces (tractions) $t$ and surface couple forces (moments) $M$ are prescribed respectively. The boundary sets have the following properties:

Γ_{u} \cap Γ_{t} = ⊘, Γ_{u} \cup Γ_{t} = Γ, Γ_{ϕ} \cap Γ_{M} = ⊘, Γ_{ϕ} \cup Γ_{M} = Γ .

Figure 1.

Body with external and internal boundaries subjected to loads.

In two dimensions, the governing equilibrium equations are

\nabla \cdot σ + b^{σ} = 0,

(8)

\nabla \cdot m + b^{m} + \tilde{σ} = 0,

(9)

where $σ$ is the Cauchy stress, $b^{σ}$ is the external forces, $m$ is the couple force (stress), $b^{m}$ is the external moment, and $\tilde{σ} = σ_{21} - σ_{12},$ which are subject to the following boundary conditions

σ \cdot n = \bar{t} on Γ_{t}, m \cdot n = \bar{M} on Γ_{M},

(10)

u = \bar{u} on Γ_{u}, ϕ = \bar{ϕ} on Γ_{ϕ},

(11)

σ \cdot n_{cr} = 0 and m \cdot n_{cr} = 0 on Γ_{cr},

(12)

where $n$ is the unit outward normal to $Γ$ and $n_{cr}$ is the normal to the crack surface $Γ_{cr} .$

The translations $u$ and the microrotations $ϕ$ are taken to be the independent variables. The strain $ε$ and the curvature $ϰ$ have the following relations with $u$ and $ϕ :$

ε = \nabla u - \tilde{ϕ},

(13)

ϰ = \nabla ϕ,

(14)

where

\tilde{ϕ} = [\begin{matrix} 0 & - ϕ_{3} \\ ϕ_{3} & 0 \end{matrix}] .

The linear constitutive equations for the coupled model are:

σ = D^{u} : ε,

(15)

m = D^{ϕ} : χ,

(16)

where $D^{u}$ and $D^{ϕ}$ are the elastic moduli tensors which for an isotropic material are defined as:

D^{u} = [\begin{matrix} 2 μ + λ & λ & 0 & 0 \\ λ & 2 μ + λ & 0 & 0 \\ 0 & 0 & μ - α & μ + α \\ 0 & 0 & μ + α & μ - α \end{matrix}]

(17)

and

D^{ϕ} = [\begin{matrix} γ + κ & 0 \\ 0 & γ + κ \end{matrix}],

(18)

where $μ,$ $λ, α,$ $γ$ and $κ$ are material constants.

3.1. Weak form

Following the standard Galerkin procedure, the following weak form of the coupled boundary value problem governing an elastic Cosserat medium can be derived.

Find $u \in U$ and $ϕ \in Φ$ such that:

\int_{Ω} \nabla w : σ d Ω = \int_{Ω} w \cdot b^{σ} d Ω + \int_{Γ} w \cdot \bar{t} d Γ, \forall w \in U_{0},

(19)

\int_{Ω} \nabla v : m d Ω = \int_{Ω} v \cdot b^{m} d Ω + \int_{Γ} v \cdot \bar{M} d Γ, \forall v \in Φ_{0},

(20)

where

U = {u \in H^{1}, u = \bar{u} on Γ_{u}, u is discontinuous on Γ_{cr}},

(21)

Φ = {ϕ \in H^{1}, u = \bar{ϕ} on Γ_{ϕ}, ϕ is discontinuous on Γ_{cr}},

(22)

and

U_{0} = {w \in H^{1}, w = 0 on Γ_{u}}

(23)

Φ_{0} = {v \in H^{1}, v = 0 on Γ_{ϕ}}

(24)

where $H^{1}$ is the corresponding Sobolev space.

4. The eXtended Finite Element Method

Consider a crack $Γ_{cr}$ passing through a discretized domain, as shown in Figure 2. Following Lakes [24], the geometry of the crack is represented by two orthogonal level sets $f (x)$ and $g (x)$ . The set $f (x)$ is defined as the signed distance to the crack surface, and so the crack surface is given by points on the zero level set, i.e. $f (x) = 0$ . $g (x)$ is defined as the signed distance to the crack front, such that $\nabla f \cdot \nabla g = 0 .$ The sign of $g (x)$ is chosen such that for all points on the crack surface $g (x) \geq 0$ . The crack surface $Γ_{cr}$ is therefore defined at all points where the following conditions are true: $f (x) = 0$ and $g (x) \geq 0$ .

Figure 2.

Level set method for the crack positioning.

Following the XFEM methodology developed by Belytschko and Black [26] and Moes et al. [27], the FEM displacement approximation is enriched by the Heaviside step function at nodes with supports cut by the crack (denoted by circles in Figure 3) and by a set of tip functions which span the asymptotic crack tip fields for nodes within a given domain about the crack tip (denoted by squares in Figure 3).

Figure 3.

Illustration of a two dimensional domain containing a crack $cr$ . The circled nodes are enriched by the Heaviside function, and the squared nodes are enriched by the crack tip enrichment function.

The XFEM approximation for the translation and microrotation for a crack in a linearly elastic Cosserat continuum are

\begin{matrix} u^{h} (x) = \sum_{I \in S} N_{I} (x) d_{I}^{u} + \sum_{J \in S_{cr}} N_{J} (x) [H (f (x)) - H (f (x_{J}))] a_{J}^{u} + \\ + \sum_{K \in S_{tip}} N_{K} (x) \sum_{m = 1}^{M} [F_{m}^{u} (r, θ) - F_{m}^{u} (r_{n}, θ_{n})] b_{mK}^{u}, \forall x \in Ω, \end{matrix}

(25)

\begin{matrix} ϕ^{h} (x) = \sum_{I \in S} N_{I}^{ϕ} (x) d_{I}^{ϕ} + \sum_{J \in S_{cr}} N_{J}^{ϕ} (x) [H (f (x)) - H (f (x_{J}))] a_{J}^{ϕ} + \\ + \sum_{K \in S_{tip}} N_{K}^{ϕ} (x) [F^{ϕ} (r, θ) - F^{ϕ} (r_{n}, θ_{n})] b_{K}^{ϕ}, \forall x \in Ω, \end{matrix}

(26)

where $S$ is the set of all nodes, $S_{cr}$ is the set of nodes with support crossed by the crack surface, $S_{tip}$ is the set of nodes inside the predefined area around the crack tip, $N_{I}$ are the standard finite element shape functions for the vector field (translations), $N_{I}^{ϕ}$ are the standard finite element shape functions for the microrotations, $d_{I}^{u}$ and $d_{I}^{ϕ}$ are the standard FEM nodal translational and rotational degrees of freedom, $a_{J}^{u}$ and $a_{J}^{ϕ}$ are the Heaviside function enriched degrees of freedom, and $b_{K}^{u}$ and $b_{K}^{ϕ}$ are the branch function enriched degrees of freedom. The sets $S_{cr}$ and $S_{t ip}$ are mutually exclusive, as illustrated in Figure 3.

The nodes in $S_{t ip}$ are identified using the level sets, i.e.

S_{tip} = {K | f^{2} (x_{K}) + g^{2} (x_{K}) \leq r_{tip}^{2}},

where $r_{tip}$ is predefined radius of a domain around the crack tip to be enriched with special tip enrichment functions. The nodes in $S_{cr}$ are those with supports is cut by $Γ_{cr}$ and is not in $S_{t ip}$ .

The Heaviside step function is defined in the usual sense as

H (z) = {\begin{matrix} 1 if z \geq 0 \\ 0 otherwise, \end{matrix}

(27)

The enrichment functions $F_{m}$ and $F_{m}^{ϕ}$ are the basis functions (3) and (4) that span the first-order terms of the near tip asymptotic translation and microrotation fields, respectively. As discussed in Section 2, the asymptotic crack tip displacements for a Cosserat continuum are spanned by the same set of functions as the asymptotic crack tip solutions for classical elasticity. Therefore, the crack tip displacements will be enriched by the standard XFEM branch functions. The asymptotic crack tip microrotation solution is spanned by a single function (4), which is used to enrich the microrotations in the vicinity of the crack tip.

4.1. Numerical integration

In order to achieve high quality results and optimal convergence rates, accurate numerical integration of the discrete equations is required since the gradient of the tip enrichment functions are singular at the crack tip and the Heaviside enrichment function is discontinuous across the crack surface [49 –52]. Integration over the unenriced elements, Heaviside enriched elements and the crack tip element is done by second-order Gauss quadrature rules, element subdivision combined with the second-order Gauss quadrature [27], and element subdivision combined with 15th-order almost polar integration Gauss quadrature [50]. Lastly, for the elements that are enriched by the branch functions but do not contain the crack tip, an 8th-order Gauss quadrature rule is used.

According to Khoei et al. [41], several aspects are crucial in order to obtain accurate results and optimal convergence rates in the XFEM modelling:

The proper integration rules in the elements cut by a crack and containing a crack tip should be considered.

The geometrical enrichment is crucial to obtain optimal convergence rates. It allows us to get the optimal convergence rate in oppose to the topological enrichment, when only particular number of nodes layers are enriched.

Blending elements allow for the smooth transition between enriched and un-enriched subdomains by using the special weight function with compact support.

4.2. Blending

Another important issue for the accuracy of the solution in the XFEM is the blending elements [52, 54, 55]. The blending elements are those that have some of the nodes enriched and some of them un-enriched. The standard enrichment can add parasitic terms to the approximation space of the blending elements, i.e. the reproducing property of the partition of unity is lost. No difficulty occurs in the blending of the shifted Heaviside enrichment, because the enrichment is only considered for the elements cut by the crack. Errors due to blending of the branch function for this problem are negligible and do not affect the convergence rate of the XFEM. Therefore, no blending treatment has been adopted in this work.

4.3. Crack tip enrichment domain

To obtain optimal convergence rates with the XFEM, the tip enrichment domain around the crack tip must be a constant (i.e. independent of the element size); this enrichment scheme is called the geometrical enrichment [49, 50]. This enrichment scheme is used in all convergence studies presented in this article.

5. Examples

In this section, three numerical examples are presented to demonstrate the XFEM simulation of the cracks in the Cosserat elastic solid. The first and second examples introduce the solutions of the Mode I and II edge crack problems, respectively. The patch tests and the convergence properties of the method are demonstrated for both cases. The third example presents the numerical calculation of the energy release rate for the edge crack problem under tension. The J-integral is calculated for the several crack sizes, considering varying Cosserat elastic parameters, and compared to the J-integral of the classical elasticity.

5.1. Mode I edge crack

A square plate of size $L$ by $L$ is considered with an edge crack of length $a$ as shown in Figure 4, where $L = 100$ mm and $a = L / 2 .$ The domain is fixed at the bottom, and the displacement $u_{2}$ = 5 mm in the $x_{2}$ direction is prescribed along the top edge. The microrotations $ϕ$ = 0 are prescribed along the bottom and top edges of the domain, which allows for the microrotations to emerge around the crack tip. The crack surface is traction-free. A plain strain condition is assumed. The material properties are: Young’s modulus $E = 10, 000$ MPa, Poisson’s ratio $ν$ = 0.3, $λ = 5796$ ${N / mm}^{2}$ , the elastic constants $μ = 3846$ ${N / mm}^{2}$ , $κ = 81$ ${N / mm}^{2}$ , $γ = 1000$ N, $β = 0$ . As mentioned $α = 0$ reduces the problem to the one of the classical elasticity. For this and the following examples, the values of $α$ were chosen to be within the range [ $10^{- 1},$ $10^{4}]$ ${N / mm}^{2}$ to be relevant with the provisions of the Cosserat elasticity. Particularly for the patch tests and the convergence studies in the first and second examples, $α$ was taken to be 103 ${N / mm}^{2}$ .

Figure 4.

The domain with an edge crack under tension.

In order to perform the convergence study for the Mode I and II edge crack problems, the normalized energy norm is employed as a convergence criterion. The employed mesh sizes are 11 $\times$ 11, 21 $\times$ 21, 31 $\times$ 31, 41 $\times$ 41, 51 $\times$ 51, 71 $\times$ 71 elements.

The energy norm is defined as

∥ u_{ex} - u_{h} ∥_{E} = (0.5 \int_{Ω} {(ε_{ex} - ε_{h})}^{T} D^{u} (ε_{ex} - ε_{h}) d Ω + 0.5 \int_{Ω} {(χ_{ex} - χ_{h})}^{T} D^{ϕ} (ε_{ex} - ε_{h}) d Ω^{1 / 2},

(28)

where $u,$ $ε$ and $χ$ with the subscript “ex” are the exact solutions of the displacement, strain and curvature, and with the superscript “h” are the corresponding numerical solutions.

The energy norm is normalized as:

\frac{∥ u_{ex} - u_{h} ∥_{E}}{∥ u_{ex} ∥_{E}},

(29)

The solution of the addressed Mode I edge crack problem by the FEM with the mesh of 40,000 elements was taken to be the exact solution. Two types of the tip enrichment zones were considered for the convergence study of the XFEM solution: the topological tip enrichment of one layer of nodes around the crack tip and the geometrical tip enrichment of the size $L / 10$ . Figure 5 shows three convergence curves and the corresponding convergence rates for the Mode I edge crack solution in the Cosserat medium.

Figure 5.

Convergence study for the Mode I edge crack problem.

As expected, the optimal convergence rate of approximately 1 in the energy norm for the XFEM with geometrical enrichment is obtained. The suboptimal convergence rates of 0.5 in energy norm are obtained for the XFEM solutions with topological tip enrichment and for the FEM solution. The geometrical tip enrichment of the size $L / 10$ was found to be the optimal choice between the accuracy of the solution and the computational time.

5.2. Mode II edge crack

In the second example, a Mode II edge crack problem is studied. Consider a square plate of size $L \times L$ with an edge crack of length $a$ , as shown in Figure 6. The length of the plate is $L = 100$ mm, and the crack size is $a = L / 2$ . The domain is fixed at the bottom. The displacement $u_{1} = 5$ mm in $x_{1}$ direction is prescribed along the top edge. The microrotations are prescribed to be $ϕ = 0$ along the bottom and top edges of the domain. The crack surface is traction free. A plain strain condition is assumed. The material constants are taken to be the same as in the previous example of the Mode I edge crack.

Figure 6.

The domain with an edge crack under shear.

The error estimation procedure, described in the first example, is performed for the convergence study for the Mode II edge crack problem (Figure 7). Again, the optimal convergence rate of 1 in energy norm was obtained for the XFEM solution with the geometrical tip enrichment, and the suboptimal convergence rate of approximately 0.5 in energy norm was obtained for the XFEM solution with the topological tip enrichment and the FEM solution. The slight decrease in the convergence rates from the optimal values is observed, but expected. In order to get higher convergence rates, the very fine mesh should be considered for the exact FEM solution. Although, the computational time will increase drastically.

Figure 7.

Convergence study for the Mode II edge crack problem.

Concluding, the results of the convergence studies for the cracks of Modes I and II in the Cosserat medium allows us to rely on the developed XFEM model with confidence.

5.3. J-integral computation

5.3.1. Verification

The calculation of the J-integral for the Cosserat elastic solid is performed using the contour integral (5). The most important concept of the contour integral—path-independence under the static condition—is satisfied in our calculations. Considering the contour as a circle, the obtained values of the J-integral are constant for the radii greater than $W / 10$ with its centre at the crack tip, where $W$ is the width of the domain. Therefore, for the following calculations, the radius of the contour integral is assumed to be $W / 2.2$ to capture the cracks of different sizes. The J-integral is zero on the contour without the discontinuity, which is satisfied by the presented numerical calculation. The number of sectors is taken to be 1000, based on the test of the J-integral invariance versus the number of sectors on the contour. And finally, the J-integral converges with the mesh refinement, as expected.

In order to verify the calculation of the J-integral, two specific tests were proposed. In the first test, the Cosserat problem was fully reduced to the classical case and the calculated J-integral was compared to the one of the corresponding analytical solution of the classical J-integral. Assuming the plane strain condition, the classical J-integral is calculated as

J_{1} = K_{I}^{2} (1 - ν^{2}) / E .

For the test, the Mode I edge crack in the domain with the prescribed analytical solution on the boundaries [49] was considered. The numerical values of the calculated J-integral was found to be less than 1% different from the analytical values.

For the second test, Mode I and II edge crack were considered in the Cosserat elastic medium with $α = 1000$ ${N / mm}^{2}$ . The J-integral can be analytically calculated as the change in the energy of the domain due to the crack length increment, i.e $Δ U / Δ a$ . Considering the mesh of 101 $\times$ 101 elements with the crack increment $Δ a$ = 3, for the Mode I crack case the deviation from the analytical value was of 4%, and for the Mode II crack case the deviation was of 5%. Taking into account the fact that this approximation is valid for the crack increment approaching zero $Δ a \to$ 0, the obtained results are accurate enough to consider the calculation of the J-integral in the Cosserat elastic medium to be correct.

The obtained results verify the most important characteristics of the contour J-integral, which allow us to proceed to the study of the energy release rate for the Cosserat medium.

5.4. Comparison of the J-integral for the Cosserat and classical elastic solids

Consider the domain illustrated in Figure 8. A plate of height $L$ and width $W$ contains an edge crack of length $a$ , where $L / W = 2$ , $W = 100$ mm, and $a$ has varying value. The bottom edge of the domain is fixed. The stress $σ = 3$ MPa is applied at the top edge of the domain. The Young’s modulus is $E = 10, 000$ MPa, the Poisson’s ratio is $ν$ = 0.3, and the material constants are: $λ$ = 5796 ${N / mm}^{2}$ , $μ$ = 3846 ${N / mm}^{2}$ , $μ$ = 81 ${N / mm}^{2}$ , $γ$ = 1000 N and $β$ = 0.

Figure 8.

The domain with an edge crack under tension loading.

In this example, the effect of the crack size on the J-integral is analyzed. Several cases are considered depending on the values of the Cosserat elastic constant, which influences the level of assymetricity of the stress and strain tensors. For this study, parameter $α$ is chosen to be $10^{- 1},$ $10^{2}$ , $10^{3}$ , $10^{4}$ , which is attributed to the Cosserat elastic theory. The mesh size is 101 $\times$ 101 elements. The elements to be enriched by the asymptotic field are chosen within the geometrical tip enrichment zone of size $W / 10$ .

The numerically calculated J-integral for the Cosserat elasticity is compared to the one of the classical elasticity. The exact stress intensity factor solution of the addressed problem for the classical elasticity is given by Pasternak et al. [44]:

K_{I} = C σ \sqrt{a π},

(30)

where $C$ is a geometry correction factor.

For $a / W \leq 0.6$ , the factor $C$ is defined as:

C = 1.122 - 0.231 (\frac{a}{W}) + 10.55 {(\frac{a}{W})}^{2} - 21.72 {(\frac{a}{W})}^{3} + 30.382 {(\frac{a}{W})}^{4} .

(31)

The reference J-integral of the classical elasticity for the Mode I crack under the plane strain condition is related to the stress intensity factors as follows:

J_{I}^{ref} = \frac{K_{I}^{2}}{E} (1 - ν^{2}) .

(32)

Figure 9 summarizes the J-integral calculations versus the crack size for the Cosserat problem, parametrized by $α$ , and the reference classical J-integral. Figure 10 illustrates the normalized J-integral of the Cosserat elasticity, which is calculated as:

\frac{| J_{I}^{h} - J_{I}^{ref} |}{J_{I}^{ref}} \cdot 100 %,

(33)

where $J_{I}^{h}$ is the J-integral for the Cosserat problem and $J_{I}^{ref}$ is the J-integral for the corresponding classical elasticity problem (32).

Figure 9.

The J-integral versus the crack length.

Figure 10.

The normalized J-integral with respect to the crack length.

It can be seen that the behavior of the J-integral for the Cosserat elasticity is similar to that of the classical elasticity for the $α \in 10^{- 1}, 10^{2}]$ ${N / mm}^{2}$ . However, once the value of the parameter $α$ has passed $10^{2}$ ${N / mm}^{2}$ , we can see the difference in the values for the J-integrals already at the level of 5% and further with the increase of $α$ one can observe a very significant difference in the values of the J-integrals for Cosserat and the classical theory, especially when the value of $α$ achieves the magnitude of $10^{4}$ ${N / mm}^{2}$ . At this stage, the value of the J-integral computed under the provisions of the micropolar theory can be up to 40% higher than the results for the J-integral computed under the classical assumptions.

The difference of up to 40 % between the values of the J-integrals when the material constant $α$ passes through the magnitude of $10^{4}$ ${N / mm}^{2}$ can be attributed to the fact that at this range of $α$ the micropolar properties of the material are the strongest. Lakes provides the table of values of the micropolar elastic constants corresponding to the materials with strong micropolar properties [16]. He introduces the parameter

Ψ = \frac{β + γ}{β + γ + α},

(34)

which is called the polar ratio. $β,$ $γ$ and $α$ from (34) are micropolar elastic constants. The measured value for $Ψ$ in Lakes’s work [16] varies from 0 to 1.5 for such materials as foams, graphite and human bones. The interval from 0 to 1.5 for $Ψ$ corresponds to the interval from $10^{3}$ to $10^{5}$ ${N / mm}^{2}$ for $α$ . Therefore, as expected on the basis of the physical assumptions, at these values of the parameter $α$ we indeed observe the largest difference in the values of the J-integrals. The results obtained herein are in good agreement obtained for stress concentrations around a crack in micropolar elasticity obtained by Shmoylova et al. [18], for torsion of micropolar bars obtained by Potapenko and Shmoylova [19] and by Potapenko et al. [56]. In addition, Mindlin [57], Weitsmann [58] and Hartranft and Sih [59] obtained approximately the same difference in the values of stresses between the approach engaging the provisions of Cosserat elasticity and the classical approach when they considered stress concentrations around rigid inclusions in micropolar elasticity.

Based on the results of these studies, the J-integral presented by Atkinson and Leppington [47] is valid for the calculations of the energy release rate in the Cosserat medium. In general, the higher values of J-integral are obtained for the problems with the microrotations involved than the J-integral of the classical elasticity. This means that, under the same boundary conditions and loading mode, a crack grows faster in the Cosserat elastic solid than in the classical elastic solid. For the higher values of the Cosserat parameter, the crack will grow up to 40% faster in the Cosserat material.

6. Conclusions

The simulation of the edge cracks of Modes I and II in a Cosserat medium using the XFEM was presented. The numerical calculation of the path-independent J-integral for the Cosserat elastic medium was introduced. The Cosserat J-integral dependence on the crack size was studied for the parameter $α$ in the range from $10^{- 1}$ to $10^{5}$ ${N / mm}^{2}$ , which is known to affect the coupling level for the Cosserat elasticity, and was compared to the classical J-integral for the analogous problem. It was discovered that values of the J-integral for the Cosserat medium can be up to 40% higher than those for the classical elastic medium for the values of $α$ corresponding to such materials as human bones and graphite. The results obtained herein are consistent with the previous results published elsewhere [18, 19, 56 –59]. This investigation gives more evidence that material microstructure really does have a significant effect on the elastic response of materials.

Based on the results obtained from the convergence studies, it is recommended that the proper branch functions for the translational and rotational degrees of freedom and geometrical tip enrichment zone are considered for the modelling of cracks in Cosserat solids in order to obtain the optimal convergence rates for both Mode I and II cracks. For the J-integral calculation the elastic constants should be carefully measured and evaluated while modelling the cracks in a Cosserat elastic media, as they affect the fracture behavior of the material.

Footnotes

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Tomita

Yuan

. Effective properties of Cosserat composites with periodic. Mech Res Commun 2001; 28: 265–270.

Green

Naghdi

. A thermomechanical theory of a Cosserat point with application to composite materials. Q J Mech Appl Math 1991; 44: 335–355.

Cosserat

. Theorie des corps deformables. Paris: A Hermann, 1909.

Schaefer

. Das Cosserat kontinuum. ZAMM—Z Angew Math Me 1967; 47: 485–498 (1967).

Schaefer

. Analysis der motorfelder im Cosserat-kontinuum. ZAMM—Z Angew Math Me 1967; 47: 319–328.

Eringen

. Linear theory of micropolar elasticity. J Math Mech 1966; 15: 909–923.

Nowacki

. Theory of micropolar elasticity. Warsaw: Polish Scientific Publishers, 1972.

Potapenko

Schiavone

Mioduchowski

. Anti-plane shear deformations in a linear theory of elasticity with microstructure. Z Angew Math Phys 2005; 56: 516–528.

Potapenko

Schiavone

Mioduchowski

. On the solution of mixed problems in anti-plane micropolar elasticity. Math Mech Sol 2003; 8: 151–160.

10.

Potapenko

. A generalized Fourier approximation in anti-plane micropolar elasticity. J Elasticity 2005; 81: 159–177.

11.

Potapenko

. Fundamental sequences of functions in the approximation of solutions to mixed boundary-value problems of Cosserat elasticity. Acta Mechanica 2005; 177: 61–70.

12.

Schiavone

. On existence theorems in the theory of extensional motion of thin micropolar plates. Int J Eng Sci 1989; 27: 1129–1133.

13.

Schiavone

. Integral equation methods in plane asymmetric elasticity. J Elasticity 1996; 43: 31–43.

14.

Shmoylova

Potapenko

Rothenburg

. Weak solutions of the interior boundary value problems of plane Cosserat elasticity. Z Angew Math Phys 2006; 57: 506–522.

15.

Shmoylova

Potapenko

Rothenburg

. Stress distribution around a crack in plane micropolar elasticity. J Elasticity 2007; 86: 19–39.

16.

Lakes

. Experimental methods for study of Cosserat elastic solids and other generalized elastic continua. In: Muhlhaus

(ed) Continuum Models for Materials with Microstructure. Chichester, UK; New York: John Wiley & Sons, 1995.

17.

Shmoylova

Potapenko

. Weak solutions to anti-plane boundary value problems in a linear theory of elasticity with microstructure. J Elasticity (accepted for publication).

18.

Shmoylova

Potapenko

Rothenburg

. Boundary element analysis of stress distribution around a crack in plane micropolar elasticity. Int J Eng Sci 2007; 45: 199–209.

19.

Potapenko

Shmoylova

. Weak solutions of the problem of torsion of micropolar elastic beams. Z Angew Math Phys 2010; 61: 529–536.

20.

Shi

. Discontinuous deformation analysis: A new numerical model for the statics and dynamics of block systems. PhD Thesis, University of California, Berkeley, CA, 1988.

21.

Riahi

Curran

. Full 3D finite element Cosserat formulation with application in layered structures. Appl Math Model 2009; 33: 3450–3464.

22.

Park

Lakes

. Cosserat micromechanics of human bone: Strain redistribution by a hydration sensitive constituent. J Biomech 1986; 19: 385–397.

23.

Lakes

Nakamura

Behiri

. Fracture mechanics of bone with short cracks. J Biomech 1990; 23: 967–975.

24.

Lakes

. Size effects and micromechanics of a porous solid. J Material Sci 1983; 18: 2572–2580.

25.

Yang

JFC

Lakes

. Experimental study of micropolar and couple stress elasticity in compact bone in bending. J Biomech 1982; 15: 91–98.

26.

Belytschko

Black

. Elastic crack growth in finite elements with minimal remeshing. Int J Num Method Eng 1999; 45: 601–620.

27.

Moes

Dolbow

Belytschko

. A finite element method for crack growth without remeshing. Int J Num Method Eng 1999; 46: 131–150.

28.

Sukumar

Moes

Moran

. Extended finite element method for three-dimensional crack modelling. Int J Num Method Eng 2000; 48: 1549–1570.

29.

Belytschko

Moes

Usui

. Arbitrary discontinuities in finite elements. Int J Num Method Eng 2001; 50: 993–1013.

30.

Sukumar

Belytschko

. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Num Method Eng 2000; 48: 1741–1760.

31.

Belytschko

Gracie

. On XFEM applications to dislocations and interfaces. Int J Plasticity, 2007; 23: 1721–1738.

32.

Sukumar

Chopp

Moes

. Modeling holes and inclusions by level sets in the extended finite element method. Comp Meth Appl Mech Eng 2001; 190: 6183–6200.

33.

Gracie

Ventura

Belytschko

. A new fast method for dislocations based on interior discontinuities. Int J Num Method Eng 2007; 69: 423–441.

34.

Belytschko

Gracie

. On XFEM applications to dislocations in problems with interfaces. Int J Plasticity 2007; 23: 1721–1738.

35.

Gracie

Belytschko

. Concurrently coupled atomistic and XFEM models for dislocations and cracks. Int J Num Method Eng 2009; 78: 354–378.

36.

Stolarska

Chopp

Moes

. Modelling crack growth by level sets in the extended finite element method. Int J Num Method Eng 2001; 51: 943–960.

37.

Diegele

Elsasser

Tsakmakis

. Linear micropolar elastic crack-tip fields under mixed mode loading conditions. Int J Fracture 2004; 129: 309–339.

38.

Nakamura

Lakes

. Finite element analysis of stress concentration around a blunt crack in a Cosserat elastic solid. Comp Meth Appl Mech Eng 1988; 66: 257–266.

39.

Lakes

Yang

JFC

. Micropolar elasticity in bone: Rotation modulus. In: Proceedings of the 18th Midwest mechanics conference, Iowa City, IA, 1983, Vol. 12, pp.239–242.

40.

Khoei

Karimi

. An enriched FEM model for simulation of localization phenomenon in Cosserat continuum theory. Comput Mater Sci 2008; 44: 733–749.

41.

Khoei

Anahid

Shahim

. An extended arbitrary Lagrangian–Eulerian finite element method for large deformation of solid mechanics. Finite Elem Anal Des 2008; 44: 401–416.

42.

Atkinson

Leppington

. The effect of couple stresses on the tip of a crack. Int J Sol Struct 1977; 13: 1103–1122.

43.

Muhlhaus

Pasternak

. Path independent integrals for Cosserat continua and application to crack problems. Int J Fracture 2002; 113: 21–26.

44.

Pasternak

Dyskin

Muhlhaus

. Cracks of higher modes in Cosserat continua. Int J Fracture 2006; 140: 189–199.

45.

Lubarda

Markenscoff

. On conservation integrals in micropolar elasticity. Philosoph Mag 2003; 83: 1365–1377.

46.

Sternberg

Muki

. The effect of couple-stresses on the stress concentration around a crack. Int J Sol Struct 1967; 3: 69–95.

47.

Atkinson

Leppington

. Some calculations of the energy-release rate for cracks in micropolar and couple-stress elastic media. Int J Fracture 1974; 10: 599–602.

48.

Eshelby

. The elastic energy-momentum tensor. J Elasticity 1975; 5: 321–335.

49.

Bechet

Minnebo

Moes

. Improved implementation and robustness study of the X-FEM for stress analysis around cracks. Int J Num Method Eng 2005; 64: 1033–1056.

50.

Laborde

Pommier

Renard

. High-order extended finite element method for cracked domains. Int J Num Method Eng 2005; 64(3): 354–381.

51.

Ventura

. On the elimination of quadrature subcells for discontinuous functions in the extended finite element method. Int J Num Method Eng, 2006; 66: 761–795.

52.

Ventura

Gracie

Belytschko

. Fast integration and weight function blending in the extended finite element method. Int J Num Method Eng 2009; 77: 1–29.

53.

Gracie

Wang

Belytschko

. Blending in the extended finite element method by discontinuous Galerkin and assumed strain methods. Int J Num Method Eng 2008; 74: 1645–1669 (2008).

54.

Fries

. A corrected XFEM approximation without problems in blending elements. Int J Num Method Eng 2008; 75: 503–532.

55.

Tada

Paris

Irwin

. The Stress Analysis of Cracks Handbook. 3rd edition. New York: ASME Press, 2000.

56.

Potapenko

Schiavone

Mioduchowski

. Generalized Fourier series solution of torsion of an elliptic beam with microstructure. Appl Math Let 2004; 17: 189–192.

57.

Mindlin

. Effect of couple stresses on stress concentrations. Exp Mech 1963; 3: 1–7.

58.

Weitsmann

. Couple-stress effects on stress concentration around a cylindrical inclusion in a field of uniaxial tension. J Appl Mech 1965; 6: 424–427.

59.

Hartranft

Sih

. The effect of couple-stresses on the stress concentration of a circular inclusion. J Appl Mech 1965; 32: 429–431.