# Spherical Indentation of a Micropolar Solid: A Numerical Investigation Using the Local Point Interpolation–Boundary Element Method

Gaël Pierson, M'barek Taghite, Pierre Bravetti, Richard Kouitat Njiwa

► To cite this version:

Gaël Pierson, M'barek Taghite, Pierre Bravetti, Richard Kouitat Njiwa. Spherical Indentation of a Micropolar Solid: A Numerical Investigation Using the Local Point Interpolation–Boundary Element Method. Applied Mechanics, 2021, 2, pp.581 - 590. 10.3390/applmech2030033. hal-04055648

HAL Id: hal-04055648

https://hal.science/hal-04055648

Submitted on 3 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Spherical Indentation of a Micropolar Solid: A Numerical Investigation Using the Local Point Interpolation-Boundary Element Method

Gaël Pierson, M'Barek Taghite, Pierre Bravetti and Richard Kouitat Njiwa *

Institut Jean Lamour-CNRS UMR 7198—Université de Lorraine, 2 Allée André Guinier-BP 50840, CEDEX, 54011 Nancy, France; gael.pierson@univ-lorraine.fr (G.P.); mbarek.taghite@univ-lorraine.fr (M.T.); pierre.bravetti@univ-lorraine.fr (P.B.)

* Correspondence: richard.kouitat@univ-lorraine.fr

**Abstract:** The load-penetration curve in elastic nanoindentation of an elastic micropolar flat by a diamond spherical punch is analyzed. The presented results are obtained by a specifically developed numerical tool based on a judicious combination of the conventional boundary element method and strong form local point interpolation method. The results show that the usual linear relationship between the material depression and the square of the radius of the contact area is also valid in this case of micropolar elastic material. It is also shown that the relation between the indentation stress (applied load over the contact surface) and the indentation strain (ratio of contact radius by the punch radius) is linear. The proportionality coefficient which is none other than the indentation stiffness varies with the coupling factor of the micropolar elastic medium. A relation between the indentation stiffness of a micropolar solid and that of a conventional solid with the same Young modulus and Poisson ratio is derived.

**Keywords:** micropolar elasticity; spherical indentation; local point interpolation; boundary elements

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Pierson, G.; Taghite, M.; Bravetti, P.; Njiwa, R.K. Spherical Indentation of a Micropolar Solid: A Numerical Investigation Using the Local Point Interpolation-Boundary Element Method. *Appl. Mech. 2021*, *2*, 581-590. https://doi.org/10.3390/aplmech2030033

Received: 21 June 2021
Accepted: 17 August 2021
Published: 21 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Depth sensing indentation is commonly adopted for the determination of local elastic and plastic properties of small size samples. In the case of homogeneous and isotropic classical elastic materials deformed at small strain, the determination of the resolved elastic modulus follows a well established analytical solution of the normal compression of two smooth solids (see, e.g., Sneddon [1], K.L Johnson [2], for example). One of the useful results relates the applied load to the relative approach of the two solids or the material depression (if one of the solids is rigid). The load and the penetration are the only quantities recorded during the test. In the classical theory, the depression is also related to the radius of the projected contact area and the radius of the spherical punch. It has been shown that, in the small strain deformation regime, the useful relation between the material depression and the radius of the contact area remains valid in the plastic regime [3].

