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

: The attractiveness of the boundary element method—the reduction in the problem dimen-sion 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


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][2][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][5][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 largestrain 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 strains with material nonlinearities (e.g., [16][17][18][19][20][21][22][23][24][25][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][28][29][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, multiphysics, and multi-field problems (e.g., [31][32][33][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 nonlinear 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 uniaxial 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.

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 Ω 0 with boundary Γ 0 .The position vector of a particle is denoted by .After deformation, the body occupies the domain Ω and the position vector of a particle is denoted by .It is assumed that there is a twice differentiable invertible map  that takes Ω 0 to Ω, that is,  = ().The displacement vector  is given by A measure of the deformation is described by the deformation gradient  relative to  given by The identity tensor is denoted as 1.The right Cauchy-Green tensor () can be found from  according to  =   .
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]) The tensor  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  is expressed in terms of the symmetric second Piola-Kirchhoff stress tensor  by   =     .In the case of a hyper-elastic law,  is defined from an elastic potential  by  = 2   .Note that Equation ( 1) is a highly nonlinear partial differential equation.

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.
In the above expression,   is a fourth-order tensor of the elastic constants of an isotropic homogeneous material.Now, let us assume that the displacement  can be viewed as the sum of two terms: a complementary term (  ) which behaves as the displacement vector in small-strain deformation of an isotropic elastic solid, and a particular term (  ) which accounts for all nonlinear effects of the deformation process.Bearing this in mind, Equation (1) is rewritten as In Equation (2), we introduce the following: where    is the complementary stress tensor,    is the particular stress tensor, and   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 Naviertype differential equation: The associated boundary traction is where   is the well-established fundamental solution of Navier's equations (Kelvin solution) and   is the corresponding traction.Their expressions are provided in many textbooks (e.g., [1][2][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 .The terms   ( 1 ,  2 ) are the interpolation functions of the selected element.Substituting these interpolations into Equation ( 6), one obtains Reorganizing the above relation, with  being the total number of boundary nodes, one obtains Using Equation ( 7) for all nodal points (), a system of the following form is obtained: where {  } is the vector of nodal displacements and {  } 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 Ω and its boundary Γ are represented by properly scattered collocation centers (nodes), as shown in Figure 1 below.In the local radial point interpolation method, a field  is approximated on a local sub-domain Ω  around a star point  as (see, e.g., [27]) with the following constraints: where   () is the selected radial basis function,  is the number of collocation centers in the neighborhood (support domain) of point , and  is the number of monomials   () on the basis of the selected augmented polynomial degree.More precisely,   () = (‖ −   ‖. Enforcing approximation (9) to be satisfied at all centers in the support domain, coefficients   and   are determined by solving a system of equations of the following form: where { / } denotes the vector of nodal values of ().Matrices R and P are defined as

It can be shown that
Then, approximation ( 9) is now rewritten in the matrix form as or more compactly as
Using the vectors and operators defined above, Equation ( 4) is rewritten in the following matrix form: By adopting interpolation (10) for each component of the vectors {  } and {} , Equation ( 11) is written for an internal collocation center in the following, obvious form: In relation (12), [Φ 1 ] and [Φ 2 ] are properly constructed matrices from interpolation functions defined by (10) and .    )  with  = 6 for {  } and  = 9 for {}.
In Voigt notation, as (10) for each component of the particular displacement vector, with [Φ 3 ] a properly constructed matrix, Equation ( 12) becomes Equation ( 13) can be rewritten in a compact form as Collecting Equation ( 14) for all internal collocation centers leads to In relation (15), indices B and I stand, respectively, for boundary and internal centers and [    ] is a rectangular matrix.The particular solution   is chosen such that it is identically zero at all boundary centers.It then follows that Matrix [  ] is now a square matrix and can be inverted.Thus, the particular integral is obtained as Let us now consider the relation between the complementary traction, the particular traction, and the true (nominal) traction.The nominal traction is defined as   =     .With respect to the introduced displacement partition, one has 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: Final equations: Considering that {  } = {} − {  } and {  } = {} − {  } − {}, Equation ( 8) is put in the following form: 16) and ( 17), it becomes 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, {} = { 0 } + {  (∇)}.
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 {  (∇)}.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 ∆ and ∆ be, respectively, the displacement increment and traction increment between two consecutives loading steps.The incremental form of Equation ( 16) reads In relation (19),   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 {∆} = [  ]{} + [  ] {∆  } and {∆} = [  ]{} + [  ] {∆  } , respectively, where {} , {∆  } , and {∆  } 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 {∇} = [] {} .After some algebraic manipulations, Equation ( 19) is written in the obvious compact form: The new matrices and vectors in Equation ( 18) are defined as follows: 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:

Numerical Examples
The generalized multi-quadrics radial basis functions   () = ( 2 +  2 )  are adopted. = ‖  − ‖ is the Euclidian distance between the field point  and the collocation center   .The constants  and  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.

Material Models
Five material models are considered.Subsequently, the following invariants of the right Cauchy-Green tensor are used: )), and  3 =  ().Given a unit vector , one also has The simplest nonlinear elastic model is known as the Saint Venant-Kirchhoff solid.The second Piola-Kirchhoff stress tensor is expressed as where E = (G − 1) /2 is the Green-Lagrange strain tensor, and  and  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  is given by In Equation ( 22),  −1 is the inverse of  and  = ().
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 where the material parameters  10 and  01 are related to the shear modulus by  = 2 ( 01 +  10 ), and  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 In the relations above,  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 ,   can be identified.
When the axis of isotropy coincides with the -axis, the five engineering material parameters are   , ,   ,  1 , and   .It can be shown that .

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.

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.

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.Under such a uniform loading state, the displacement field is given by   = ( 1 − 1) ,   = ( 1 − 1) , and   = ( 3 − 1) .Following this, the non-zero components of the deformation gradient tensor, the Green strain tensor () and its inverse ( −1 ), and the Green-Lagrange tensor () are gathered in Table 1 below.Table 1.Unidirectional loading: non-zero components of the deformation gradient, Green strain, inverse Green strain, and Green-Lagrange tensors.
The specimen is loaded either by applying the traction at the top surface (  known) or by imposing the displacement of points of this surface ( 3 known).From the numerical results, one can assess the values of  3 or   and  1 .Their analytical values can be calculated using the relations below where  =   ,  =   2(1+) , and  are, respectively, Young's modulus, the shear modulus, and Poisson's ratio.In the transversely isotropic case, these data are complemented by   ,  1 , and   : Saint Venant-Kirchhoff Solid: Neo-Hookean Solid: Mooney-Rivlin Solid: In our calculations, we used  = 5 MPa,  = 0.2, 10 Transversely Isotropic Saint Venant-Kirchhoff Solid: Transversely Isotropic neo-Hookean Solid: 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   at the top surface of the cube with respect to the applied stretch  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  1 with respect to  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  = 10 −3 .For the transversely isotropic case, two sets of engineering material parameters are adopted.They are presented in Table 2  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  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  3 were practically analytical solutions.The computed values of λ 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.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  and  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.

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 ninenode 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 (, , ) ≡ (, , ).
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.

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   = ( 1 − 1) +  ;   = ( 2 − 1);   = ( 3 − 1)  1 ,  2 , and  3 are the stretches and  is the shear factor.When  1 =  2 =  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  = ( 0  0  0 0 0 0 0 ) 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   =  ;   = 0; and   = 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  = 0.25 are given in Table 5 below.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.

