1 Introduction

In the past three decades, there have been lots of works using the weak-form based methods to solve problems involving the geometric nonlinearity. Traditionally, the finite element method (FEM) is known to be the most popular one^[1-3], but numerical issues, such as the mesh distortion and reconstruction, have made it difficult to be applied to large deformation problems, which remains to be overcome. As a consequence, the meshfree methods with the point discretization have attracted great attention in recent years. The methods used are as follows: the smooth particle hydrodynamics (SPH)^[4], the diffuse element method (DEM)^[5], the element free Galerkin method (EFG)^[6-8], the reproducing kernel particle method (RKPM)^[9-16], and the $h$-$p$ clouds^[17]. In the aforementioned weak-form formulation, the numerical integration is needed for establishing the equations in the Galerkin meshfree methods. To account for Gauss quadrature's added complexity and direct nodal integration's instability, Chen et al.^[18] first proposed the stabilized conforming nodal integration (SCNI) to achieve the Galerkin linear exactness by using the strain smoothing stabilization for the nodal integration. Wang and Sun^[19] extended the SCNI for geometrically nonlinear analysis of shear deformable plates by a Lagrangian curvature smoothing within a nodal representative domain in the initial configuration. Chen et al.^[20] further formulated an arbitrary order variationally consistent integration by letting the test functions be variationally consistent with the integration scheme up to a desired order. For the buckling analysis of thin plates, Wang and Peng^[21] introduced the Hermite reproducing kernel meshfree approximation with the strain smoothing technique consistently utilized to construct the smoothed rotation and curvature fields. Very recently, Wang and Wu^[22] proposed a nesting sub-domain gradient smoothing integration algorithm, in which the quadratic exactness and efficiency are established through optimal two-level nesting sub-domains. Even with improvements in the numerical integration, the Dirichlet boundary conditions cannot be imposed directly, however, they can be done with the aid of techniques such as the Lagrange multiplier method^[6], the transformation method^[12], or the penalty method^{[14, 15]}.

Moreover, the meshfree methods formulated in the strong form do not require the domain integration, while enforcing the boundary conditions directly, which makes the solution procedure computationally efficient. The commonly used strong form collocation methods include the radial basis collocation method (RBCM)^[23] and reproducing kernel collocation method (RKCM)^[24-25], both of which have been successfully used to solve Poisson's problems and elasticity problems lately. On the subject of the RBCM, the collocation method with the radial basis function (RBF) approximation offers the exponential convergence rate, but leads to ill-conditioned systems^{[23, 26]}. As the shape parameter of the RBF governs the linear dependency of the discrete equations, it is closely related to the condition number of the system^[27]. A weighted RBCM was proposed to enhance the solution accuracy and convergence rate by imposing a higher weight on the Dirichlet boundary to balance the errors^[23].

To our knowledge, however, studies for the geometrically nonlinear analysis by the strong form collocation methods are very few. In this connection, there exist only iteration schemes in conjunction with the collocation methods for solving the nonlinear elliptic problems^[28]. Inspired by the aforementioned works and aiming at filling the existing gap, we propose a strong-form formulated incremental-iterative scheme using the total Lagrangian formulation with the aid of the RBCM for solving problems involving the geometric nonlinearity, with the nonlinear procedure presented in Ref. ^[29].

The arrangement of the paper is as follows. Section 2 gives the exposition of the problem of concern involving the geometric nonlinearity. In Section 3, the RBCM is introduced to discretize the nonlinear boundary value problem. In Section 4, the derivation of the discrete collocation equations is presented, and the algorithm for an incremental-iterative analysis is elaborated. In Section 5, the numerical experimentation is presented to demonstrate the reliability of the present approach. The concluding remarks are given in Section 6.

2 Geometric nonlinearity

For problems involving strong nonlinearities due to the change in the geometry during the deformation process, we consider a step-by-step analysis, in which all the state variables satisfying the equilibrium are known up to the last calculated time step $(t-1) $, and any further applied force and deformation will satisfy the equilibrium in the neighboring configuration at the time $t$. A body in the equilibrium at the current time step $t$ is described as follows:

(1)

where ${ }_0^t P_{ji} $ is the first Piola-Kirchhoff stress tensor, and ${ }_0^t B_i $ is the corresponding body force at the time $t$. As a reference, the subscript "0" denotes the initial undeformed configuration of the body. Let ${ }_0^t P_{ji} $ denote the first Piola-Kirchhoff stress tensor, which is related to the second Piola-Kirchhoff stress tensor ${ }_0^t S_{jk} $, i.e.,

(2)

where $F_{ik} $ is the deformation gradient. Due to lack of symmetry in the constitutive equation expressed by ${ }_0^t P_{ji} $, the second Piola-Kirchhoff stress tensor is introduced in the incremental analysis and expressed in the following incremental way:

(3)

where $\Delta({ }_0^t S_{ij})$ is the increment of the second Piola-Kirchhoff stress tensor, and ${ }^{t-1}{ }_0S_{ij} $ is the second Piola-Kirchhoff stress at the time $t-1$. The incremental constitutive law states

(4)

in which ${ }_0C_{ijkl} $ is the tangent modulus tensor, and $\Delta ({ }_0^t E_{kl})$ is the increment of the Green strain tensor.

By introducing Eqs. (2) -(4) to Eq. (1) , we can arrive at the following incremental equation:

(5)

where ${ }_0D_{ijkl} $ is the material response tensor, ${ }_0^t b_i $ is the corresponding body force measured in the undeformed configuration, and $F_{ij} { }_0^t B_j ={ }_0^t b_i $. Letting $\Delta({ }^{t-1}{ }_0e_{kl})$ and $\Delta({ }^{t-1}{ }_0 \eta _{kl})$ denote the linear and nonlinear parts of the increment of the Green strain tensor $\Delta({ }^{t-1}{ }_0E_{kl})$, respectively, we have

(6)

in which $\Delta({ }^{t-1}{ }_0e_{kl}) =(u_{k, l} +u_{l, k} )/2$ and $\Delta ({ }^{t-1}{ }_0 \eta _{kl}) ={u_{k, l} u_{l, k} }/2$. Substituting Eq. (6) into Eq. (5) leads to

(7)

To account for the geometric nonlinearity of the body, the geometric response tensor ${ }_0T_{ijkl} $ is introduced herein. Further, assuming that the strains are small, an approximate solution to Eq. (7) can be sought for the incremental step. Namely, the above equation can be approximated as

(8)

where the geometric response tensor ${ }_0T_{ijkl} $ is expressed as ${ }_0T_{ijkl} =\delta _{ik}\; { }^{t-1}{ }_0 S_{jl} $, and the material response tensor ${ }_0D_{ijkl} $ is defined as ${ }_0D_{ijkl} =F_{ip} F_{kq} { }_0C_{pjql} $. For isotropic Kirchhoff materials, the tangent modulus tensor ${ }_0C_{ijkl} $ is given by ${ }_0C_{ijkl} =\lambda \delta _{ij} \delta _{kl} +\mu ( {\delta _{ik} \delta _{jl} +\delta _{il} \delta _{jk} })$, where $\lambda $ and $\mu $ are Lame's constants.

Without loss of generality, the two terms on the right side of Eq. (8) can be interpreted as the external forces applied to the body at the time $t$ and those at the time $t-1$, respectively, and can be denoted by

(9)

Using the expressions in Eq. (9) , Eq. (8) can be rewritten as

(10)

Correspondingly, the boundary conditions for the nonlinear incremental problem are given by

(11)

where $n_j $ is the surface normal, $\Delta({ }_0^t h_i)$ is the prescribed traction on ${ }_0\Gamma _h $, $\Delta({ }_0^t u_i)$ is the incremental displacement, and $\Delta ({ }_0^t g_i) $ is the prescribed displacement on ${ }_0\Gamma _g $.

The total displacement ${ }_0^t u_i $ of the body can be obtained by accumulation of all the displacement increments generated in each previous step, i.e.,

(12)

As of now, the equation of equilibrium for the body has been formulated in an incremental form along with the boundary conditions.

