Keywords

1 Introduction

In the recent years, many authors have shown the interest in developing mesh-free methods to overcome some difficulties. These methods provide the possibility to avoid mesh generation and connectivity between the nodes which require a considerable computation time. Nevertheless, not all mesh-free methods have this advantage. They are divided in two categories. The first one concerns the mesh-free methods based on weak formulation that still require a background mesh like the Diffuse Element Method (DEM) [1], Element Free Galerking (EFG) [2] and Point Interpolation Method (PIM) [3] etc. The second category of mesh-free methods is based on strong formulation. This category don’t need a background mesh and the numerical integration [4, 5].

The Weighted Least Squares (WLS) approximation is wildly used in data fitting [6]. Wallstedt et al. [7] took advantage of the concepts of weighted least squares surface estimation and implicit surface definition to more precisely define the region of integration for solid mechanic’s simulations that are performed within a PIC framework. Onate et al. [8] presented an approach termed generically the ‘Finite Point Method’ based on a Weighted Least Square Interpolation of point data and point collocation for evaluating the approximation integrals. Baeza et al. [9] have presented a technique based on the application of a variant of the Lagrange extrapolation through the computation of weights capable of detecting regions with discontinuities. Recently several authors developed a high order continuation using the Moving Least Squares approximation [10] as a discretization method under a strong formulation to solve various non linear problems [11,12,13,14,15,16,17,18,19]. Other authors used a high order continuation using Radial Point Interpolation method under a strong formulation [20].

In this paper, a new numerical approach is proposed. This approach combines a high order continuation [21] and the Weighted Least Squared (WLS) as a discretization method. Thanks to the development of the main variables into Taylor series, non-linear elastic problem is transformed into a succession of linear differential equations with the same tangent operator followed of a continuation technique is used to obtain the whole solution. In this work, the proposed approach has been referred to it as HOC-WLS. The performance of this proposed approach is illustrated on examples of non-linear elasticity problems. A comparison has been held between the proposed approach, the High Order Continuation coupled with Moving Least Squares (HOC-MLS) and the High Order Continuation coupled with Finite Element Method (HOC-FEM).

2 Description of Weighted Least Squares (WLS) and Moving Least Squares (MLS) Methods

In this section, we present in detail the description of the approximations by the Weighted Least Squares and by the Mobile Least Squares methods.

2.1 Weighted Least Squares Approximation

The Weighted Least Squares (WLS) approximation is highly used for data fitting [6]. The approximation \(u^{h}(x)\) of a field u(x) at a point x by WLS in an influence domain \(\varOmega \) is given by:

(1)

where \({<}p(x){>}={<}p_1(x),p_2(x),\cdots ,p_m(x){>}\) is a vector of monomial basis functions, with m is the number of monomials and \({}^t\{a\}={<}a_1,a_2,\cdots ,a_m{>}\) is a vector of the unknown constant coefficients. The most used vector of monomial basis functions is defined as follows:

(2)

To determine the vector of constant coefficients \(\{a\}\) in Eq. (1), we select n nodes in the local support domain for the approximation. We can solve the Eq. (1) for all nodes in the support domain using a Weighted Least Square (WLS) method by minimizing the weighted discrete \(L_{2}\) norm:

(3)

where \(w_{i}\) is the weight function associated to the node i and n is the number of nodes in the support domain. We note that weight functions plays an important role in constructing the WLS shape function. We should note that the weights used here are considered constant i.e they do not depend on x. They can be computed from any weight function with the bell shape [6].

For an arbitrary point x, the vector \(\{a\}\) is chosen to minimize the weighted residual. The stationary condition of the quadrature form J is given by:

(4)

After the calculation, the stationarity condition (4) can take the following form:

(5)

where the matrices [A] and [B] are defined by:

(6)

where \([P_m]\) is a matrix of order \(n\times n\) evaluated at all nodes of influence domain, [w] is the diagonal matrix constructed from the weight constants. The vector \(\{U\}\) collects all nodal unknowns. The resolution of Eq. (5) gives the following expression of the vector \(\{a\}\):

$$\begin{aligned} \{a\}=[A]^{-1}[B]\{U\} \end{aligned}$$
(7)

Using Eq. (7), the local approximation \(u^{h}(x)\) of Eq. (1) is written as follows:

(8)

where the vector of shape functions \({<}\phi (x){>}\) and its derivatives are defined as follows:

(9)

