# On Solving Nonlinear Elasticity Problems Using a Boundary-Elements-Based Solution Method

Aly Rachid Korbeogo, Bernard Kaka Bonzi, Richard Kouitat Njiwa

► To cite this version:

Aly Rachid Korbeogo, Bernard Kaka Bonzi, Richard Kouitat Njiwa. On Solving Nonlinear Elasticity Problems Using a Boundary-Elements-Based Solution Method. Applied Mechanics, 2023, 4 (4), pp.1240 - 1259. 10.3390/aplmech4040064. hal-05041391

HAL Id: hal-05041391

https://hal.univ-lorraine.fr/hal-05041391v1

Submitted on 21 Apr 2025

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

# On Solving Nonlinear Elasticity Problems Using a Boundary-Elements-Based Solution Method

Aly Rachid Korbeogo<fnref n="1" />, Bernard Kaka Bonzi<fnref n="2" /> and Richard Kouitat Njiwa<fnref n="3" />,*

**Abstract:** The attractiveness of the boundary element method—the reduction in the problem dimension by one—is lost when solving nonlinear solid mechanics problems. The point collocation method applied to strong-form differential equations is appealing because it is easy to implement. The method becomes inaccurate in the presence of traction boundary conditions, which are inevitable in solid mechanics. A judicious combination of the point collocation and the boundary integral formulation of Navier’s equation allows a pure boundary element method to be obtained for the solution of nonlinear elasticity problems. The potential of the approach is investigated in some simple examples considering isotropic and anisotropic material models in the total Lagrangian framework.

**Keywords:** elasticity; anisotropy; large deformation; BEM; meshfree strong form; radial basis function

## 1. Introduction

The boundary element method (BEM) is entrenched as a powerful numerical method for the solution of linear problems with well-established fundamental solutions including structural problems (see, e.g., [1–3]). Its efficiency surpasses that of the well-known and versatile finite element method (FEM) in the case of linear fracture mechanics problems. The consideration of a structure undergoing large strain elastic and/or plastic deformations is relevant in many engineering and scientific fields including the study of soft biological tissues. In the case of geometrical and/or material nonlinearities, the integral formulation of the field equations involves a domain integral term. The method is thus known as the field boundary element method (FBEM) or the domain boundary element method (e.g., [4–6]). The solution’s procedure uses volume cells in the same manner as in the FEM. Thus, in this common approach, the main appeal of the BEM (reduction in the domain dimension by one) is tarnished to some extent. A number of strategies have been developed in order to convert domain integrals into boundary integrals. One can mention the dual reciprocity method [7], the radial integration method [8], the generalized boundary element method [9], and the multiple reciprocity method which has been recently applied in the context of large deformation [10]. Other boundary-element-based approaches have been proven effective for a variety of problems: the local integral equation [11] and the analogue equation method [12]. A boundary-element-only formulation for the large-strain incremental elasticity was presented by Brun et al. [13,14]. The approach requires the use of a special anisotropic Green fundamental solution, the implementation of which is tedious. Recent years have seen an increasing interest in the development of so-called meshless or meshfree methods. They can be classified into two main categories: meshless methods based on the weak form of the partial differential equations and those based on the strong form [15]. The first category was successfully applied to structural problems at large

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

### check for updates

**Citation:** Korbeogo, A.R.; Bonzi, B.K.; Kouitat Njiwa, R. On Solving Nonlinear Elasticity Problems Using a Boundary-Elements-Based Solution Method. *Appl. Mech. 2023*, *4*, 1240–1259. https://doi.org/10.3390/aplmech4040064

Received: 21 October 2023

Revised: 26 November 2023

Accepted: 6 December 2023

Published: 12 December 2023

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

**Copyright:** © 2023 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/).

<footnote n="1">¹ Laboratoire des Matériaux et Environnement (LAME), Ecole Polytechnique de Ouagadougou, Ouaga 03 B.P. 7021, Burkina Faso; alykorbeogo@yahoo.fr</footnote>

<footnote n="2">² Laboratory of Mathematics and Computer Science (LAMI), Department of Mathematics and Computer Science, Université Joseph Ki-Zerbo & INSEM, Ouaga 03 B.P. 7021, Burkina Faso; bernardkbonzi@gmail.com</footnote>

<footnote n="3">³ Institut Jean Lamour, CNRS UMR 7198, Université de Lorraine, Campus Artem, 2 allée André Guinier, 54000 Nancy, France</footnote>

<footnote>* Correspondence: richard.kouitat@univ-lorraine.fr</footnote>
strains with material nonlinearities (e.g., [16–26]). The point collocation method belongs to the class of meshfree methods based on the strong form of the differential equations. It has been successful for small-strain problems [27–30]. The method is accurate as long as the boundary conditions are of the Dirichlet type. In the presence of Neuman-type boundary conditions, which is practically inevitable in solid mechanics, the solution’s accuracy deteriorates. One way to overcome this difficulty is to use the weak form at boundary points with Neumann-type boundary conditions. This method is known as the meshfree weak-strong form and was proposed by Liu and Gu [28]. To our knowledge, this attractive point collocation method has not yet been applied to the solution of problems involving large strains.

The purpose of the present study is to develop a BEM-based solution for 3D nonlinear elasticity problems. The proposed strategy judiciously combines the strong-form point collocation method with the conventional BEM for the small-strain isotropic elasticity. The method which is herein called the LPI-BEM (local point interpolation-boundary element method) has been successful for the solution of small-strain anisotropic elasticity, multi-physics, and multi-field problems (e.g., [31–34]). The method uses an additive partition of the displacement field into complementary and particular parts. The complementary field satisfies a Navier-type equation. The particular integral is obtained by a point collocation solution of strong-form differential equations. The final system of nonlinear equations to be solved is similar in form to that obtained by other methods such as the field boundary element method.

This paper is organized as follows. In Section 2, the main steps of the LPI-BEM for non-linear elastic problems are presented. The section starts with a reminder of the field equations followed by the main steps of the BEM. The section continues with the description of the point collocation method until the final nonlinear system of equations is obtained. Section 3 is devoted to numerical examples. The section begins with a presentation of the considered solids: a Saint Venant–Kirchhoff solid, a Neo–Hookean solid, a Mooney–Rivlin solid, a transversely isotropic Saint Venant–Kirchhoff, and a transversely isotropic Neo–Hookean solid. The validity of this easy-to-implement approach is investigated by considering uni-axial loading states for which analytical solutions are available. Results of the application of the solution to the problems of a constrained tension of a cubic specimen and of the bending of a bar are also presented.

## 2. Governing Equations and Solution Method

### 2.1. Basic Definitions and Governing Equations

Let us consider a nonlinear elastic solid undergoing large elastic deformation. We adopt the total Lagrangian formulation and express the equilibrium equations in the undeformed (stress-free) initial configuration of the body. The homogeneous isotropic solid occupies the initial domain $\Omega_0$ with boundary $\Gamma_0$. The position vector of a particle is denoted by $\underline{X}$. After deformation, the body occupies the domain $\Omega$ and the position vector of a particle is denoted by $\underline{x}$. It is assumed that there is a twice differentiable invertible map $\chi$ that takes $\Omega_0$ to $\Omega$, that is, $\underline{x} = \chi(\underline{X})$. The displacement vector $\underline{u}$ is given by

$$ \underline{u}(\underline{X}) = \underline{x} - \underline{X} = \chi(\underline{X}) - \underline{X}. $$

A measure of the deformation is described by the deformation gradient $\underline{F}$ relative to $\underline{X}$ given by

$$ \underline{F} = \frac{\partial\chi(\underline{X})}{\partial\underline{X}} = \frac{\partial\underline{u}(\underline{X})}{\partial\underline{X}} + \underline{1}. $$

The identity tensor is denoted as $\underline{1}$. The right Cauchy–Green tensor ($\underline{G}$) can be found from $\underline{F}$ according to $\underline{G} = \underline{F}^T \underline{F}$.