It is more and more evident that a finer description of the mechanical behavior of some material should take into account their microstructure. This can be done in the framework of generalized continuum theories which extend the conventional continuum mechanics for incorporating intrinsic microstructural effects in the mechanical behavior of materials. Amongst the various approaches, there is the so-called micromorphic medium of Eringen [4] which is widely accepted as the most successful phenomenological top-down approach. In this theory, the impact of the microstructure of the medium is expressed at the macroscopic scale through an incompatible microdeformation tensor. The latter can be specialized depending on the prominent microscopic effect at the macro scale. Accordingly, if the prominent microscopic effect is the individual rigid rotation of the material points, the microdeformation tensor is specialized to represent this micro-rotation. The corresponding
theory is that of micropolar materials also known as Cosserat materials [5]. The model was applied to porous and granular media (e.g., [6] (Ehlers, 1997) [7,8]) and even to living tissues such as bones (e.g., [9–13]). The identification of the material parameters of the micropolar elastic model has also been the subject of a number of works (see e.g., [14,15]. Depth sensing indentation is, among others, a valuable experimental tool for the mechanical characterization of some of this materials. The question arises as to whether the well-known relationships for homogeneous and isotropic materials remain valid in the case of elastic isotropic Cosserat materials.

The present numerical study intends to shed some light on this point. The numerical studies consider a diamond punch indenting a micropolar medium. The analysis uses the identified engineering parameters for micropolar elasticity, namely: E (the Young modulus), G (the shear modulus), v (the Poisson ratio), l_t (the characteristic length in torsion), l_b (the characteristic length in bending), χ (the polar ratio), and N (the coupling number).

The specifically developed numerical tool combines the advantages of the boundary elements method and a strong form point collocation method. The method called Local point interpolation-boundary element method has been initiated by Kouitat [16] in the context of anisotropic materials. It has since proved efficient for various fields including micropolar materials (see [17–19]).

The governing equations of micropolar media are recalled in Section 2 below. The main steps of the LPI-BEM and the global flowchart of the numerical tool are presented in Section 3. The numerical results are shown and discussed in Section 4.

## 2. Governing Equations

In the theory of micropolar medium occupying the domain Ω with boundary Γ, the material point x is attached to a triad of directors that can rotate and stretch. The material point possesses six degrees of freedom: the three components of the traditional displacement vector u and the three components of a microrotation vector φ. The field equations governing this type of medium when under quasi-static evolution without external body loads are ([20,21]):

$$ \sigma_{ji,j}(x) = 0 \qquad (1) $$

$$ m_{ji,j}(x) + \epsilon_{ijk} \sigma_{jk}(x) = 0 \qquad (2) $$

In these equations, σ represents the force stress tensor and m the couple stress tensor. Next, the case of a quasi-homogeneous and isotropic solid is considered, producing the following corresponding constitutive relations:

$$ \sigma_{ij}(x) = \lambda \epsilon_{rr} \delta_{ij} + 2\mu \epsilon_{ij} + \kappa u_{j,i} - \kappa \epsilon_{ijk} \varphi_k \qquad (3) $$

$$ m_{ij} = \alpha \omega_{rr} \delta_{ij} + (\beta + \gamma) \omega_{ij} + (\beta - \gamma) \bar{\omega}_{ij} \qquad (4) $$

$$ \epsilon_{ij} = (u_{i,j} + u_{j,i}) / 2 $$

$$ \omega_{ij} = (\varphi_{i,j} + \varphi_{j,i}) / 2 $$

$$ \bar{\omega}_{ij} = (\varphi_{i,j} - \varphi_{j,i}) / 2 $$

λ and μ are the Lamé constants; α, β, γ, and κ the micropolar constants. More precisely, α, β, and γ are rotation gradient moduli. α, (γ + β)/2 and (γ - β)/2 are also called Cosserat twist coefficients. κ is the rotation modulus is also known as the Cosserat couple modulus. The main symbols used in the paper are grouped together in the Appendix A.

With n_j as the outward normal vector on the boundary, the tractions acting at a regular point of the boundary are given by:

$$ t_i = \sigma_{ji} n_j \text{ and } m_i = m_{ji} n_j \qquad (5) $$
The material parameters must fulfill the following constraints:

$$ (3\lambda + 2\mu + \kappa) \geq 0; \ 2\mu + \kappa \geq 0; \ \kappa \geq 0; \ 3\alpha + \beta + \gamma \geq 0; \ \beta + \gamma \geq 0; \\ \gamma - \beta \geq 0 \qquad (6) $$

### 3. Solution Method

The boundary element method has already proven highly efficient for the solution of linear problems with well-established analytical fundamental solution. The boundary integral equation method in micropolar elasticity has been proposed by Sladek and Sladek [22]. An extension of the study to microstretch media is intended, thus this formulation will not be adopted. On application of the conventional BEM in this study, with Somigliana fundamental solution, the boundary element method (BEM) loses its principal appeal, namely the reduction of the problem dimension by one, due to traditional volume cells being needed in the so-called “field boundary element method” or “domain boundary element method”. This obstacle can be overcome by a number of strategies such as the dual reciprocity method [23] or radial integration method [24], which enable the conversion of volume integrals into surface ones.

In recent years, a large number of researchers have invested in the development of so-called meshless or meshfree methods. Among the various meshless approaches, the local point interpolation method is highly appealing because it is simple to implement. This approach falls in accuracy in the presence of Neumann boundary conditions, which are almost an inevitability when solving solid mechanic problems. Liu et al. have suggested a way to circumvent this difficulty by adopting the “weak-strong-form local point interpolation” method [25]. Recently, Kouitat [16] proposed a novel strategy that combines the best elements of both the conventional BEM and local point interpolation methods. This LPI-BEM approach has proved efficient in a number of contexts including micropolar elasticity [17,19]. This method is adopted in the present study and the main steps of the approach in the context of a micropolar elasticity are recalled below.

Let $v_u = \frac{\lambda - \kappa}{2(\lambda + \mu)}$, $\mu_u = \mu + \kappa$, then one could write the force stress tensor in the form:

$$ \sigma_{ij} = \frac{2 \mu_u v_u}{1 - 2v_u} \varepsilon_{rr} \delta_{ij} + 2\mu_u \varepsilon_{ij} - \kappa u_{i,j} - \kappa \varepsilon_{ijk} \varphi_k $$

Similarly, set $\mu_w = (\beta + \gamma)/2$ and $v_w = \frac{\alpha}{2\alpha + \beta + \gamma}$; then, the couple stress tensor reads:

$$ m_{ij} = \frac{2\mu_w v_w}{1 - 2v_w} \varphi_{r,r} \delta_{ij} + 2\mu_w \omega_{ij} + (\beta - \gamma)(\varphi_{i,j} - \varphi_{j,i})/2 $$

The calculations were based on the assumption that the kinematical primary variables are the sum of a complementary part and a particular term, namely: $u_i = u_i^H + u_i^P$ and $\varphi_i = \varphi_i^H + \varphi_i^P$.

The complementary fields satisfied the following homogeneous equations:

$$ \frac{\mu_u}{1 - 2v_u} u_{r,rj}^H + \mu_u u_{j,ii}^H = 0 $$

$$ \frac{\mu_w}{1 - 2v_w} \varphi_{r,rj}^H + \mu_w \varphi_{j,ii}^H = 0 $$

These equations which are of the Navier type were solved by the conventional boundary element method, thus producing the following systems of equations:

$$ [H_u] \{u^H\} = [G_u] \{t^H\} \text{ and } [H_\varphi] \{\varphi^H\} = [G_\varphi] \{m^H\} \qquad (7) $$

where $\{u^H\}$, $\{\varphi^H\}$ are vectors of nodal kinematical fields and $\{t^H\}$, $\{m^H\}$ are the vectors of nodal traction.
The particular fields solve:

$$ \frac{\mu_u}{1 - 2\nu_u} u_{r,ri}^P + \mu_u u_{i,jj}^P - \kappa u_{j,ij} - \kappa \epsilon_{jik} \varphi_{k,j} = 0 \quad (8) $$

$$ \frac{\mu_w}{1 - 2\nu_w} \varphi_{r,ri}^P + \mu_w \varphi_{i,jj}^P + \frac{(\beta - \gamma)(\varphi_{j,ji} - \varphi_{i,jj})}{2} - 2\kappa \varphi_i + \kappa \epsilon_{ijk} u_{k,j} = 0 \quad (9) $$

The tractions at a regular point on the boundary are written as:

$$ \begin{gathered} t_i = t_i^H + t_i^P + \delta t_i \quad \text{with} \quad t_i^A = \left( \frac{2 \mu_u v_u}{1 - 2\nu_u} \epsilon_{rr}^A \delta_{ij} + 2\mu_u \epsilon_{ij}^A \right) n_j \quad (A = H \text{ or } P) \\ \text{and} \quad \delta t_i = (\kappa u_{i,j} - \kappa \epsilon_{jik} \varphi_k) n_j \end{gathered} \quad (10) $$

$$ \begin{gathered} m_i = m_i^H + m_i^P + \delta m_i \quad \text{with} \quad t_i^A = \left( \frac{2\mu_w v_w}{1-2\nu_w} \omega_{rr}^A \delta_{ij} + 2\mu_w \omega_{ij}^A \right) n_j \quad (A = H \text{ or } P) \\ \text{and} \quad \delta m_i = (\beta - \gamma)(\varphi_{j,i} - \varphi_{i,j}) n_j / 2 \end{gathered} \quad (11) $$

Following this, the solution of Equations (7) and (8) were considered, using a local radial point collocation method. In this method [5], a field $\omega(x)$ is approximated as:

$$ \omega(x) = \sum_{i=1}^{N} R_i(r)a_i + \sum_{j=1}^{M} p_j(x)b_j $$

with the following constraints: $\sum_{i=1}^{N} p_j(x) a_i = 0, j = 1 - M$ and $i = 1 - N$.

Here, $R_i(r)$ is the selected radial basis functions, $N$ the number of nodes in the neighborhood (support domain) of point $x$, and $M$ the number of monomial terms in the selected polynomial basis $P_j(x)$.

Coefficients $a_i$ and $b_j$ can be determined by enforcing the approximation to be satisfied at the $N$ centers. Following some algebraic manipulations, coefficients $a_i$ and $b_j$ are expressed in terms of the field nodal values, and the interpolation is written in the following compact form:

$$ \omega^h(x) = [\Phi(x)]\{\omega_{/L}\} \quad (12) $$

When adopting interpolation (12) for all kinematical fields, at a given collocation center, the following were obtained:

$$ \begin{gathered} [B(\nabla)]^T [C_u] [B(\nabla)] [\tilde{\Phi}] \{u_{/L}^P\} + \kappa [\tilde{\Phi}_1] \{u_{/L}\} + \kappa [\tilde{\Phi}_2] \{\varphi_{/L}\} = \{0\} \\ [B(\nabla)]^T [C_\varphi] [B(\nabla)] [\tilde{\Phi}] \{\varphi_{/L}^P\} + \kappa [\tilde{\Phi}_0] \{\varphi_{/L}\} + (\beta - \gamma) [\tilde{\Phi}_1] \{\varphi_{/L}\} + \kappa [\tilde{\Phi}_3] \{u_{/L}\} = \{0\} \end{gathered} $$

In the above,

$$ \{\nabla\} = (\partial/\partial x \quad \partial/\partial y \quad \partial/\partial z)^T, \quad \{z\} = (z_1 \quad z_2 \quad z_3)^T, \quad [B(z)] = \begin{bmatrix} z_1 & 0 & 0 & z_2 & z_3 & 0 \\ 0 & z_2 & 0 & z_1 & 0 & z_3 \\ 0 & 0 & z_3 & 0 & z_1 & z_2 \end{bmatrix}^T $$

Matrix $B$ was given in terms of vector $z = (z_1, z_2, z_3)^T$ and matrix $C$ was the Voigt representation of the elasticity tensor.

On collection of the above equations for all the internal collocation centers, taking the assumption that the particular integrals are identically zero at all boundary points, the following forms of systems of equations were obtained:

$$ \{u^P\} = [B_{uu}]\{u\} + [B_{u\varphi}]\{\varphi\} \quad (13) $$

$$ \{\varphi^P\} = [A_{\varphi\varphi}]\{\varphi\} + [B_{\varphi u}]\{u\} \quad (14) $$
Following a similar strategy, the tractions at the boundary points could be written in the following forms:

$$ \begin{aligned} \{t\} &= \{t^H\} + [AK_{uu}]\{u\} + [AK_{u\varphi}]\{\varphi\} \\ \{m\} &= \{m^H\} + [AK_{\varphi u}]\{u\} + [AK_{\varphi\varphi}]\{\varphi\} \end{aligned} $$

After conducting some algebraic manipulations, the final coupled systems of equations were of the following forms:

$$ [\overline{H}_u] \{u\} + [H_{u\varphi}] \{\varphi\} = [G_u] \{t\} \quad (15) $$

$$ [\overline{H}_{\varphi}] \{\varphi\} + [H_{\varphi u}] \{u\} = [G_{\varphi}] \{m\} \quad (16) $$

It is particularly worthy of mention that the final equations contained similar boundary primary variables and internal kinematic unknowns to those of a traditional BEM. Boundary conditions can be taken into account as in standard practice.

When the problem at hand is concerned with more than one medium, then the above equations are valid for each sub-domain. Let $\Gamma^I$ denote a perfectly bonded interface between body A and body B. Then, for point $X^A$ of body A and point $X^B$ of body B sharing the same geometrical location in $\Gamma^I$, the following conditions must be fulfilled:

$$ t_i^A + t_i^B = 0 (i = 1, 2, 3), \quad u_i^A - u_i^B = 0 (i = 1, 2, 3) \quad (17) $$

$$ m_i^A + m_i^B = 0 (i = 1, 2, 3), \quad \varphi_i^A - \varphi_i^B = 0 (i = 1, 2, 3) \quad (18) $$

In the case of a non-conforming contact problem, the equations must be supplemented by the unilateral contact conditions between the contacting bodies. Focusing on the case of an elastic punch, these conditions are written only for the macro displacement and macro traction. Indeed, it is believed that it is not possible to apply micro torque at the specimen boundary. The micro torque is assumed to take zero value on all boundaries. On the boundary of the contacting body A (resp. B), a potential contact area $\Gamma^A$ (resp. $\Gamma^B$) is defined a priori. For a node $X^A$ of $\Gamma^A$ and a target node $X^B$ on $\Gamma^B$, the following relations apply:

$$ t_i^A + t_i^B = 0 (i = 1, 2, 3) \quad (19) $$

$$ d_n = u_n^B - u_n^A \le g_0, \quad t_n^B \le 0 \text{ and } t_n^B = 0 \quad (20) $$

In the above relations, $g_0$ denotes the initial normal gap between $X^A$ of $\Gamma^A$ and $X^B$ on $\Gamma^B$. $u_n^A$ and $u_n^B$ stand for the normal displacements of points A and B in the direction of the common normal, which is usually taken as that of the flat specimen. $t_n^B$ is the normal traction at point B and $t_n^A$ is the tangential traction. It assumed here that the contact is frictionless.

Indentation by a spherical punch is a non-conforming contact problem in the sense that, for a given load, the contact surface is not known in advance. It is then advised to solve the above coupled systems of equations incrementally. Let $\Delta f$ be the increment of the field $f$ from the previously converged solution to the actually sought solution. Then, the incremental forms of systems (15) and (16) are:

$$ [\overline{H}_u] \{\Delta u\} + [H_{u\varphi}] \{\Delta \varphi\} = [G_u] \{\Delta t\} \quad (21) $$

$$ [\overline{H}_{\varphi}] \{\Delta \varphi\} + [H_{\varphi u}] \{\Delta u\} = [G_{\varphi}] \{\Delta m\} \quad (22) $$

Applying the boundary conditions as usual, accounting eventually for Equation 16, global systems of the following forms are obtained:

$$ [A_{uu}]\{Y_u\} + [A_{u\varphi}]\{\Delta\varphi\} = \{F_u\} \quad (23) $$

$$ [A_{\varphi\varphi}]\{Y_{\varphi}\} + [A_{\varphi u}]\{\Delta u\} = \{F_{\varphi}\} \quad (24) $$
In the case of a non-conforming contact problem, Equation (23) must be supplemented with relations (19) and (20). The resulting system of equations is nonlinear and non-differentiable. The algorithm introduced by Christensen [26] has been adopted. The pseudo code implemented for the solution of the problem is as follows:

For the actual load increment:

1. Set $k = 0, \{Y_u^k\} = 0, \{\Delta\phi^k\} = 0$;

2. Solve $[A_{uu}]\{Y_u^{k+1}\} + [A_{u\phi}]\{\Delta\phi^k\} = \{F_u\}$;

3. Solve $[A_{\phi\phi}]\{Y_{\phi}^{k+1}\} + [A_{\phi u}]\{\Delta u^{k+1}\} = \{F_{\phi}\}$;

4. If the $Y_u^{k+1} - Y_u^k / Y_u^{k+1} > \varepsilon$ (a given tolerance) go to 2.

## 4. Numerical Results

In this work, the multi-quadrics radial basis functions are applied as follows: $R_i(r) = (r_i^2 + c^2)^q$, where $r_i = x - x_i$ and $c$ and $q$ were known as shape parameters. Shape parameter $c$ was taken proportionally to minimum distance $d_0$, defined as the maximum value among the minimum distances in the $x, y$, and $z$ directions between collocation points.

Let us focus on the elastic indentation problem. For conventional materials, a famous relation of the Hertzian theory of contact states that: $\delta = \frac{a^2}{R}$, with ($\delta$) the indenter depression, ($a$) the radius of the projected contact disc, and ($R$) the radius of the spherical punch. Another famous relationship of the Hertzian theory of contact is the relation between the applied load ($P$), the punch radius, and the contact radius, $P = K_0 \frac{a^3}{R}$ with $K_0 = \frac{4}{3} \frac{E}{1-\nu^2}$. Knowing the value of $K_0$ and postulating the value of the Poisson ratio, the Young modulus can be calculated.

A part of a diamond ball with radius 100 µm is considered. The boundary of the ball piece was subdivided into 100 elements with 402 boundary nodes including 153 nodes of the potential contact area. The considered micropolar solid is a cuboid with depth 5 mm and squared cross section with size 600 µm. Its boundary was subdivided into 236 elements with 1022 boundary nodes. These latter were supplemented with 1069 internal collocation centers.

The material constants of micropolar elasticity (see Equations (3) and (4)) are related to technical material parameters that have been identified individually via experiments [6]. These material parameters are:

The Young modulus $E = \frac{(2\mu+\kappa)(3\lambda+2\mu+\kappa)}{2\lambda+2\mu+\kappa}$;

The shear modulus $G = \mu + \kappa/2$;

The Poisson ratio $\nu = \frac{\lambda}{2\lambda+2\mu+\kappa}$;

The characteristic length in torsion $l_t = \sqrt{\frac{\beta+\gamma}{2\mu+\kappa}}$;

The characteristic length in bending $l_b = \sqrt{\frac{\gamma}{2(2\mu+\kappa)}}$;

The polar ratio $\chi = \frac{\beta+\gamma}{\alpha+\beta+\gamma}$;

The coupling number $N = \sqrt{\frac{\kappa}{2(\mu+\kappa)}}$.

Let us remind the reader that the values of $E, G$, and $\nu$ are those that can be assessed by a simple tension test.

### 4.1. Influence of the Coupling Number

In a first set of simulation, five sets of material constants were selected (see Table 1 below) leading to the following values of material parameters: $E = 25,789.47$ MPa, $G = 8750$ MPa, $\nu = 0.474$, $l_t = 0.1095$ mm, $l_b = 0.0775$ mm, and $\chi = 1$. It should be noted that, $E, G$ and $\nu$ have the same values for all considered cases.
<table><caption>**Table 1.** Micropolar material constants adopted for the study of the influence of the coupling number.</caption><thead><tr><th>λ (MPa)</th><th>μ (MPa)</th><th>κ (MPa)</th><th>α (N)</th><th>β (N)</th><th>γ (N)</th><th>N<sup>2</sup></th></tr></thead><tbody><tr><td>157,500</td><td>8166.67</td><td>1166.67</td><td>0</td><td>0</td><td>210</td><td>1/16</td></tr><tr><td>157,500</td><td>7000</td><td>3500</td><td>0</td><td>0</td><td>210</td><td>1/6</td></tr><tr><td>157,500</td><td>5833.33</td><td>5833.33</td><td>0</td><td>0</td><td>210</td><td>1/4</td></tr><tr><td>157,500</td><td>3141.026</td><td>11,217.95</td><td>0</td><td>0</td><td>210</td><td>25/64</td></tr><tr><td>157,500</td><td>-2500</td><td>22,500</td><td>0</td><td>0</td><td>210</td><td>9/16</td></tr></tbody></table>

As can be observed, the square of the coupling number varies from 1/16 to 9/16.

A load of 50 mN was applied on the punch in five load steps. The indenter depression ($\delta$), the contact pressure distribution ($p_n$), and the interval containing the contact radius ($a$) were calculated. Let $R$ denote the radius of the spherical punch. It is found that the famous Hertzian relation $\delta = \frac{a^2}{R}$ is still valid for the tested cases of micropolar solids. It is then possible to say that, having measured the punch depression and knowing its radius of curvature, the contact radius can be reached.

In Figure 1 below, the plots of $P$ versus $\frac{a^3}{R}$ are shown for the considered material parameters. It can be observed that all the curves are linear. A simple linear fit then allowed to determine the parameters of the fitting lines. The results are collected in Table 2 below.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Plots of load ($P$) versus $(a^3/R)$ for different values of the coupling number $N^2$.</figcaption></figure>

<table><caption>**Table 2.** Parameters of the line fit of the load ($P$) versus $(a^3/R)$.</caption><thead><tr><th>$N^2$</th><th>Equation</th></tr></thead><tbody><tr><td>1/16</td><td>$P = 46,374 \frac{a^3}{R} + 10^{-4}$</td></tr><tr><td>1/6</td><td>$P = 52,286 \frac{a^3}{R} + 8 \ 10^{-6}$</td></tr><tr><td>1/4</td><td>$P = 57,780 \frac{a^3}{R} + 10^{-4}$</td></tr><tr><td>25/64</td><td>$P = 71,325 \frac{a^3}{R} + 6 \ 10^{-5}$</td></tr><tr><td>9/16</td><td>$P = 99,582 \frac{a^3}{R} + 2 \ 10^{-4}$</td></tr></tbody></table>

It should be noted that, according to Hertzian theory of contact, for the given values of Young modulus and Poisson ratio, the slope of the fitting line should be $K_0 = 44,333.33$ MPa. It is clear from data in Table 2 that a relation of the form $P = K_N \frac{a^3}{R}$ is still valid (index $N$ is
introduced to represent the coupling number). However, the value of the parameter $K_N$ varies with the coupling number.

The ratios of the calculated $K_N$ over $K_0$ are collected in Table 3 below.

<table><caption>**Table 3.** Comparison of the ratios $K_N/K_0$ and $1/(1-N^2)$ for different values of $N$.</caption><thead><tr><th>N<sup>2</sup></th><th>1/16</th><th>1/6</th><th>1/4</th><th>25/64</th><th>9/16</th></tr></thead><tbody><tr><td>K<sub>N</sub>/K<sub>0</sub></td><td>1.046</td><td>1.179</td><td>1.303</td><td>1.609</td><td>2.246</td></tr><tr><td>1/(1-N<sup>2</sup>)</td><td>1.07</td><td>1.2</td><td>1.33</td><td>1.64</td><td>2.285</td></tr></tbody></table>

The results in this table show that: $K_N \approx \frac{K_0}{1-N^2}$.

## 4.2. Influence of the Characteristic Length in Torsion and the Polar Ratio

Other sets of material constants were adopted (see Table 4). They led to the same coupling numbers as in the former case. Now, $l_t = 0.1342$ and $\chi = 1.5$.

<table><caption>**Table 4.** Micropolar material constants adopted for the case $l_t = 0.1342$ and $\chi = 1.5$.</caption><thead><tr><th>λ (MPa)</th><th>μ (MPa)</th><th>κ (MPa)</th><th>α (N)</th><th>β (N)</th><th>γ (N)</th><th>N<sup>2</sup></th></tr></thead><tbody><tr><td>157,500</td><td>8166.67</td><td>1166.67</td><td>-105</td><td>105</td><td>210</td><td>1/16</td></tr><tr><td>157,500</td><td>7000</td><td>3500</td><td>-105</td><td>105</td><td>210</td><td>1/6</td></tr><tr><td>157,500</td><td>5833.33</td><td>5833.33</td><td>-105</td><td>105</td><td>210</td><td>1/4</td></tr><tr><td>157,500</td><td>3141.026</td><td>11,217.95</td><td>-105</td><td>105</td><td>210</td><td>25/64</td></tr><tr><td>157,500</td><td>-2500</td><td>22,500</td><td>-105</td><td>105</td><td>210</td><td>9/16</td></tr></tbody></table>

The values obtained for $K_N$ were those of Table 2. This means that the results are not disturbed by parameters $l_t$ and $\chi$.

## 4.3. Influence of the Characteristic Length in Bending

Let us now consider the potential effect of the characteristic length in bending. For this purpose, the material parameters in Table 5 are adopted.

<table><caption>**Table 5.** Micropolar material constants adopted for the case $l_b = 0.1095$ and $\chi = 0.75$.</caption><thead><tr><th>λ (MPa)</th><th>μ (MPa)</th><th>κ (MPa)</th><th>α (N)</th><th>β (N)</th><th>γ (N)</th><th>N<sup>2</sup></th></tr></thead><tbody><tr><td>157,500</td><td>8166.67</td><td>1166.67</td><td>105</td><td>-105</td><td>420</td><td>1/16</td></tr><tr><td>157,500</td><td>7000</td><td>3500</td><td>105</td><td>-105</td><td>420</td><td>1/6</td></tr><tr><td>157,500</td><td>5833.33</td><td>5833.33</td><td>105</td><td>-105</td><td>420</td><td>1/4</td></tr><tr><td>157,500</td><td>3141.026</td><td>11,217.95</td><td>105</td><td>-105</td><td>420</td><td>25/64</td></tr><tr><td>157,500</td><td>-2500</td><td>22,500</td><td>105</td><td>-105</td><td>420</td><td>9/16</td></tr></tbody></table>

The characteristic length in torsion is that of the preceding case. Now, $l_b = 0.1095$ and $\chi = 0.75$.

Once more, the results of the calculation show that values of $K_N$ are still those in the table.

# 5. Conclusions

This paper concerns the numerical simulation of elastic indentation of a sphere into a flat micropolar elastic specimen. First, details of specifically developed numerical tool are given. Then, it is shown that the relation $\delta = \frac{a^2}{R}$, extremely useful for direct computation of the real contact area as a function on the depression $\delta$ is equally valid for a micropolar elastic medium. It was then possible to plot the load ($P$) as a function of $\frac{a^3}{R}$. It then appears
that a relation of the Hertzian type $P = K_N \frac{a^3}{R}$ remains valid. The coefficient $K_N$ increases with the value of the coupling number $N$ and is unaltered by other micropolar parameters, namely the characteristic lengths and the polar ratio. It is also found that $K_N \approx \frac{K_0}{1-N^2}$ where $K_0$ is the corresponding value as given by the Hertzian theory of contact. In conclusion, the only technical parameter affecting the indentation response of a micropolar medium is the coupling number.

Suppose that the values of parameters $E$, $G$, and $\nu$ have been obtained by a simple tension test. Then, the value of $K_0$ can be calculated. Knowing $K_0$, one can determine the coupling number $N$ and hence the Cosserat couple modulus $\kappa$.

Further studies are needed to allow identification of the other parameters of the model.

**Author Contributions:** Resources, M.T. and R.K.N.; Software, G.P., R.K.N.; funding acquisition and project supervision, P.B. and R.K.N.; writing—original draft preparation, G.P.; writing—review and editing, G.P. and R.K.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partly funded by the IMPLANT operation (Positive Medical Personalized innovants Implantables) co-financed by the European Union through ERDF-ESF Operational "Programme Lorraine et Massif des Vosges 2014–2020".

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no part in the design of the study; in the production, evaluation, or interpretation of results; in the drafting of the manuscript or the final decision to publish the data.

## Appendix A

List of notations

<table><tr><td>σ<sub>ji</sub></td><td>Non symmetric macro stress tensor</td></tr><tr><td>m<sub>ji</sub></td><td>Couple stress tensor</td></tr><tr><td>u<sub>i</sub></td><td>Macro displacement vector</td></tr><tr><td>φ<sub>i</sub></td><td>Micro rotation vector</td></tr><tr><td>u<sub>i,j</sub></td><td>Displacement gradient tensor</td></tr><tr><td>ε<sub>ij</sub></td><td>Small strain tensor</td></tr><tr><td>φ<sub>i,j</sub></td><td>Micro rotation gradient tensor</td></tr><tr><td>ω<sub>ij</sub></td><td>Symmetric part of φ<sub>i,j</sub></td></tr><tr><td>ω&#773;<sub>ij</sub></td><td>Anti-symmetric part of φ<sub>i,j</sub></td></tr><tr><td>t<sub>i</sub></td><td>Traction vector</td></tr><tr><td>m<sub>i</sub></td><td>Micro torque vector</td></tr><tr><td>λ</td><td>Lamé coefficient</td></tr><tr><td>μ</td><td>Lamé coefficient</td></tr><tr><td>κ</td><td>Cosserat couple modulus</td></tr><tr><td>α, β, γ</td><td>Rotation gradient moduli</td></tr><tr><td>δ</td><td>Material depression</td></tr><tr><td>a</td><td>Radius of the projected contact area</td></tr><tr><td>R</td><td>Radius of the indenter</td></tr><tr><td>P</td><td>Applied load</td></tr><tr><td>E</td><td>Young modulus</td></tr><tr><td>G</td><td>Shear modulus</td></tr><tr><td>ν</td><td>Poisson ratio</td></tr><tr><td>l<sub>t</sub></td><td>Characteristic length in torsion</td></tr><tr><td>l<sub>b</sub></td><td>Characteristic length in bending</td></tr><tr><td>χ</td><td>Polar ratio</td></tr><tr><td>N</td><td>Coupling number</td></tr><tr><td>K<sub>N</sub></td><td>Contact stiffness</td></tr></table>

## References

1. Sneddon, I.N. *Fourier Transforms*; McGraw-Hill: New York, NY, USA, 1951.

2. Johnson, K.L. *Contact mechanics*, 9th ed.; Cambridge University Press: Cambridge, UK, 2003.

3. Oumarou, N.; Jehl, J.P.; Kouitat, R.; Stempfle, P. On the variation of mechanical parameters obtained from spherical depth sensing indentation. *Int. J. Surf. Sci. Eng.* **2010**, *4*, 416. [CrossRef]
4. Eringen, A.C.; Suhubi, E.S. Nonlinear theory of simple micro-elastic solids—I. *Int. J. Eng. Sci.* **1964**, *2*, 189–203. [CrossRef]

5. Cosserat, E.; Cosserat, F. *Théorie Des Corps Déformables*; Herman et Fils: Paris, France, 1909.

6. Lakes, R.S. Experimental microelasticity of two porous solids. *Int. J. Solids Struct.* **1986**, *22*, 55–63. [CrossRef]

7. Diebels, S. Micropolar mixture models on the basis of the Theory of Porous Media. In *Porous Media*; Ehlers, P.D.-I.W., Bluhm, P.D.D.-I.J., Eds.; Springer: Berlin/Heidelberg, Germany, 2002; pp. 121–145. [CrossRef]

8. Goddard, J.D. From Granular Matter to Generalized Continuum. In *Mathematical Models of Granular Matter*; Capriz, G., Mariano, P.M., Giovine, P., Eds.; Springer: Berlin/Heidelberg, Germany, 2008; pp. 1–22. [CrossRef]

9. Park, H.C.; Lakes, R.S. Cosserat micromechanics of human bone: Strain redistribution by a hydration sensitive constituent. *J. Biomech.* **1986**, *19*, 385–397. [CrossRef]

10. Cowin, S.C. Bone poroelasticity. *J. Biomech.* **1999**, *32*, 217–238. [CrossRef]

11. Rosenberg, J.; Cimrman, R. Microcontinuum approach in biomechanical modelling. *Math. Comput. Simul.* **2003**, *61*, 249–260. [CrossRef]

12. Ramtani, S. Electro-mechanics of bone remodelling. *Int. J. Eng. Sci.* **2008**, *46*, 1173–1182. [CrossRef]

13. Eringen, A.C. Electromagnetic theory of microstretch elasticity and bone modelling. *Int. J. Eng. Sci.* **2004**, *42*, 231–242. [CrossRef]

14. Diebels, S.; Geringer, A. Micromechanical and macromechanical modelling of foams: Identification of Cosserat parameters: Micromechanical and macromechanical modelling of foams. *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik* **2014**, *94*, 414–420. [CrossRef]

15. Neff, P.; Jeong, J.; Fischle, A. Stable identification of linear isotropic Cosserat parameters: Bounded stiffness in bending and torsion implies conformal invariance of curvature. *Acta Mech.* **2010**, *211*, 237–249. [CrossRef]

16. Njiwa, R.K. Isotropic-BEM coupled with a local point interpolation method for the solution of 3D-anisotropic elasticity problems. *Eng. Anal. Bound. Elem.* **2011**, *35*, 611–615. [CrossRef]

17. Pierson, G.; Bravetti, P.; Njiwa, R.K. Interaction implant-bone as a micropolar elastic medium: Porosity impact of the hard living media. *Int. J. Theor. Appl. Multiscale Mech.* **2020**, *3*, 229. [CrossRef]

18. Njiwa, R.K.; Pierson, G.; Voignier, A. Coupling BEM and the Local Point Interpolation for the Solution of Anisotropic Elastic Nonlinear, Multi-Physics and Multi-Fields Problems. *Int. J. Comput. Methods* **2019**, *17*, 1950067. [CrossRef]

19. Pierson, G.; Kouitat-Njiwa, R.; Bravetti, P. A boundary elements only solution method for 3D micropolar elasticity. *Eng. Anal. Bound. Elem.* **2021**, *123*, 84–92. [CrossRef]

20. Eringen, A.C. *Microcontinuum Field Theories*; Springer: New York, NY, USA, 1999. [CrossRef]

21. Iesan, D.; Pompei, A. On the equilibrium theory of microstretch elastic solids. *Int. J. Eng. Sci.* **1995**, *33*, 399–410. [CrossRef]

22. Sládek, V.; Sládek, J. Boundary integral equation method in micropolar elasticity. *Appl. Math. Model.* **1983**, *7*, 433–440. [CrossRef]

23. Nardini, D.; Brebbia, C.A. A new approach to free vibration analysis using boundary elements. *Appl. Math. Model.* **1983**, *7*, 157–162. [CrossRef]

24. Gao, X.-W. The radial integration method for evaluation of domain integrals with boundary-only discretization. *Eng. Anal. Bound. Elem.* **2002**, *26*, 905–916. [CrossRef]

25. Liu, G.R.; Gu, Y.T. A local radial point interpolation method (LRPIM) for free vibration analyses of 2-D solids. *J. Sound Vib.* **2001**, *246*, 29–46. [CrossRef]

26. Christensen, P.W. A semi-smooth Newton method for elasto-plastic contact problems. *Int. J. Solids Struct.* **2002**, *39*, 2323–2341. [CrossRef]