where the symbol \((\bullet )_{,i}\) denotes the derivative with respect to the ith coordinate and the symbol \((\bullet )_{,ij}\) denotes the derivative with respect to ith and to jth coordinates. It should be noted that the vector \(\{a\}\) is constant and the derivatives of the matrices [A] and [B] are null. This leads to a considerable reduction in computation time. Under a strong formulation, the shape function can be derived twice or more and we notice when the matrix [A] is differentiable and it’s derivatives increase the instability of the solution. If the matrix [A] is constant, the degree of instability of the solution decreases. This effect has been investigated in the numerical examples.

The shape functions so constructed do not have the Kronecker delta function property, which can cause difficulties in imposing the boundary conditions, [6]. In this case, the collocation methods are necessary [6] to overcome this difficulty. Note also that the WLS shape functions are compatible only in the influence domain rather than in the global domain. This is not a problem when the WLS shape functions are used in the local weak-form methods or collocation methods. But, care needs to be taken when it is used in a global weak-form formulation for which the Moving Least Squares method is chosen of a smart manner [6]. Fortunately in this paper we are using the WLS under a strong formulation (mesh-free collocation methods) and a local support domain.

2.2 Moving Least Squares Approximation

The Moving Least Squares (MLS) method is an efficient numerical method that is classed as a meshless approach that have a highly accurate approximation. Due to the method flexibility, numerous authors have used it to solve a large number of problems [6, 10]. For the implementation of this method, we use the same procedure as presented in Sect. 2.1. In this case, the approximation \(u^{h}(x)\) is defined as follows:

(10)

In the same way as before, the weighted residual is constructed and minimized with respect to the vector \(\{a(x)\}\) which depends on the point x in this case. In addition, in this case the weighted function depends on the variable x. In the same manner as in Sect. 2.1, the vector \(\{a(x)\}\) is the solution of the following system:

$$\begin{aligned}{}[A(x)]\{a(x)\}=[B(x)]\{U\} \end{aligned}$$
(11)

where the matrices [A(x)] and [B(x)] are defined by:

(12)

and \(\{U\}\) is the vector that collects the nodal unknowns of the influence domain. We should note that the matrix [A(x)] is not always invertible and depends on the nodal distribution inside the support domain. Assuming that [A(x)] is invertible, the vector \(\{a(x)\}\) is given by:

$$\begin{aligned} \{a(x)\}=[A(x)]^{-1}[B(x)]\{U\} \end{aligned}$$
(13)

Substituting the Eq. (13) in Eq. (10) leads to:

$$\begin{aligned} u^{h}(x)={<}\phi (x){>}\, \{U\} \end{aligned}$$
(14)

In this case, the vector of shape functions \({<}\phi (x){>}\) and its derivatives are defined as follows:

(15)

Properly constructed MLS shape functions are compatible and consistent with the order of polynomials included in the formulation. These shape functions do not have the Kronecker delta function property because they are not interpolation functions.

3 Strong Form Formulation of Two-Dimensional Nonlinear Elastic Problems

We consider a two dimensional elastic problem described by the equilibrium equations and boundary conditions. Let \(\varOmega \) be the domain occupied by a two dimensional solid, \(\partial \varOmega _{u}\) and \(\partial \varOmega _{f}\) are complementary portions of \(\partial \varOmega \) in which the Dirichlet and Newman conditions are applied. The equilibrium of this solid is governed by the following matrix equations:

(16)

where the matrices [div] and [N] are given by:

(17)

\(\{T\}\) is the vector containing all the components of the first Piola-Kirchhoff stress tensor, \(N_{x}\) and \(N_{y}\) are the components of the normal vector outward the boundary \(\partial \varOmega _{f}\), \(\lambda \) is a control parameter, \(\{u_{d}\}\) is the imposed displacement on \(\varOmega _{u}\) and \(\{f\}\) is the stress vector applied on \(\varOmega _{f}\). The relation between the first and second Piola-Kirchhoff stress tensors and the constitutive law are represented by the following equations:

(18)

The matrices [D], [H], \([A(\theta )]\), [III] and \([B(\theta )]\) are defined as follows [20]:

(19)

where u and v describe respectively the horizontal and vertical displacements, [D] is the elasticity matrix and \(\{\gamma \}\) is the strain vector. For a plane stress \(E^{\prime }=E\) and \(\nu ^{\prime }=\nu \) and for a plane strain \(E^{\prime }=\frac{E}{(1-\nu ^2)}\) and \(\nu ^{\prime }=\frac{\nu }{1-\nu ^{2}}\) with E is the Young modulus and \(\nu \) is the Poisson ratio.