Subsequently, a Cartesian system of coordinates and an indicial notation with the associated summation convention over repeated indices are adopted. When body forces
are neglected, the equilibrium equations with respect to the reference configuration are (see, e.g., [35,36])

$$ \frac{\partial P_{kj}(\underline{X})}{\partial X_j} = 0 \qquad (1) $$

The tensor $\underline{P}$ is the first Piola–Kirchhoff stress tensor which is non-symmetric. Equation (1) must be supplemented by properly defined boundary conditions in terms of known displacement and traction. The stress tensor $\underline{P}$ is expressed in terms of the symmetric second Piola–Kirchhoff stress tensor $\underline{S}$ by $P_{kj} = F_{ki}S_{ij}$. In the case of a hyper-elastic law, $\underline{S}$ is defined from an elastic potential $W$ by $\underline{S} = 2 \frac{\partial W}{\partial \underline{G}}$. Note that Equation (1) is a highly nonlinear partial differential equation.

## 2.2. Solution Method Using the LPI-BEM

In this section, we present the main lines of the method as applied for the solution of elastic large deformation problems.

First, an isotropic stress tensor is added and subtracted to Equation (1) as follows:

$$ \frac{\partial(C_{kjm} u_{m,n} + (P_{kj} - C_{kjm} u_{m,n}))}{\partial X_j} = 0. $$

In the above expression, $C_{kjm}$ is a fourth-order tensor of the elastic constants of an isotropic homogeneous material. Now, let us assume that the displacement $u$ can be viewed as the sum of two terms: a complementary term ($u^c$) which behaves as the displacement vector in small-strain deformation of an isotropic elastic solid, and a particular term ($u^p$) which accounts for all nonlinear effects of the deformation process. Bearing this in mind, Equation (1) is rewritten as

$$ \frac{\partial(\sigma_{kj}^c + \sigma_{kj}^p + T_{kj})}{\partial X_j} = 0 \qquad (2) $$

In Equation (2), we introduce the following:

$$ \sigma_{jk}^c = C_{jkmn} u_{m,n}^c, \quad \sigma_{jk}^p = C_{jkmn} u_{m,n}^p \quad \text{and} \quad T_{jk} = P_{jk} - C_{jkmn} u_{m,n} $$

where $\sigma_{jk}^c$ is the complementary stress tensor, $\sigma_{jk}^p$ is the particular stress tensor, and $T_{jk}$ is the stress difference tensor.

Let us remind the reader that the adopted partition of the kinematical field is not new in the BEM community. Indeed, the analog equation method uses a similar partition with a complementary solution which satisfies the Laplace equation [12]. It is worth mentioning that the method of particular integrals is now widely used in many problems [37].

In our approach, the complementary solution is defined as the solution of a Navier-type differential equation:

$$ \frac{\partial(C_{kjm} u_{m,n}^c)}{\partial X_j} = 0 \qquad (3) $$

The associated boundary traction is $t_k^c = C_{kjm} u_{m,n}^c N_j$ with $N_j$ being the normal vector to the boundary in the undeformed configuration. The relation between this complementary traction vector and the true boundary traction is made clear below.

The particular integral is obtained by solving Equation (4):

$$ \frac{\partial(C_{kjm} u_{m,n}^p + T_{kj})}{\partial X_j} = \frac{\partial(\sigma_{kj}^p + T_{kj})}{\partial X_j} = 0 \qquad (4) $$
Complementary field:

Equation (3) is similar to that of the classical small-strain isotropic elastostatics, and its well-established boundary integral formulation reads

$$ \int_{\Gamma_0} U_{mi}(x, \xi) t_i^c(x) d\Gamma(x) = \int_{\Gamma_0} T_{mi}(x, \xi) [u_i^c(x) - u_i^c(\xi)] d\Gamma(x) \quad (5) $$

where $U_{mi}$ is the well-established fundamental solution of Navier's equations (Kelvin solution) and $T_{mi}$ is the corresponding traction. Their expressions are provided in many textbooks (e.g., [1–3]). Note that Equation (5) is valid for boundary points as well as for internal points. The boundary of the solid is subdivided into N elements and Equation (5) becomes

$$ \sum_{N} \int_{\Gamma^{el}} T_{mi}(x, \xi) [u_{i}^{c}(x) - u_{i}^{c}(\xi)] d\Gamma(x) = \sum_{N} \int_{\Gamma^{el}} t_{i}^{c}(x) U_{mi}(x, \xi) d\Gamma(x) \quad (6) $$

Let $N_i$ be the number of nodes of a boundary element. The usual interpolation over an element leads to $u_i^c(x) = \sum_{a=1}^{N_i} P_a(\eta_1, \eta_2)(u_i^c)^a$ and $t_i^c(x) = \sum_{a=1}^{N_i} P_a(\eta_1, \eta_2)(t_i^c)^a$. The terms $P_a(\eta_1, \eta_2)$ are the interpolation functions of the selected element. Substituting these interpolations into Equation (6), one obtains

$$ c_{mi} u_i^c(\xi) + \sum_{N} \sum_{a=1}^{N_i} \int_{\Gamma^{el}} T_{mi}(x, \xi) P_a(\eta_1, \eta_2) d\Gamma(x) (u_i^c)^a = \sum_{N} \sum_{a=1}^{N_i} \int_{\Gamma^{el}} U_{mi}(x, \xi) P_a(\eta_1, \eta_2) d\Gamma(x) (t_i^c)^a $$

Reorganizing the above relation, with $M$ being the total number of boundary nodes, one obtains

$$ \sum_{b=1}^{M} H_{mi}^{\xi b}(u_{i}^{c})^{b} = \sum_{b=1}^{M} G_{mi}^{\xi b}(t_{i}^{c})^{b} \quad (7) $$

Using Equation (7) for all nodal points $(\xi)$, a system of the following form is obtained:

$$ [H] \{u^c\} = [G] \{t^c\} \quad (8) $$

where $\{u^c\}$ is the vector of nodal displacements and $\{t^c\}$ is the vector of nodal boundary traction. Particular Integral:

We now consider the solution to Equation (4) using the local point collocation method applied to the strong-form differential equations.

### Interpolation

Let us start with a brief reminder of the local point interpolation with a radial basis function. The domain $\Omega$ and its boundary $\Gamma$ are represented by properly scattered collocation centers (nodes), as shown in Figure 1 below.

<figure><img src="image_5.png" /><figcaption>Figure 1. Sample of domain $\Omega$ with scattered collocation centres. Star point with its local influence domain $\Omega_r$.</figcaption></figure>
In the local radial point interpolation method, a field $v$ is approximated on a local sub-domain $\Omega_r$ around a star point $x$ as (see, e.g., [27])

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

with the following constraints:

$$\sum_{i=1}^{N} p_j(x)a_i = 0, \quad j = 1 - M \text{ and } i = 1 - N$$

where $R_i(r)$ is the selected radial basis function, $N$ is the number of collocation centers in the neighborhood (support domain) of point $x$, and $M$ is the number of monomials $p_j(x)$ on the basis of the selected augmented polynomial degree. More precisely, $R_i(r) = R(||x - x^i||)$.

Enforcing approximation (9) to be satisfied at all centers in the support domain, coefficients $a_i$ and $b_j$ are determined by solving a system of equations of the following form:

$$\left\{ \begin{matrix} v_L \\ 0 \end{matrix} \right\} = \begin{bmatrix} R & P \\ P^T & 0 \end{bmatrix} \left\{ \begin{matrix} a \\ b \end{matrix} \right\}$$

where $\{v_L\}$ denotes the vector of nodal values of $v(x)$. Matrices $R$ and $P$ are defined as