Constrained Tension of a Cubic Specimen
The cross-section of the considered sample is presented in Figure 5.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.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  3 , the smaller the value of stretch  1 , the greater the lateral contraction.This result shows that the proposed approach is valid in this case too.

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

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.

Figure 1 .
Figure 1.Sample of domain Ω with scattered collocation centres.Star point with its local influence domain   .

Figure 2 .
Figure 2. Geometry of the simply supported (triangles) unit cube under uniform load (arrows) with associated boundary discretization.

Figure 5 .
Figure 5. Cross-section of the constrained tension specimen.

Figure 8 .
Figure 8. Geometry of the bending specimen in the initial configuration.

Figure 9 .
Figure 9. Displacement of the center line of a Saint Venant-Kirchhoff material under bending.

Figure 10 .
Figure 10.Displacement of the center line of a Mooney-Rivlin material under bending.

Table 2 .
below.Engineering material parameters of the transversely isotropic solids.

Table 3 .
Comparison of analytical and numerical stretch  1 for various values of  3 in the cases of transversely isotropic cube under tension load.

Table 4 .
Comparison of analytical and numerical stretches  1 and  3 in the cases of compression of a cylindrical solid.

Table 5 .
Non-zero components of the Cauchy stress tensor under simple shear of a cubic specimen (shear factor of 0.25).