Finally, we have obtained a nonlinear problem that requires a resolution numerical method. For this effect, we have proposed a high order continuation coupled with the Weighted Least Squares method HOC-WLS and with Moving Least Squares method HOC-MLS to solve this nonlinear problem.

4 Resolution Strategy

The solution of the problem governed by Eqs. (16) and (18) is sought under the Taylor series expansion truncated at order P with respect to the path parameter “r” [21, 22]:

(20)

where the terms with the subscript 0 represent the initial solutions of the starting point of each solution branch and the terms with the subscript i are the solutions at each order i. Injecting the development (20) in Eqs. (16) and (18), we have obtained a succession of linear problems which have the same tangent operator. This succession of linear problems is written as follows: Problem at order \(i=1\)

(21)

Problem at order \(2\le i\le P\)

(22)

where the terms with the superscript nl are found using the solutions at the previous orders. The approximation of the unknown \(\{u_i\}\) using WlS approximation is given by:

(23)

where n is the number of nodes inside the support domain, \([\phi (x)]\) is the matrix of the shape functions and \(\{U_{n}\}\) is the vector collecting the nodal unknown displacements in the support domain. Taking into consideration the approximation of the unknown in (23) and after substitution and assembly technique, we get the discrete problems of Eqs. (21) and (22) completed by the pilotage equation in the following compact form:

(24)

where \([K_{T}]\) is the tangent stiffness matrix evaluated at the initial solution (\(\{T_{0}\}\), \(\{S_{0}\}\), \(\{U_{0}\}\), \(\{\gamma _{0}\}\)), \(\{F^{nl}_{i}\}\) is a right hand side that depends on the previous orders and \(\{U_i\}\) is the vector of unknown nodal displacements at each order i. The Taylor series expansion of Eq. (20) has a validity radius “\(r_{max}\)” given by:

$$\begin{aligned} r_{max}=\left( \varepsilon \frac{\parallel \{U_{1}\}\parallel }{\parallel \{U_{P}\}\parallel }\right) ^{\frac{1}{P-1}} \end{aligned}$$
(25)

where \(\varepsilon \) is a given tolerance parameter.

5 Numerical Examples and Discussions

In this section, the robustness and efficiency of the proposed approach has been putting to the test. To compute the shape functions in each examples, the truncated Gauss-type distribution is used as a weight function.

(26)

the influence domain for all nodes is a circle with a varying radius \(r=h\times d_{max}\), \(q={d_{i}}/{r}\), \(d_{i}=\parallel x-x_{i} \parallel \), h is controlled by the user and \(d_{max}\) is the distance between two successive nodes (see Fig. 1). In the case of a nodes irregular distribution inside a varying influence domain, we measure \(d_{max}\) only for the points inside this domain.

The obtained solutions of each problem will be verified by computing the residual logarithmic norm \(log_{10}(\parallel Res\parallel )\) of problem (16) via inserting the obtained solution of Eq. (20).

Fig. 1.
figure 1

Influence domain of a point

Full size image

5.1 Two Dimensional Elastic Plate in Bending

To investigate the convergence of the two approaches, geometrical nonlinearity is considered in this example. A two dimensional plate is clamped on the left hand side and subjected to loading \(\lambda F\) on the right hand side, with \(F=1\) N (Fig. 2). The mechanical and geometrical characteristics of the plate are: length \(L=100\) mm, height \(l=10\) mm, Young’s modulus \(E=210\) GPa and Poison’s ratio \(\nu =0.3\). The comparison has been executed with same control parameters of the both approaches for various values of parameter h, the truncation order \(P=15\), the tolerance parameter \(\varepsilon =10^{-6}\), a shape parameter \(c=0.2\) and fixed nodal distribution 605 and for HOC-FEM approach as the reference solution, we take 825 T6 elements. In Table 1, we present the validity range “\(r_{max}\)” and the solution quality measured by “\(log_{10}(\parallel Res\parallel )\)” obtained by HOC-WLS approach for different numbers of monomial basis and parameter h. In this table, we remark that the proposed approach converges from \(h=2.3\) for the quadratic basis while for cubic and quartic basis it converges from \(h=2.8\). When this approach converges, the validity range is almost constant. In Table 2, we present the validity range “\(r_{max}\)” and the solution quality “\(log_{10}(||Res||)\)” obtained by HOC-MLS approach for different numbers of monomial basis and parameter h. In this table, we remark that the HOC-MLS approach converges for the values \(2.3\le h\le 3\) for the quadratic basis while for cubic basis it converges only for \(h=3\). The HOC-MLS approach diverges always for quartic basis. When this approach converges, the validity range is almost constant.

Fig. 2.
figure 2