$$R = \begin{bmatrix} R(||x^1 - x^1||) & \cdots & R(||x^N - x^1||) \\ \vdots & \ddots & \vdots \\ R(||x^1 - x^N||) & \cdots & R(||x^N - x^N||) \end{bmatrix} \quad \text{and} \quad P = \begin{bmatrix} P_1(x^1) & \cdots & P_M(x^1) \\ \vdots & \ddots & \vdots \\ P_1(x^N) & \cdots & P_M(x^N) \end{bmatrix}^{-1}$$

It can be shown that $\{b\} = [[P]^T[R]^{-1}[P]]^{-1}[P]^T[R]^{-1}\{v_L\} = [F_b]\{v_L\}$ and $\{a\} = [R]^{-1}[[I] - [P][F_b]]\{v_L\} = [F_a]\{v_L\}$

Then, approximation (9) is now rewritten in the matrix form as

$$v(x) = [R_1 \quad R_2 \quad \dots \quad R_N][F_a]\{v_L\} + [P_1 \quad P_2 \quad \dots \quad P_M][F_b]\{v_L\}$$

or more compactly as

$$v(x) = [\Phi(x)] \{v_L\} \qquad (10)$$

### Application

Let us introduce the following vectors:

$$\begin{aligned} \{\sigma^P\}^T &= (\sigma_{11}^P \quad \sigma_{22}^P \quad \sigma_{33}^P \quad \sigma_{12}^P \quad \sigma_{13}^P \quad \sigma_{23}^P) \\ \{T\}^T &= (T_{11} \quad T_{22} \quad T_{33} \quad T_{12} \quad T_{21} \quad T_{13} \quad T_{31} \quad T_{23} \quad T_{32}) \\ \{\varepsilon^P\}^T &= (\varepsilon_{11}^P \quad \varepsilon_{22}^P \quad \varepsilon_{33}^P \quad 2\varepsilon_{12}^P \quad 2\varepsilon_{13}^P \quad 2\varepsilon_{23}^P) \\ \{\nabla\} &= (\partial/\partial X_1 \quad \partial/\partial X_2 \quad \partial/\partial X_3)^T \end{aligned}$$

Note that in terms of particular displacement, $\epsilon_{ij}^p = \frac{1}{2}(u_{i,j}^p + u_{j,i}^p)$.

Given a vector $\{z\} = (z_1 \quad z_2 \quad z_3)^T$, let us define the operators $[B_1(z)] = \begin{pmatrix} z_1 & 0 & 0 & z_2 & z_3 & 0 & 0 & 0 \\ 0 & z_2 & 0 & z_1 & 0 & z_3 & 0 & 0 \\ 0 & 0 & z_3 & 0 & z_1 & z_2 & 0 & 0 \end{pmatrix}$ and $[B_2(z)] = \begin{pmatrix} z_1 & 0 & 0 & z_2 & 0 & z_3 & 0 & 0 \\ 0 & z_2 & 0 & 0 & z_1 & 0 & 0 & z_3 \\ 0 & 0 & z_3 & 0 & 0 & 0 & z_1 & z_2 \end{pmatrix}$.

Using the vectors and operators defined above, Equation (4) is rewritten in the following matrix form:

$$[B_1(\nabla)] \{\sigma^P\} + [B_2(\nabla)] \{T\} = 0 \qquad (11)$$
By adopting interpolation (10) for each component of the vectors {$\sigma^P$} and {$T$}, Equation (11) is written for an internal collocation center in the following, obvious form:

$$[B_1(\nabla)][\Phi_1]\{\sigma_{/el}^P\} + [B_2(\nabla)][\Phi_2]\{T_{/el}\} = 0 \quad (12)$$

In relation (12), $[\Phi_1]$ and $[\Phi_2]$ are properly constructed matrices from interpolation functions defined by (10) and {$z_{/el}$} = ($z_1^1$ $z_2^1$... $z_n^1$ ... $z_1^N$ $z_2^N$... $z_n^N$)T with N = 6 for {$\sigma^P$} and N = 9 for {$T$}.

In Voigt notation, as {$\sigma^P$} = [C] {$\epsilon^P$}, {$\epsilon^P$} = [B_1(\nabla)]{$u^P$}, and {$\sigma^P$} = [C] [B_1(\nabla)]{$u^P$}. Now, adopting interpolation (10) for each component of the particular displacement vector, with $[\Phi_3]$ a properly constructed matrix, Equation (12) becomes

$$[B_1(\nabla)][\Phi_1][C][B_1(\nabla)]^T[\Phi_3]\{u_{/el}^P\} + [B_2(\nabla)][\Phi_2]\{T_{/el}\} = 0 \quad (13)$$

Equation (13) can be rewritten in a compact form as

$$[Q_c]\{u_{/el}^P\} + [\delta Q]\{T_{/el}\} = 0 \quad (14)$$

Collecting Equation (14) for all internal collocation centers leads to

$$[A_B \quad A_I] \begin{Bmatrix} u_B^P \\ u_I^P \end{Bmatrix} + [B] \begin{Bmatrix} T^{NL}(\nabla u) \end{Bmatrix} = 0 \quad (15)$$

In relation (15), indices B and I stand, respectively, for boundary and internal centers and $[A_B \quad A_I]$ is a rectangular matrix. The particular solution $u^P$ is chosen such that it is identically zero at all boundary centers. It then follows that

$$[A_I]\{u_I^P\} + [B]\{T^{NL}(\nabla u)\} = 0$$

Matrix $[A_I]$ is now a square matrix and can be inverted. Thus, the particular integral is obtained as

$$\{u_I^P\} = -[A_I]^{-1}[B]\{T^{NL}(\nabla u)\} = [US]\{T^{NL}(\nabla u)\} \quad (16)$$

Let us now consider the relation between the complementary traction, the particular traction, and the true (nominal) traction. The nominal traction is defined as $P_k = P_{kj} N_j$. With respect to the introduced displacement partition, one has

$$P_k = \sigma_{kj}^C N_j + \sigma_{kj}^P N_j + T_{kj} N_j = t_k^c + t_k^P + T_k$$

By adopting the radial point interpolation (Equation (10)) for the displacement field, after some algebraic manipulations, the following vectors of traction-like densities can be constructed:

$$\{t^P\} = [K_1]\{T^{NL}(\nabla u)\} \quad \text{and} \quad \{T\} = [K_2]\{T^{NL}(\nabla u)\} \quad (17)$$

Final equations:

Considering that {$u^c$} = {$u$} - {$u^P$} and {$t^c$} = {$P$} - {$t^P$} - {$T$}, Equation (8) is put in the following form:

$$[H]\{u - u^P\} = [G]\{P - t^P - T\}$$

Using relations (16) and (17), it becomes
$$[H]\{u\} - [H][US]\{T^{NL}(\nabla u)\} = [G]\{P\} - [G]([K_1] + [K_2])\{T^{NL}(\nabla u)\}$$

More compactly,

$$[H]\{u\} - [R]\{T^{NL}(\nabla u)\} = [G]\{P\} \quad (18)$$

Note that all the matrices in the nonlinear system of Equation (18) are computed only once and stored. What remains now is to apply the real boundary conditions as usual. This is why it was not necessary to specify the boundary conditions associated with Equations (3) and (4). The traction boundary condition may be fixed in magnitude and direction. This type of nonlinear traction is called a follower load and has a known value in the current configuration. Pressure loads are always perpendicular to the boundary and are called rotating loads. Their consideration can be made by decomposing the traction vector as suggested in reference [38], that is, $\{P\} = \{P_0\} + \{P^{NL}(\nabla u)\}$.

Remarkably, the final system of Equation (18) does not depend explicitly on the selected constitutive equation. The latter is involved in the formulation through the nonlinear vector $\{T^{NL}(\nabla u)\}$. In addition, Equation (18) has a form similar to the displacement boundary integral obtained using the FBEM. The results presented in this work are based on the sole use of Equation (18).