3 Discretization of nonlinear boundary value problem by RBCM 3.1 RBF approximation

The collocation method with the RBF approximation is proposed to solve the nonlinear problems with large deformations. For multi-dimensional problems, the displacement and the increment of displacement can be respectively approximated at $N_{\rm s}$ source points as follows:

(13)

(14)

where ${ }^{t-1}{ }_0 u_i $ is the displacement evaluated at the time $t-1$, and $\Delta({ }_0^t u_i)$ is the increment of displacement evaluated at the current time step $t$. ${ }^{t-1}{ }_0 a_{iI} $ is the corresponding generalized coefficient evaluated at the time $t-1$, and $\Delta({ }_0^t a_{iI})$ is the corresponding increment of the generalized coefficient evaluated at the current time step $t$. The term $g_I ({ x})$ consists of RBFs. In this study, the following reciprocal multiquadrics (MQ) RBF is adopted:

(15)

where $c$ is the shape parameter of the RBF, $r_I $ denotes the radial distance between the collocation point ${x}=(x, y)$ and source point $ {x} _I =(x_I, y_I)$ of the RBF in two dimensions, which is given by $r_I =\sqrt {( {x-x_I } )^2+( {y-y_I } )^2} $.

3.2 Strong form discretization

For clarity, the left subscript "0" is omitted from now on. With reference to Eqs. (10) and (11) , the general form of a nonlinear boundary value problem in an incremental form is given by

(16)

where ${ }^{t-1}{L}^{\rm M}$ and ${ }^{t-1}{L}^{\rm G}$ are the differential operators in the domain $\Omega $ describing the material response and the geometric nonlinearity, respectively. ${ }^{t-1}{B}^h$ is the differential operator on the Neumann boundary $\Gamma _h $. ${ }^{t-1}{B}^g$ is the operator on the Dirichlet boundary $\Gamma _g $. The source terms are respectively defined as follows: ${ }^t{f}-{ }^{t-1}{f}$ in $\Omega $, $\Delta { }^t{h}$ on $\Gamma _h $ and $\Delta { }^t{g}$ on $\Gamma _g $.

For a two-dimensional elasticity problem, the operators ${ }^{t-1}{L}^{\rm M}$, ${ }^{t-1}{L}^{\rm G}$, ${ }^{t-1}{B}^h$, ${ }^{t-1}{B}^g$, and vectors ${ }^t{f}$, ${ }^{t-1}{f}$, $\Delta { }^t{h}$, $\Delta { }^t{g}$ corresponding to Eq. (16) , are derived as follows:

(17)

in which each entry is given by

(18)

(19)

(20)

Here, $F_{ij} $ is the component of the deformation gradient at the time $t-1$, and

(21)

with each entry given by

(22)

where $E_{ij} $ is the component of the Green strain tensor at the time $t-1$.

The boundary operators are given as follows:

(23)

in which each entry is listed as follows:

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

and

(32)

Here, ${I}$ is the identity matrix, and

(33)

(34)

(35)

in which the vector ${ }^{t-1}{f}$ is calculated as

(36)

with ${\bf\Psi}^{\rm T}$ denoting the shape function matrix containing the RBFs.

3.3 Nonlinear elasticity

In the discretization of the equation for the nonlinear problem in the strong form, we have introduced the following expressions for the operators of the displacement gradient $u_{i, j} $, the deformation gradient $F_{ij} $, and the Green strain tensor $E_{ij} $:

(37)

(38)

(39)

4 RBCM with iterative scheme 4.1 Derivation of discrete collocation equations

To derive the discrete equations for the iterative scheme concerning the geometric nonlinearity, we rewrite Eq. (16) in the following abridged incremental form:

(40)

where the source term in the differential equation, ${ }^t{f}-{ }^{t-1}{f}$, represents the unbalanced forces arising from the external forces acting at the time $t$ and the time $t-1$ as a result of geometric nonlinearity in the domain $\Omega $.

The approximation of $\Delta { }^t{u}$ is introduced as follows:

(41)