Two dimensional elastic plate and nodal distribution

Full size image
Table 1. Effects of parameters of HOC-WLS on the validity range and the solution quality
Full size table
Table 2. Effects of parameters of HOC-MLS on the validity range and the solution quality
Full size table

In Fig. 3, we plot the load-displacement curve obtained by HOC-WLS, HOC-MLS and HOC-FEM algorithms. The horizontal displacement u and the vertical displacement v are plotted at the interest point. The HOC-WLS and HOC-MLS algorithms are in good agreement with the reference solution obtained by HOC-MEF algorithm. However, the HOC-WLSM algorithm shows a better solution quality than that of HOC-MLS algorithm.

Fig. 3.
figure 3

Load-displacement and residual-load curves at the registration point \((x=100,y=10)\) for a quadratic basis and \(h=2.5\)

Full size image

5.2 Two Dimensional Plate Under Traction Load

We consider a two dimensional plate with a length \(L=100\) mm and a height \(l=10\) mm. The structure is clamped at the left hand side \(x=0\) and subjected to a loading at the right hand side \(x=100\) (see Fig. 4a). The mechanical properties are: Young’s modulus \(E=210\) GPa and Poison’s ratio \(\nu =0.3\). We adopt the following numerical parameters: a truncation order \(P=15\), a tolerance parameter \(\varepsilon =10^{-8}\) and a shape parameter \(c=0.2\). The domain occupied by the plate is replaced by 583 nodes irregularly distributed (see Fig. 4b). For the reference solution, 3825 T6 elements are used.

Fig. 4.
figure 4

Two dimensional elastic plate under traction load and irregular nodal distribution

Full size image

The Fig. 5 represents the load-displacements curves at the interest point \((x=100, y=50)\) obtained the both mesh-free approaches and by the reference algorithm. The horizontal displacement u and the vertical displacement v are plotted at the interest point. Even when using an irregular nodal distribution, the load-displacement curve computed by the HOC-WLS is in a good agreement with the reference solution. Whereas the load-displacement curve obtained by HOC-MLS algorithm could not reach the same load as the proposed approach with the same number of steps. This shows that the MLS approximation does not give the same solution when the nodes are distributed irregularly.

Fig. 5.
figure 5

Load-displacement and residual-load curves at the registration point \((x=100,y=50)\)

Full size image

5.3 Two Dimensional Plate with Three Holes

In this numerical study, the effect of complex geometry on the proposed approach is tested. Two dimensional plate with three holes is considered which is clamped on the left hand side \(x=0\) and subjected to a load \(\lambda F\) on the top hand side \(y=l\) (see Fig. 6a). The mechanical and geometrical properties are: Young’s modulus \(E=210\) GPa, Poison’s ratio \(\nu =0.3\), Length \(L=100\) mm, height \(l=10\) mm and radius of each hole is \(a=3\) mm. We adopt the following numerical parameters: a truncation order \(P=15\), tolerance parameter \(\varepsilon =10^{-8}\) and a shape parameter \(c=0.2\). The domain occupied by the plate is replaced by 1030 nodes (see Fig. 6b). For the reference solution a 1784 T6 elements are used. In Fig. 7, the load-displacement curve is plotted at the interest point \((x=100, y=5)\).The horizontal displacement u and the vertical displacement v are plotted at the interest point. The obtained solution is compared to HOC-MLS approach and the reference solution. From Fig. 7a, the obtained solution is in a good agreement with the reference solution, while the one obtained by HOC-MLS diverges. In Fig. 7b, the proposed approach shows a good solution quality than the HOC-MLS approach.

Fig. 6.
figure 6

Two dimensional plate with three holes and nodal distribution

Full size image
Fig. 7.
figure 7

Load-displacement and residual-load curves at the interest point \((x=100, y=5)\)

Full size image

6 Conclusion

A new approach for non-linear elasticity is presented in this work. This approach is based on the WLS approximation and a High Order Continuation under a strong formulation (HOC-WLS). The HOC-WLS requires no mesh generation and it avoids numerical integration which is computationally expensive. The proposed approach is tested and verified in various numerical examples of 2D elasticity with different nodal distributions and geometries. These tests were done also by the HOC-MLS approach which is another mesh-free approach of the same category as the proposed one and HOC-FEM that represents a reference solution. These numerical tests are done to find the best numerical parameters for the proposed approach. Both mesh-free methods show a good quality results when the nodal distribution is regular. However, the HOC-WLS approach shows a better convergence and efficiency than HOC-MLS regarding irregular nodal distributions and complex geometries.