Usually, the solution to a geometrically nonlinear mechanical problem is obtained by an incremental loading strategy. Thus, let $\Delta u$ and $\Delta P$ be, respectively, the displacement increment and traction increment between two consecutive loading steps. The incremental form of Equation (16) reads

$$[H]\{\Delta u\} - [R]\{T^{NL}(\nabla u_n + \nabla(\Delta u)) - T^{NL}(\nabla u_n)\} = [G]\{\Delta P\} \quad (19)$$

In relation (19), $u_n$ stands for the displacement obtained at the preceding loading step. The nonlinear system of Equation (19) should be solved by an appropriate iterative method. The conventional Newton-Raphson method with its quadratic rate of convergence can be adopted. However, it is well known the convergence of this method is not guaranteed. It is therefore necessary to write Equation (19) in a form that allows the use of other iterative methods such as the Levenberg-Marquardt method for nonlinear equations [39]. Thus, let us write the vectors of displacement and traction increments as $\{\Delta u\} = [A_u]\{X\} + [B_u]\{\Delta u_d\}$ and $\{\Delta p\} = [A_p]\{X\} + [B_p]\{\Delta p_d\}$, respectively, where $\{X\}$, $\{\Delta u_d\}$, and $\{\Delta p_d\}$ stand, respectively, for the vectors of nodal unknowns, known nodal displacement increments, and given nodal traction increments. Adopting the radial point interpolation (8), the vector of displacement gradient can be written as $\{\nabla u\} = [B]\{u\}$. After some algebraic manipulations, Equation (19) is written in the obvious compact form:

$$[HG]\{X\} - [R]\{T^{NI}(\{K\} + [D]\{X\})\} = \{F\} \quad (20)$$

The new matrices and vectors in Equation (18) are defined as follows:

$$\{K\} = [B](\{u_n\} + [B_u]\{\Delta u_d\});$$

$$[D] = [B][A_u];$$

$$[HG] = [H][A_u] - [G][A_p]$$

$$\{F\} = -[H][B_u]\{\Delta u_d\} - [R]\{T^{NL}(\nabla u_n)\} + [G][B_p]\{\Delta p_d\}$$
A Newton-Raphson strategy or a Levenberg-Marquardt algorithm can be adopted for the solution of (18).

A specific numerical tool is developed using the FORTRAN programming language. The associated pseudocode is as follows:

1/ *Initialisation*: material data, loading steps, etc.

2/ *Conventional BEM phase*

* Construct matrices H and G (cf. Equation (8))

3/ *Radial point interpolation phase*

* Loop over all collocation centers

* Obtain the number of centers in the influence region of point

* Construct interpolation function and its derivatives

* Obtain local matrices of Equation (14)

* Update Equation (15) and prepare Equation (17)

End loop

* ✓ Construct matrices [US], [K₁], [K₂] and [R]

4/ *Apply boundary conditions:*

* Construct matrices [D] and [HG]

5/ *Nonlinear solution phase*

Do until the final load increment

* Construct vectors {K} and {F}

* Solve the nonlinear system

* Update the displacement vector and nominal traction vector

6/ *End program*

What remains now is to investigate the effectiveness of the proposed approach.

## 3. Numerical Examples

The generalized multi-quadrics radial basis functions $R_i(r) = (r^2 + c^2)^q$ are adopted. $r = ||x^i - x||$ is the Euclidian distance between the field point $x$ and the collocation center $x^i$. The constants $c$ and $q$ are known as the shape parameters. The reader is reminded that these parameters are known to affect the accuracy of the solution in methods such as the LRPIM [27,28], the analog equation method [12], and even the MLPG [40]. Accordingly, results from these methods are usually presented with the associated optimal values of these shape parameters.

In order to validate the proposed approach, instead of considering complex deformation states, we choose to use simple geometries and loading cases with different material models (three isotropic materials and two transverse isotropic materials). In fact, we believe that the consideration of complex cases does not allow a clear view of the applicability of a method. In such situations, the solution can only be compared with other available numerical solutions which themselves may be affected by errors. For all the results presented in this work, the nonlinear system of equations is solved by a damped Newton-Raphson method.

### 3.1. Material Models

Five material models are considered.

Subsequently, the following invariants of the right Cauchy-Green tensor are used: $I_1 = \text{tr}(\underline{G})$, $I_2 = \frac{1}{2}(I_1^2 - \text{tr}(\underline{G}^2))$, and $I_3 = \det(\underline{G})$. Given a unit vector $\underline{A}$, one also has $I_4 = \underline{A} \underline{G} \underline{A}$ and $I_5 = \underline{A} \underline{G}^2 \underline{A}$. The derivatives of these invariants with respect to $\underline{G}$ are given as
$$ \frac{\partial I_1}{\partial G} = \underline{1} \frac{\partial I_2}{\partial G} = I_1 \underline{1} - \underline{G} \frac{\partial I_3}{\partial G} = I_3 \underline{G}^{-1} $$

$$ \frac{\partial I_4}{\partial G} = \underline{A} \otimes \underline{A} $$

$$ \frac{\partial I_5}{\partial G} = \underline{A} \otimes (\underline{G} \underline{A}) + (\underline{G} \underline{A}) \otimes \underline{A} $$

### Isotropic Saint Venant–Kirchhoff Model

The simplest nonlinear elastic model is known as the Saint Venant–Kirchhoff solid. The second Piola–Kirchhoff stress tensor is expressed as

$$ \underline{S} = \lambda \operatorname{tr}(\underline{E}) \underline{1} + 2\mu \underline{E} \qquad (21) $$

where $\underline{E} = (\underline{G} - \underline{1})/2$ is the Green–Lagrange strain tensor, and $\lambda$ and $\mu$ are the Lamé constants.

This model has been considered by Prieto et al. [6] in a plane strain context to prove the effectiveness of a field boundary element implementation. Kompis et al. [41] used the same model to validate their implementation of a Treff boundary element.

### Isotropic Neo–Hookean Model

We adopt the Neo–Hookean solid considered in reference [42] for which the stress tensor $\underline{S}$ is given by

$$ \underline{S} = \lambda \ln J \underline{G}^{-1} + \mu (\underline{1} - \underline{G}^{-1}) \qquad (22) $$

In Equation (22), $\underline{G}^{-1}$ is the inverse of $\underline{G}$ and $J = \det(F)$.

The considered Neo–Hookean model has been used by Gu et al. [42] to show the efficiency of a meshless local Kriging method.

### Isotropic Mooney–Rivlin Model

In this work, a modified Mooney–Rivlin material [43] is considered and the second Piola–Kirchhoff stress tensor is expressed as

$$ \underline{S} = C_1 \underline{1} + C_2 \underline{G} + C_3 \underline{G}^{-1} \qquad (23) $$

$$ C_1 = 2 A_{10} I_3^{-1/3} + 2 A_{01} I_3^{-2/3} I_1 $$

$$ C_2 = -2 A_{01} I_3^{-2/3} $$

$$ C_3 = - \left( \frac{2}{3} A_{10} I_3^{-1/3} I_1 + \frac{4}{3} A_{01} I_3^{-2/3} I_2 - K (J-1) I_3^{1/2} \right) $$

where the material parameters $A_{10}$ and $A_{01}$ are related to the shear modulus by $\mu = 2(A_{01} + A_{10})$, and K is the bulk modulus. The model has been applied in the solution of a large-deformation contact analysis of elastomers by Hu et al. [44].

### The Transversely Isotropic Saint Venant–Kirchhoff Model

The first transversely isotropic model considered in this work is of the Saint Venant–Kirchhoff type. The model is presented in the work by Bonet and Burton [45]. The second Piola–Kirchhoff stress tensor is expressed as the sum of two terms. The first term corresponds to the isotropic parts and is similar to that in Equation (21). The second term accounts for the anisotropy and is given by