where ${\bf\Psi} ^{\rm T}$ is the matrix containing the RBFs, and $\Delta { }^t{a}$ is the corresponding vector containing generalized coefficients expressed as

(42)

and

(43)

Using the collocation method, the residuals are enforced to be zeros at collocation points. Letting ${P}$ be a set of $N_{\rm p} $ collocation points in $\Omega $, ${Q}$ be a set of $N_{\rm q} $ collocation points on $\Gamma _h $, and ${R}$ be a set of $N_{\rm r} $ collocation points on $\Gamma _g $, we have

(44)

Substitution of Eqs. (41) and (44) into Eq. (40) yields the following discrete collocation equation:

(45)

with

(46)

in which

(47)

(48)

4.2 Method of solution

In this study, the Newton-Raphson method is adopted to solve Eq. (40) , by which the iteration is performed at each incremental step under constant loads. The nonlinear equation for the $n$th iteration of the $t$th incremental step can be written as

(49)

The initial conditions for $n=1$ are given by

(50)

where the subscript "l" is used to indicate the last iteration in the $( {t-1} )$th incremental step. The total displacement for the $n$th iteration can be accumulated via

(51)

The stopping criterion for the iteration is chosen as

(52)

4.3 Algorithm for incremental-iterative nonlinear analysis

The step-by-step procedure for the geometric nonlinear analysis in this study is described in detail as follows:

(i) Select the total number of incremental steps that is allowed. Set the initial conditions as follows: ${ }^1{f}_{0} ={\rm {\bf 0}}$, ${ }^1{u}_{0} ={\rm {\bf 0}}$, ${ }^1{L}_{0}^{\rm G} ={\rm {\bf 0}}$, in which ${ }^1{L}_0^{\rm M} $ in Eqs. (17) -(20) and ${ }^1{B}_0^h $ in Eqs. (23) -(31) are reduced to the operators derived in the linear elasticity, given by

(53)

(54)

(ii) For the first iteration ($n=1) $ at each incremental step $t$,

(a) Calculate ${ }^t{L}_0^{\rm M} $, ${ }^t{L}_0^{\rm G} $, ${ }^t{f}_0 $, and ${ }^t{B}_0^h $ using Eq. (50) , in which ${ }^t{f}_0 $ is given by

(b) Calculate ${ }^t{f}_1 $ using Eq. (35) .

(c) Calculate ${ }^t{B}_0^g $ using Eq. (32) .

(iii) For the remaining iterations ($n\ge 2) $,

(a) Compute $u_{i, j} $, $F_{ij} $, and $E_{ij} $ using Eqs. (37) -(39) .

(b) Update ${ }^t{L}_{n-1}^{\rm M} $ and ${ }^t{L}_{n-1}^{\rm G} $ using Eqs. (17) -(20) and Eqs. (21) -(22) , respectively.

(c) Calculate ${ }^t{f}_n $ using Eq. (35) , and calculate ${ }^t{f}_{n-1} $ using

(d) Calculate ${ }^t{B}_{n-1}^h $ and ${ }^t{B}_{n-1}^g $ using Eqs. (23) -(31) and Eq. (32) , respectively.

(iv) Solve Eq. (49) for the generalized coefficient increment $\Delta { }^t{a}_n $. For the current iterative step, compute the displacement increment $\Delta { }^t{u}_n $ and the total displacement ${ }^t{u}_n $ using Eq. (51) .

(v) Check whether the error norm of the displacement given in Eq. (52) is satisfied. If the condition is not satisfied, let $n=n+1$ and return to the step (iii). Otherwise, go to the step (vi).

(vi) If the total number of incremental steps is less than the preset value, let $t:=t+1$ and go to the step (ii). Otherwise, stop the procedure.

5 Numerical examples

In the following numerical study, the proposed framework using the RBCM with reciprocal MQ RBFs will be validated through the solution of two benchmark problems that contain the geometric nonlinearity.

Example 1 Tensile bar problem

Consider the two-dimensional panel subjected to tension along the axial direction as shown in Fig. 1. The governing equation for the problem is incrementally expressed as

(55)

with the boundary conditions as

(56)