$$ \underline{S}_{itr} = 2b(I_4 - 1) \underline{1} + 2(a + b(I_1 - 3) + 2c(I_4 - 1)) \underline{A} \otimes \underline{A} - c(\underline{G}\underline{A} \otimes \underline{A} + \underline{A} \otimes \underline{G}\underline{A}) \qquad (24) $$
In the relations above, $\underline{A}$ is the unit vector in the direction perpendicular to the plane of isotropy.

Considering that the elasticity tensor must be that of a transversely isotropic medium in the small-strain regime, the material parameters $a, b, c$ can be identified.

When the axis of isotropy coincides with the Z-axis, the five engineering material parameters are $E_x, \nu, E_z, \nu_1$, and $G_{xz}$. It can be shown that

$$a = \frac{E_x}{2(1 + \nu)} - G_{xz}, \quad b = \frac{E_x}{4(1 + \nu)\Delta}(\nu + m\nu_1(\nu_1 - \nu - 1)),$$

and

$$c = \frac{E_x}{8(1 + \nu)\Delta}(1 - m(1 + \nu_1^2 - \nu)) + \frac{a}{2} - b.$$

In these expressions, $\Delta = \nu - 1 + 2mv_1^2$ and $m = \frac{E_z}{E_x}$.

## Neo-Hookean Transversely Isotropic Model

The second transversely isotropic model considered is a Neo-Hookean solid also presented in the work by Bonet and Burton [45]. The second Piola-Kirchhoff stress tensor is again expressed as a sum of two terms. The first term corresponds to the isotropic parts and is similar to Equation (22). The second term accounts for the anisotropy and is given by Equation (25). Considering that in the small-strain regime, the right Cauchy-Green tensor is close to identity, approximate expressions of the model parameters in this term are those of the preceding model.

$$\underline{S}_{itr} = 2b(I_4 - 1)\underline{1} + 2(a + 2b \ln J + 2c(I_4 - 1))\underline{A} \otimes \underline{A} - c(\underline{G}\underline{A} \otimes \underline{A} + \underline{A} \otimes \underline{G}\underline{A}) \quad (25)$$

### 3.2. Case of Uni-Axial Loading

Let us first consider the case of uniaxial loading. This academic case is considered in order to assess the effectiveness of the proposed solution method and to investigate the impact of the multi-quadrics shape parameters on the solution's accuracy.

#### 3.2.1. Unit Cube under Tension

A unit cube is simply supported at its lower end ($Z = 0$). The upper surface ($Z = 1$) is uniformly loaded (see Figure 2). The other faces of the cube are traction-free. The boundary of the cube is subdivided into twenty-four nine-noded quadrilateral elements, that is, four elements per face.

<figure><img src="image_6.png" /><figcaption>Figure 2. Geometry of the simply supported (triangles) unit cube under uniform load (arrows) with associated boundary discretization.</figcaption></figure>
Under such a uniform loading state, the displacement field is given by $u_X = (\lambda_1 - 1)X$, $u_Y = (\lambda_1 - 1)Y$, and $u_Z = (\lambda_3 - 1)Z$. Following this, the non-zero components of the deformation gradient tensor, the Green strain tensor ($G$) and its inverse ($G^{-1}$), and the Green–Lagrange tensor ($E$) are gathered in Table 1 below.

<table><caption>**Table 1.** Unidirectional loading: non-zero components of the deformation gradient, Green strain, inverse Green strain, and Green–Lagrange tensors.</caption><thead><tr><th>F<sub>XX</sub> = F<sub>YY</sub></th><th>F<sub>ZZ</sub></th><th>G<sub>XX</sub> = G<sub>YY</sub></th><th>G<sub>ZZ</sub></th><th>G<sub>XX</sub><sup>-1</sup> = G<sub>YY</sub><sup>-1</sup></th><th>G<sub>ZZ</sub><sup>-1</sup></th><th>E<sub>XX</sub> = E<sub>YY</sub></th><th>E<sub>ZZ</sub></th></tr></thead><tbody><tr><td>λ<sub>1</sub></td><td>λ<sub>3</sub></td><td>λ<sub>1</sub><sup>2</sup></td><td>λ<sub>3</sub><sup>2</sup></td><td>1/λ<sub>1</sub><sup>2</sup></td><td>1/λ<sub>3</sub><sup>2</sup></td><td>(λ<sub>1</sub><sup>2</sup> - 1)/2</td><td>(λ<sub>3</sub><sup>2</sup> - 1)/2</td></tr></tbody></table>

The specimen is loaded either by applying the traction at the top surface ($P_{ZZ}$ known) or by imposing the displacement of points of this surface ($\lambda_3$ known). From the numerical results, one can assess the values of $\lambda_3$ or $P_{ZZ}$ and $\lambda_1$. Their analytical values can be calculated using the relations below where $E = E_x$, $\mu = \frac{E_x}{2(1+\nu)}$, and $\nu$ are, respectively, Young's modulus, the shear modulus, and Poisson's ratio. In the transversely isotropic case, these data are complemented by $E_z$, $\nu_1$, and $G_{xz}$:

**Saint Venant–Kirchhoff Solid:**

$$P_{ZZ} = \frac{E}{2} \lambda_3 (\lambda_3^2 - 1) \\ (\lambda_1^2 - 1) = -\nu (\lambda_3^2 - 1)$$

**Neo–Hookean Solid:**

$$\frac{2\nu}{1-2\nu} \ln\left(1-\frac{S_{ZZ}}{\mu}\right) + \left(1-\frac{S_{ZZ}}{\mu}\right)\lambda_3^2 + \frac{3\nu}{1-2\nu}\ln(\lambda_3^2) = 1 \\ \lambda_1^2 = \left(1-\frac{S_{ZZ}}{\mu}\right)\lambda_3^2$$

**Mooney–Rivlin Solid:**

$$\lambda_3^2 S_{ZZ} = C_1(\lambda_3^2 - \lambda_1^2) + C_2(\lambda_3^4 - \lambda_1^4) \\ (C_1 + C_2 \lambda_1^2) \lambda_1^2 + C_3 = 0$$

In our calculations, we used $\mu = 5$ MPa, $\nu = 0.2$, $A_{10} = 0$, $A_{01} = 2.5$ MPa, and $K = 20/3$.

**Transversely Isotropic Saint Venant–Kirchhoff Solid:**

$$\lambda_1^2 = 1 + \frac{\lambda + 4b}{2(\lambda + \mu)} (1 - \lambda_3^2) \\ P_{ZZ} = (\lambda + 4b)(\lambda_1^2 - 1) + (\lambda + 2\mu + 8b + 8c - 4a) \frac{(\lambda_3^2 - 1)}{2}$$

**Transversely Isotropic Neo–Hookean Solid:**

$$2\lambda \ln(\lambda_1) + (\mu + 2b(\lambda_3^2 - 1))\lambda_1^2 - \mu = 0 \\ \lambda_3^2 S_{ZZ} = \mu \ln(\lambda_1^2 \lambda_3) + \mu(\lambda_3^2 - 1) + (2a - 8b - 8c)\lambda_3^2 + 4b\lambda_1^2 \lambda_3^2 + 2(2b + 2c - a)\lambda_3^4$$
For the results presented hereafter, the 150 boundary nodes of the unit cube are supplemented with 27 internal collocation centers. In Figure 3, the evolutions of the calculated nominal tractions $P_{zz}$ at the top surface of the cube with respect to the applied stretch $\lambda_3$ are shown. It is observed that the load required to achieve a certain amount of elongation for the Saint Venant–Kirchhoff solid is almost seven times that required for the two other isotropic models. The variations in the lateral stretch $\lambda_1$ with respect to $\lambda_3$ are shown in Figure 4. The results are practically analytical solutions. Let us point out that in all cases, the entire load is applied in one step. The average number of Newton iterations is less than 10. The iterative process is stopped when the infinity norm of the difference between two consecutive iterates is less than $\varepsilon = 10^{-3}$.

<figure><img src="image_7.png" /><figcaption>Figure 3. Calculated nominal traction versus applied axial stretch.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. Calculated lateral stretch versus axial stretch.</figcaption></figure>

For the transversely isotropic case, two sets of engineering material parameters are adopted. They are presented in Table 2 below.
<table><caption>**Table 2.** Engineering material parameters of the transversely isotropic solids.</caption><thead><tr><th>Material</th><th><i>E</i><sub>x</sub></th><th><i>E</i><sub>z</sub></th><th><i>ν</i></th><th><i>ν</i><sub>1</sub></th><th><i>G</i><sub>xz</sub></th></tr></thead><tbody><tr><td>Set 1</td><td>7.22</td><td>23</td><td>0.14</td><td>0.08789</td><td>4.20</td></tr><tr><td>Set 2</td><td>4.71</td><td>13.7</td><td>0.14</td><td>0.0997</td><td>1.49</td></tr></tbody></table>

These data represent the mechanical properties of a unidirectional nano-reinforced laminate composite. The first set contains the mechanical properties of 11% carbon nanotubes in a polymeric matrix and the second is for 11% carbon nanotubes and 11% carbon fibers [46]. Using set 1, and assuming that under large strains, the material is of the Saint-Venant type, the sample is denoted SVK_T1. If it is assumed to behave like a Neo-Hookean solid, it is denoted as NH_T1. Similar notation is adopted for the second set of engineering constants (SVK_T2 and NH_T2).

The traction loads corresponding to $\lambda_3 = 1.225, 1.45, 1.675$, and $1.9$ were determined analytically. This traction was then applied at the top surface of the specimen. The calculated numerical values of $\lambda_3$ were practically analytical solutions. The computed values of $\lambda_1$ are compared to the analytical values in Table 3 below. It can be observed that, in the case of anisotropic material, the method gives accurate results.

<table><caption>**Table 3.** Comparison of analytical and numerical stretch $\lambda_1$ for various values of $\lambda_3$ in the cases of transversely isotropic cube under tension load.</caption><thead><tr><th rowspan="3">λ<sub>3</sub></th><th colspan="8">λ<sub>1</sub></th></tr><tr><th colspan="2">SVK_T1</th><th colspan="2">SVK_T2</th><th colspan="2">NH_T1</th><th colspan="2">NH_T2</th></tr><tr><th>Anal</th><th>Num</th><th>Anal</th><th>Num</th><th>Anal</th><th>Num</th><th>Anal</th><th>Num</th></tr></thead><tbody><tr><td>1.225</td><td>0.9273</td><td>0.92726</td><td>0.92456</td><td>0.9242</td><td>0.9359</td><td>0.93592</td><td>0.93356</td><td>0.93356</td></tr><tr><td>1.45</td><td>0.8315</td><td>0.83144</td><td>0.82479</td><td>0.8226</td><td>0.8682</td><td>0.868178</td><td>0.8632</td><td>0.86324</td></tr><tr><td>1.675</td><td>0.7032</td><td>0.70318</td><td>0.6902</td><td>0.6884</td><td>0.7952</td><td>0.79524</td><td>0.7874</td><td>0.78746</td></tr><tr><td>1.9</td><td>0.5189</td><td>0.51884</td><td>0.4931</td><td>0.4894</td><td>0.7162</td><td>0.716116</td><td>0.7054</td><td>0.7054</td></tr></tbody></table>

In the case of uniaxial loading, Equation (18) is efficiently solved by the Newton-Raphson method.

When radial basis functions are adopted, the numerical results are usually supplemented with optimal values of the shape parameters. The multi-quadrics shape parameters $c$ and $q$ vary, respectively, in the ranges [0.0005, 0.05] and [0.5; 1.5]. The numerical results remain unchanged for any couple of parameters in the mentioned ranges. In view of the foregoing, it is concluded that the proposed approach is effective and accurate in the case of unidirectional loading.

### 3.2.2. Cases of a Cylindrical Specimen under Uniaxial Loading

Now, we aim to check the effectiveness of the proposed approach in the case of curved boundaries. The considered geometry is a cylinder of height 2 mm with a circular cross-section of radius 1 mm. The boundary of the cylinder is subdivided into 80 nine-node quadrilateral elements and its boundary nodes are supplemented by 225 internal collocation centers. In the cylindrical coordinates system, the analytical solution is similar to that presented in the preceding paragraph with $(X, Y, Z) \equiv (R, \theta, Z)$.

The cylinder specimen was loaded in compression in order to shorten its length by 22.5% on the one hand and by 45% on the other hand. The total load was applied in one step. Some of the calculated results are compared to their analytical counterparts in Table 4 below. It can be observed that highly accurate results are obtained for 22.5% compression. For 45% compression, the high accuracy is maintained only for the Neo-Hookean solid.
<table><caption>**Table 4.** Comparison of analytical and numerical stretches $\lambda_1$ and $\lambda_3$ in the cases of compression of a cylindrical solid.</caption><thead><tr><th></th><th colspan="2">Saint Venant-Kirchhoff</th><th colspan="2">Neo-Hookean</th><th colspan="2">SVK_T1</th><th colspan="2">NH_T1</th></tr></thead><tbody><tr><td>P<sub>33</sub></td><td colspan="2">-1.8571</td><td colspan="2">-3.24745</td><td colspan="2">-3.09577</td><td colspan="2">-4.018277</td></tr><tr><td></td><td>Anal</td><td>Num</td><td>Anal</td><td>Num</td><td>Anal</td><td>Num</td><td>Anal</td><td>Num</td></tr><tr><td>λ<sub>3</sub></td><td>0.775</td><td>0.775</td><td>0.775</td><td>0.775</td><td>0.775</td><td>0.775</td><td>0.775</td><td>0.7625</td></tr><tr><td>λ<sub>1</sub></td><td>1.03917</td><td>1.03918</td><td>1.0507</td><td>1.0507</td><td>1.054</td><td>1.05436</td><td>1.063</td><td>1.0543</td></tr><tr><td>P<sub>33</sub></td><td colspan="2">-2.30175</td><td colspan="2">-8.612425</td><td colspan="2">-3.837</td><td colspan="2">-7.9176</td></tr><tr><td></td><td>Anal</td><td>Num</td><td>Anal</td><td>Num</td><td>Anal</td><td>Num</td><td>Anal</td><td>Num</td></tr><tr><td>λ<sub>3</sub></td><td>0.55</td><td>0.6044</td><td>0.55</td><td>0.55</td><td>0.55</td><td>0.6044</td><td>0.55</td><td>0.5152</td></tr><tr><td>λ<sub>1</sub></td><td>1.0675</td><td>1.06159</td><td>1.1179</td><td>1.1179</td><td>1.0933</td><td>1.08523</td><td>1.1298</td><td>1.1238</td></tr></tbody></table>

### 3.3. Simple Shear of a Cubic Specimen

Special attention was given to the problem of simple shear in [47,48].

Under a tri-axial stretch and simple shear, the displacement field is given by

$$U_X = (\lambda_1 - 1)X + k Y; \quad U_Y = (\lambda_2 - 1)Y; \quad U_Z = (\lambda_3 - 1)Z$$

$\lambda_1$, $\lambda_2$, and $\lambda_3$ are the stretches and $k$ is the shear factor. When $\lambda_1 = \lambda_2 = \lambda_3 = 1$, the deformation is that of the simple shear. The pure shear corresponds to the situation when the Cauchy stress tensor in the deformed configuration has the form

$$\sigma = \begin{pmatrix} 0 & s & 0 \\ s & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}$$

A number of simulations were carried out on the unit cube with Dirichlet boundary conditions modeling the simple shear. More specifically, the displacements of all boundary nodes are specified using $U_X = k Y$; $U_Y = 0$; and $U_Z = 0$. The accuracy of the proposed approach was measured by comparing the analytical and numerical values of the Cauchy stress tensor.