The panel is made of linearly elastic material. For $\nu =0$, the analytical solution is given by^[30]

(57)

To solve this tensile bar problem, the corresponding operators in Eq. (45) are listed in detail herein.

For $t=1$ and $n=1$ in the incremental-iterative analysis, the operators in Eqs. (53) and (54) are given by

(58)

(59)

where $E$ and $A$ denote Young's modulus and the cross-sectional area of the bar, respectively, and $n_{\rm s} $ denotes the unit outward normal.

For the remaining iterations ($n\ge 2) $, the operators ${ }^t{L}_{n-1}^{\rm M} $ in Eq. (17) , ${ }^t{L}_{n-1}^{\rm G} $ in Eq. (21) , and ${ }^t{B}_{n-1}^h $ in Eq. (23) are reduced to

(60)

(61)

(62)

where $F_{ij} $ denotes the components of the deformation gradient, and $E_{ij} $ denotes the components of the Green strain tensor.

For this problem, the following material parameters are adopted: $L=5$, $D=1$, $E=1$, $A=1$, and $\sigma =0.5$. As depicted in Fig. 2, the domain is discretized uniformly by using the same number of source points and collocation points given by $N_{\rm s} =N_{\rm c} =5\times 21$. The shape parameter of the RBF is chosen to be $c=2.60$. By using 720 total incremental steps, there are less than 3 iterative steps needed for each incremental step. In Fig. 3, the numerical result obtained by the RBCM is compared with the analytical solution, from which we observe that the RBCM result agrees very well with the analytical one.

Example 2 Cantilever beam problem

Consider a two-dimensional cantilever beam under pure bending as shown in Fig. 4. The parameters for this problem are given as follows: $L=3~000$~mm, $D=300$~mm, $E=210$~GPa, $\nu =0.3$, and $\sigma =0.15$~GPa. The corresponding boundary conditions for the problem are given by

(63)

As illustrated in Fig. 5, the same numbers of uniformly distributed source points and collocation points, $N_{\rm s} =N_{\rm c} =5\times 41$, are used for this problem. The result is obtained by using the shape parameter $c=0.185$ and the total incremental step $t=132$ (to meet the data available in the literature) as shown in Fig. 6, where $u_{\rm A}$ and $v_{\rm A}$ denote the horizontal and vertical displacements, respectively. For each incremental step, there are less than 2 iterations performed, as a demonstration of the convergence and effectiveness of the nonlinear algorithm. Since there is no analytical solution available for the present problem, the numerical result obtained by the FEM using ANSYS, as provided in Ref. ^[16] for the use of comparison of different numerical methods, is readily adopted as the reference solution herein. From Fig. 6, it is observed that the RBCM result matches the finite element solution very well.

6 Concluding remarks

In this study, a strong-form framework is proposed for solving the boundary value problems involving the geometric nonlinearity, in which a detailed formulation of the incremental-iterative algorithm in solving geometrically nonlinear elasticity problems is provided. The total Lagrangian formulation is first adopted to describe the equilibrium of a body involving large deformation in an incremental form. To describe the effect of strong geometric nonlinearity of the body, in addition to the material response tensor, the geometric response tensor is introduced. Then, the approximation with the RBF is introduced to discretize the nonlinear boundary value problems such that the discrete system of the collocation equations can be established. To incrementally seek the approximate solution, the Newton-Raphson method is adopted in developing the incremental-iterative algorithm for the geometric nonlinear analysis. Compared with the weak-form Galerkin meshfree methods, the proposed strong-form framework has the following advantages. Since the strong form with direct collocation enforces the residuals to be zeros straightforwardly at the collocation points, obviously, no integration quadrature is needed in constructing the governing equation, and boundary conditions can be imposed directly. Furthermore, the point-based discretization facilitates geometrically nonlinear analysis during the incremental-iterative process with no requirement for mesh reconstruction. Consequently, the proposed strong-form algorithm provides a simple way for the large deformation analysis. The performance of the proposed strong-form framework is investigated through the two benchmark problems, and the numerical results demonstrate the effectiveness and capability of the present method for the geometrically nonlinear analysis.