Results of the non-zero Cauchy stress components for shear factor $k = 0.25$ are given in Table 5 below.

<table><caption>**Table 5.** Non-zero components of the Cauchy stress tensor under simple shear of a cubic specimen (shear factor of 0.25).</caption><thead><tr><th></th><th colspan="8">Model</th></tr><tr><th></th><th></th><th>Saint-Venant</th><th>Neo-Hookean</th><th>Mooney-Rivlin</th><th>SV_T1</th><th>SV_T2</th><th>NH_T1</th><th>NH_T2</th></tr></thead><tbody><tr><th rowspan="4">Stress component</th><th>σ<sub>11</sub></th><td>0.755</td><td>0.3125</td><td>0.1042</td><td>0.436</td><td>0.2854</td><td>0.1979</td><td>0.1291</td></tr><tr><th>σ<sub>22</sub></th><td>0.417</td><td>0</td><td>-0.208</td><td>0.238</td><td>0.1563</td><td>0</td><td>0</td></tr><tr><th>σ<sub>33</sub></th><td>0.104</td><td>0</td><td>0.1042</td><td>0.078</td><td>0.0532</td><td>0.0377</td><td>0.0260</td></tr><tr><th>σ<sub>12</sub></th><td>1.354</td><td>1.25</td><td>1.25</td><td>0.792</td><td>0.5164</td><td>0.7917</td><td>0.5164</td></tr></tbody></table>

These results are practically analytical solutions. They show that in order to obtain a state of pure shear, it is necessary to stretch the specimen. With this case, it is shown that the proposed approach allows an accurate prediction of the Poynting effect. This nonlinear effect corresponds to the development of stress normal to the sheared faces of a cubic specimen.
### 3.4. Constrained Tension of a Cubic Specimen

The cross-section of the considered sample is presented in Figure 5.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Cross-section of the constrained tension specimen.</figcaption></figure>

The specimen is constrained against displacement at its lower end and subjected to a uniform displacement in the z-direction at its upper end. The remaining faces of the sample are free of traction. According to Foroutan et al. [49], the solution of this problem using the finite element method is practically impossible due to large mesh distortions. These authors solved the 2D version of the problem using the RKPM method. The LPI-BEM was applied to the solution of the 3D version of the problem, adopting the Saint Venant–Kirchhoff and Neo–Hookean models in their isotropic and transverse isotropic versions. The boundary of the unit cube was subdivided into 96 quadrilateral elements. The boundary nodes were supplemented with 27 internal collocation centers. Elongation in the z-direction was applied incrementally.

The calculated deformed shapes at 22.5% and 67.5% elongations are shown in Figures 6 and 7, respectively, for the four cases considered.

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Deformed shape of the unit cube after 22.5% elongation for different material behavior. (a) Neo–Hookean, (b) transversely isotropic Neo–Hookean, (c) Saint Venant–Kirchhoff, (d) transversely isotropic Saint Venant–Kirchhoff.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Deformed shape of the unit cube after 67.5% elongation for different material behavior. (a) Neo-Hookean, (b) transversely isotropic Neo-Hookean, (c) Saint Venant-Kirchhoff, (d) transversely isotropic Saint Venant-Kirchhoff.</figcaption></figure>

It can be observed that the contraction of the sample follows the results obtained in the case of simple tension. Indeed, for the same value of stretch $\lambda_3$, the smaller the value of stretch $\lambda_1$, the greater the lateral contraction. This result shows that the proposed approach is valid in this case too.

### 3.5. Bending of a Prismatic Bar

As a final example, let us consider the case of the bending of a bar. In this case, there is a region in the bar with a high variation in displacement gradient. Furthermore, the displacement gradient is not uniform along the bar. The prismatic bar considered in our calculation has a square cross-section (1 mm × 1 mm) and 10 mm length. It is clamped at one end and submitted at the other end to a tangential loading, as shown in Figure 8. The boundary of the bar is subdivided into 48 quadrilateral elements. The results presented hereafter are obtained with 171 regularly spaced internal collocation centers.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Geometry of the bending specimen in the initial configuration.</figcaption></figure>
In this case, with a localized region of strong geometrical nonlinearities, the Newton-Raphson method sometimes fails at a certain load increment during the loading process. The presented results are obtained by using the Levenberg-Marquardt method. The deformation of the mean line of the specimen is represented in Figures 9 and 10 for the Saint Venant-Kirchhoff solid and the Mooney-Rivlin solid.

<figure><img src="image_13.png" /><figcaption>Figure 9. Displacement of the center line of a Saint Venant-Kirchhoff material under bending.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10. Displacement of the center line of a Mooney-Rivlin material under bending.</figcaption></figure>

In this case, the validity of the approach is assessed by comparing the results obtained with two different loading steps. As can be observed, in this case also, the method provides stable results.
## 4. Conclusions

Considering a small-strain isotropic elastic solid, the conventional BEM leads to highly accurate results. In the case of geometrically nonlinear elasticity, the field boundary element method is effective. In this case, the main appeal of the BEM (reduction in dimension by 1) is tarnished to some extent. In this work, the approach called LPI-BEM, which has already been proven accurate for anisotropic, nonlocal, and piezoelectric elasticity, is extended to the case of geometrically nonlinear elasticity. The method couples the conventional isotropic BEM to local radial point interpolation applied to strong-form differential equations. The implementation of the method requires only few modifications to the existing BEM code. The solution to nonlinear elasticity problems using the LPI-BEM is demonstrated on some examples that have analytical solutions. Note that in most of the considered examples, the boundary mesh is coarse and the number of internal collocation centers is small. The method has successfully solved the problem of constrained tension. This simple-to-implement approach seems promising. The next step will be the consideration of more complex geometries as well as dynamic loading.

**Author Contributions:** Conceptualization, R.K.N., B.K.B. and A.R.K.; methodology, R.K.N.; software, A.R.K.; writing—original draft preparation, A.R.K.; writing—review and editing, R.K.N. and B.K.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available in the article.

**Acknowledgments:** The author Aly Rachid Korbeogo thanks the French Embassy in Burkina Faso for funding an internship program which allowed him to travel to France and work in a team on this work currently being published.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Balas, J.; Sladek, J.; Sladek, V. *Stress Analysis by Boundary Element Method*; Elsevier: Berlin/Heidelberg, Germany, 1989.

2. Bonnet, M. *Boundary Integral Equation Methods for Solids and Fluids*; John Wiley & Sons: New York, NY, USA, 1999. Available online: http://books.google.fr/books?id=3FBjQgAACAAJ (accessed on 20 October 2023).

3. Brebbia, C.A.; Dominguez, J. *Boundary Elements: An Introductory Course, Computational Mechanics Publications*; WIT Press: Southampton, UK, 1992; Available online: http://books.google.fr/books?id=YBEIJxdldHYC (accessed on 20 October 2023).

4. Foerster, A.; Kuhn, G. A field boundary element formulation for material nonlinear problems at finite strains. *Int. J. Solids Struct.* **1994**, *31*, 1777–1792. [CrossRef]

5. Al-Gahtani, H.J.; Altiero, N.J. Application of the boundary element method to rubber-like elasticity. *Appl. Math. Model.* **1996**, *20*, 654–661. [CrossRef]

6. Prieto, I.; Ibán, A.L.; Garrido, J.A. 2D analysis for geometrically non-linear elastic problems using the BEM. *Eng. Anal. Bound. Elem.* **1999**, *23*, 247–256. [CrossRef]

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

8. 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]

9. Chen, L.; Kassab, A.J.; Nicholson, D.W.; Chopra, M.B. Generalized boundary element method for solids exhibiting nonhomogeneities. *Eng. Anal. Bound. Elem.* **2001**, *25*, 407–422. [CrossRef]

10. Ochiai, Y. Meshless large plastic deformation analysis considering with a friction coefficient by triple-reciprocity boundary element method. *Int. J. CMEM* **2018**, *6*, 989–999. [CrossRef]

11. Zhu, T.; Zhang, J.-D.; Atluri, S.N. A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach. *Comput. Mech.* **1998**, *21*, 223–235. [CrossRef]

12. Katsikadelis, J.T.; Nerantzaki, M.S. The boundary element method for nonlinear problems. *Eng. Anal. Bound. Elem.* **1999**, *23*, 365–373. [CrossRef]

13. Brun, M.; Capuani, D.; Bigoni, D. A boundary element technique for incremental, non-linear elasticity: Part I: Formulation. *Comput. Methods Appl. Mech. Eng.* **2003**, *192*, 2461–2479. [CrossRef]

14. Brun, M.; Bigoni, D.; Capuani, D. A boundary element technique for incremental, non-linear elasticity: Part II: Bifurcation and shear bands. *Comput. Methods Appl. Mech. Eng.* **2003**, *192*, 2481–2499. [CrossRef]

15. Liu, G.R. *Meshfree Methods: Moving Beyond the Finite Element Method*, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2010.
16. Belytschko, T.; Lu, Y.Y.; Gu, L. Element-free Galerkin methods. *Int. J. Numer. Methods Eng.* **1994**, *37*, 2. [CrossRef]

17. Chen, J.-S.; Pan, C.; Wu, C.-T.; Liu, W.K. Reproducing Kernel Particle Methods for large deformation analysis of non-linear structures. *Comput. Methods Appl. Mech. Eng.* **1996**, *139*, 195–227. [CrossRef]

18. Chen, J.-S.; Pan, C.; Wu, C.-T. Large deformation analysis of rubber based on a reproducing kernel particle method. *Comput. Mech.* **1997**, *19*, 3. [CrossRef]

19. Zhang, X.; Yao, Z.; Zhang, Z. Application of MLPG in Large Deformation Analysis. *Acta Mech. Sin.* **2006**, *22*, 331–340. [CrossRef]

20. Liew, K.M.; Ng, T.Y.; Wu, Y.C. Meshfree method for large deformation analysis—a reproducing kernel particle approach. *Eng. Struct.* **2002**, *24*, 543–551. [CrossRef]

21. Areias, P.; Fernandes, J.L.M.; Rodrigues, H.C. Galerkin-based finite strain analysis with enriched radial basis interpolation. *Comput. Methods Appl. Mech. Eng.* **2022**, *394*, 114873. [CrossRef]

22. Areias, P.; Fernandes, L.; Melicio, R. Moving least-squares in finite strain analysis with tetrahedra support. *Eng. Anal. Bound. Elem.* **2022**, *139*, 1–13. [CrossRef]

23. Röthlin, M.; Klippel, H.; Wegener, K. Meshless Methods for Large Deformation Elastodynamics. *arXiv* **2018**. [CrossRef]

24. Wu, C.T.; Hu, W.; Chen, J.S. A meshfree-enriched finite element method for compressible and near-incompressible elasticity. *Numer. Methods Eng.* **2012**, *90*, 882–914. [CrossRef]

25. Zhou, B.; Zhang, C.; Zhao, F. A Finite Element-Meshless Hybrid Method (FEMLHM) of Elasticity Problem and Its Applications. *Mech. Solids* **2023**, *58*, 852–871. [CrossRef]

26. Yu, S.Y.; Peng, M.J.; Cheng, H.; Cheng, Y.M. The improved element-free Galerkin method for three-dimensional elastoplasticity problems. *Eng. Anal. Bound. Elem.* **2019**, *104*, 215–224. [CrossRef]

27. 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]

28. Liu, G.R.; Gu, Y.T. A meshfree method: Meshfree weak?strong (MWS) form method, for 2-D solids. *Comput. Mech.* **2003**, *33*, 2–14. [CrossRef]

29. Lee, S.-H.; Yoon, Y.-C. Meshfree point collocation method for elasticity and crack problems. *Int. J. Numer. Methods Eng.* **2004**, *61*, 22–48. [CrossRef]

30. Zheng, H.; Lai, X.; Hong, A.; Wei, X. A Novel RBF Collocation Method Using Fictitious Centre Nodes for Elasticity Problems. *Mathematics* **2022**, *10*, 3711. [CrossRef]

31. 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*, 4. [CrossRef]

32. 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]

33. Thurieau, N.; Njiwa, R.K.; Taghite, M. A simple solution procedure to 3D-piezoelectric problems: Isotropic BEM coupled with a point collocation method. *Eng. Anal. Bound. Elem.* **2012**, *36*, 11. [CrossRef]

34. 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. [CrossRef]

35. Ciarlet, P.G. *Mathematical Elasticity: Three-Dimensional Elasticity*, Society for Industrial & Applied Mathematics (SIAM): Philadelphia, PA, USA, 2021. [CrossRef]

36. Maugin, G.A. On canonical equations of continuum thermomechanics. *Mech. Res. Commun.* **2006**, *33*, 705–710. [CrossRef]

37. Park, K.-H.; Banerjee, P.K. Two- and three-dimensional transient thermoelastic analysis by BEM via particular integrals. *Int. J. Solids Struct.* **2002**, *39*, 2871–2892. [CrossRef]

38. Karami, G.; Derakhshan, D. A field boundary element for large deformation analysis of hyperelastic problems. *Sci. Iran.* **2001**, *2*, 110–122.

39. Yu, Z. On the global convergence of a Levenberg-Marquardt method for constrained nonlinear equations. *J. Appl. Math. Comput.* **2004**, *16*, 183–194. [CrossRef]

40. Hu, D.; Sun, Z. The Meshless Local Petrov-Galerkin Method for Large Deformation Analysis of Hyperelastic Materials. *ISRN Mech. Eng.* **2011**, *2011*, 967512. [CrossRef]

41. Kompiš, V.; Toma, M.; Žmindák, M.; Handrik, M. Use of Trefftz functions in non-linear BEM/FEM. *Comput. Struct.* **2004**, *82*, 2351–2360. [CrossRef]

42. Gu, Y.T.; Wang, Q.X.; Lam, K.Y. A meshless local Kriging method for large deformation analyses. *Comput. Methods Appl. Mech. Eng.* **2007**, *196*, 1673–1684. [CrossRef]

43. Crisfield, M.A. *Non-Linear Finite Element Analysis of Solids and Structures*, 2nd ed.; John Wiley & Sons, Ltd.: Chichester, NY, USA, 1991.

44. Hu, D.A.; Long, S.Y.; Han, X.; Li, G.Y. A meshless local Petrov-Galerkin method for large deformation contact analysis of elastomers. *Eng. Anal. Bound. Elem.* **2007**, *31*, 657–666. [CrossRef]

45. Bonet, J.; Burton, A.J. A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations. *Comput. Methods Appl. Mech. Eng.* **1998**, *162*, 151–164. [CrossRef]

46. Kulkarni, M.; Carnahan, D.; Kulkarni, K.; Qian, D.; Abot, J.L. Elastic response of a carbon nanotube fiber reinforced polymeric composite: A numerical and experimental study. *Compos. Part B Eng.* **2010**, *41*, 414–421. [CrossRef]

47. Destrade, M.; Murphy, J.; Saccomandi, G. Simple shear is not so simple. *Int. J. Non-Linear Mech.* **2012**, *47*, 210–214. [CrossRef]
48. Mihai, L.A.; Goriely, A. Numerical simulation of shear and the Poynting effects by the finite element method: An application of the generalised empirical inequalities in non-linear elasticity. *Int. J. Non-Linear Mech.* **2013**, *49*, 1–14. [CrossRef]

49. Foroutan, M.; Dalayeli, H.; Sadeghian, M. Simulation of Large Deformations of Rubbers by the RKPM Method. *Int. J. Math. Comput. Phys. Electr. Comput. Eng.* **2007**, *1*, 136–141.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.