Refined Meshless Local Strong Form solution of Cauchy-Navier equation on an irregular domain

TL;DR

This paper introduces a meshless local strong form method for solving the Cauchy-Navier equation in linear elasticity on irregular domains, demonstrating flexibility in node placement and refinement techniques with promising accuracy and convergence.

Contribution

It presents a novel meshless local strong form approach that handles irregular domains and node placement flexibility for solving elasticity problems.

Findings

Accurate solutions for irregular domain problems.

Effective node placement and refinement strategies.

Comparable or improved convergence rates.

Abstract

This paper considers a numerical solution of a linear elasticity problem, namely the Cauchy-Navier equation, using a strong form method based on a local Weighted Least Squares (WLS) approximation. The main advantage of the employed numerical approach, also referred to as a Meshless Local Strong Form method, is its generality in terms of approximation setup and positions of computational nodes. In this paper, flexibility regarding the nodal position is demonstrated through two numerical examples, i.e. a drilled cantilever beam, where an irregular domain is treated with a relatively simple nodal positioning algorithm, and a Hertzian contact problem, where again, a relatively simple h-refinement algorithm is used to extensively refine discretization under the contact area. The results are presented in terms of accuracy and convergence rates, using different approximations and refinement setups, namely Gaussian and monomial based approximations, and a comparison of execution time for each block of the solution procedure.

Equations

$$\hat{u}(\vec{p})=\sum _{i=1}^m{\alpha }_i{b}_i(\vec{p})=\mathbf{b}(\vec{p}{)}^{\mathsf{T}}\mathbf{\alpha },$$

$$u_j:=u({\vec{p}}_j)=\mathbf{b}({\vec{p}}_j{)}^{\mathsf{T}}\mathbf{\alpha },$$

$$\mathbf{u}=\mathbf{B}\mathbf{\alpha },$$

$$R^2=\sum _j^n{w({\vec{p}}_j)(u({\vec{p}}_j)-\hat{u}({\vec{p}}_j))}^2=(\mathbf{B\alpha }-\mathbf{u}{)}^{\mathsf{T}}{\mathbf{W}}^2(\mathbf{B\alpha }-\mathbf{u}),$$

$$w(\vec{p})=\exp \left(-{\left(\frac{\Vert {\vec{p}}_0-\vec{p}\Vert }{\sigma p_{\text{min}}}\right)}^2\right),$$

$${\mathbf{B}}^{\mathsf{T}}{\mathbf{W}}^2\mathbf{B}\mathbf{\alpha }={\mathbf{B}}^{\mathsf{T}}{\mathbf{W}}^2\mathbf{u}.$$

$$\mathbf{\alpha }=(\mathbf{W}\mathbf{B}{)}^{+}\mathbf{W}\mathbf{u},$$

$$\hat{u}(\vec{p})=\mathbf{b}(\vec{p}{)}^{\mathsf{T}}(\mathbf{W}\mathbf{B}{)}^{+}\mathbf{W}\mathbf{u}=\mathbf{\chi }(\vec{p})\mathbf{u},$$

$$(\mathcal{L}\hat{u})(\vec{p})=(\mathcal{L}\mathbf{\chi })(\vec{p})\mathbf{u}.$$

$${\mathbf{\chi }}^{\partial x}(\vec{p})$$

$${\mathbf{\chi }}^{\partial y}(\vec{p})$$

$${\mathbf{\chi }}^{\partial x\partial x}(\vec{p})$$

$${\mathbf{\chi }}^{\partial x\partial y}(\vec{p})$$

$${\mathbf{\chi }}^{\partial y\partial y}(\vec{p})$$

$$(\lambda +\mu )\nabla (\nabla \cdot \vec{u})+\mu {\nabla }^2\vec{u}=0,$$

$$(\lambda +\mu )\frac{\partial }{\partial x}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)+\mu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$$

$$(\lambda +\mu )\frac{\partial }{\partial y}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)+\mu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)$$

$$\sigma =\left[\begin{array}{cc}{\sigma }_{xx}& {\sigma }_{xy}\\ {\sigma }_{xy}& {\sigma }_{yy}\end{array}\right]$$

$${t_0}_1$$

$${t_0}_2$$

$$\left[\begin{array}{cc}U1& V1\\ U2& V2\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right]=\left[\begin{array}{c}{\mathbf{b}}_1\\ {\mathbf{b}}_2\end{array}\right],$$

$$\begin{array}{rl}U{1}_{i,\mathcal{N}(i{)}_j}& ={\left[(\lambda +2\mu ){\mathbf{\chi }}^{\partial x\partial x}({\vec{p}}_i)+\mu {\mathbf{\chi }}^{\partial y\partial y}({\vec{p}}_i)\right]}_j\\ V{1}_{i,\mathcal{N}(i{)}_j}& ={\left[(\lambda +\mu ){\mathbf{\chi }}^{\partial x\partial y}({\vec{p}}_i)\right]}_j\\ {{\mathbf{b}}_1}_i& =0\end{array}\},$$

$$\begin{array}{rl}U{2}_{i,\mathcal{N}(i{)}_j}& ={\left[(\lambda +\mu ){\mathbf{\chi }}^{\partial x\partial y}({\vec{p}}_i)\right]}_j\\ V{2}_{i,\mathcal{N}(i{)}_j}& ={\left[\mu {\mathbf{\chi }}^{\partial x\partial x}({\vec{p}}_i)+(\lambda +2\mu ){\mathbf{\chi }}^{\partial y\partial y}({\vec{p}}_i)\right]}_j\\ {{\mathbf{b}}_2}_i& =0\end{array}\},$$

$$\begin{array}{rl}U{1}_{i,\mathcal{N}(i{)}_j}& ={\left[\mu n_2{\mathbf{\chi }}^{\partial y}({\vec{p}}_i)+(2\mu +\lambda )n_1{\mathbf{\chi }}^{\partial x}({\vec{p}}_i)\right]}_j\\ V{1}_{i,\mathcal{N}(i{)}_j}& ={\left[\lambda n_1{\mathbf{\chi }}^{\partial y}({\vec{p}}_i)+\mu n_2{\mathbf{\chi }}^{\partial x}({\vec{p}}_i)\right]}_j\\ {{\mathbf{b}}_1}_i& =t_0({\vec{p}}_i{)}_1\end{array}\},$$

$$\begin{array}{rl}U{2}_{i,\mathcal{N}(i{)}_j}& ={\left[\mu n_1{\mathbf{\chi }}^{\partial y}({\vec{p}}_i)+\lambda n_2{\mathbf{\chi }}^{\partial x}({\vec{p}}_i)\right]}_j\\ V{2}_{i,\mathcal{N}(i{)}_j}& ={\left[\mu n_1{\mathbf{\chi }}^{\partial x}({\vec{p}}_i)+(2\mu +\lambda )n_2{\mathbf{\chi }}^{\partial y}({\vec{p}}_i)\right]}_j\\ {{\mathbf{b}}_2}_i& =t_0({\vec{p}}_i{)}_2\end{array}\},$$

$$\begin{array}{rl}U{1}_{i,i}& =1\\ {{\mathbf{b}}_1}_i& =u_0({\vec{p}}_i{)}_1\end{array}$$

$$\delta \vec{p}(\vec{p})=-{\sigma }_k\sum _{i=1}^{N_S}\nabla w(\vec{p}-\overrightarrow{p_i}),$$

$$\vec{p}\leftarrow \vec{p}+\delta \vec{p}(\vec{p}).$$

$${\vec{p}}_j^{new}=\frac{\vec{p}+{\vec{p}}_j}{2},$$

$${\sigma }_{xx}$$

$$u$$

$$v$$

$$u(L,y)$$

$$v(L,y)$$

Full text

††journal: Engineering Analysis with Boundary Elements

1 Introduction

Linear elasticity problems, governed by the Cauchy-Navier equation, are typically addressed in their weak form with the Finite Elements Method (FEM) [1]. However, the problem has also been addressed in its strong form in the past, e.g. component-wise iterative solution with the Finite Differences Method (FDM) [2] and with the Finite Volumes Method (FVM) [3]. Besides mesh based methods, meshless methods have also been employed for solving solid mechanics problems in strong and weak form [4, 5]. The conceptual difference between meshless methods and mesh based methods is in the treatment of relations between nodes. In the mesh based methods the nodes need to be structured into polygons (mesh) that cover the whole computational domain, while on the other hand, meshless methods fully define relations between nodes through the relative inter nodal positions [6], with an immediate consequence of greater generality of the meshless methods.

Strong form meshless methods can be understood as generalizations of FDM, where instead of predetermined interpolation over a local support, a more general approach with variable support and basis is used to evaluate partial differential operators [7], e.g. collocation using Radial Basis Functions [5, 8] or approximation with monomial basis  [9]. There are many other methods with more or less similar methodology introducing new variants of the strong form meshless principle [10]. On the other hand, weak form meshless methods are generalizations of FEM. Probably the most known method among weak form meshless methods is the Meshless Local Petrov Galerkin Method (MLPG) [11], where for each integration point a local support is used to evaluate field values and weight functions of a Moving Least Squares (MLS) approximation are used as test functions. In last few decades there have been many variants of MLPG introduced to mitigate numerical instabilities and to improve accuracy and convergence rate, etc. [10].

In general, recent developments in meshless community are vivid, ranging from analyses of computer execution on different platforms [6, 12], reducing computational cost by introducing a piecewise approximation [13] to implementation of more complex multi-phase flow [14], and many more.

This paper extends the spectra of published papers with a generalized formulation of a local strong form meshless method, termed Meshless Local Strong Form Method (MLSM) enriched with h-refinement [15] and ability to discretize arbitrary domains [7].

The introduced meshless approach is demonstrated on a solution of a benchmark cantilever beam case [16] and a Hertzian contact problem [17]. The results are presented in terms of displacement and stress plots, comparison against closed form solutions, convergence analyses, and execution time analyses.

The goal of this paper is to demonstrate generality of MLSM that is driven by the fact that all the building blocks of the method depend only on the relative positions between the computational nodes. This is a very useful feature, especially when dealing with problems in multidimensional spaces, complex geometries, and moving boundaries. This feature can be also exploited to write elegant generic code [18].

The rest of the paper is organized as follows: in section 2 the MLSM principle is explained, in section 3 the governing problem is introduced, section 4 is focused on solution procedure, section 5 focuses on discussing the results, and finally, the paper offers some conclusions and guidelines for future work in the last section.

2 MLSM formulation

The core of the spatial discretization used in this paper is a local approximation of a considered field over the overlapping local support domains, i.e. in each node we use approximation over a small local subset of neighbouring $n$ nodes. The trial function $\hat{u}$ is thus introduced as

[TABLE]

with $m$, $\boldsymbol{\alpha}$, $\boldsymbol{b}$ and $\vec{p}$ standing for the number of basis functions, approximation coefficients, basis functions and the position vector, respectively. In cases when the number of basis functions and the number of nodes in the support domain are the same, $n=m$, the determination of coefficients $\boldsymbol{\alpha}$ simplifies to solving a system of $n$ linear equations, resulting from evaluating equation (1) in each support node and setting it equal to a true value $u(\vec{p}_{j})$, for $j$ from 1 to $n$:

[TABLE]

where $\vec{p}_{j}$ are positions of support nodes and $u_{j}$ is the actual value of considered field in the support node $\vec{p}_{j}$. The above system can be written in matrix form as

[TABLE]

where $\boldsymbol{B}$ stands for coefficient matrix with elements $B_{ji}=b_{i}(\vec{p}_{j})$. The most known method that uses such an approach is the Local Radial Basis Function Collocation Method (LRBFCM) that has been recently used in various problems [5, 8].

In cases when the number of support nodes is higher than the number of basis functions ($n>m$) a WLS approximation is chosen as a solution of equation (3), which becomes an overdetermined problem. An example of this approach is DAM [9] that was originally formulated to solve fluid flow in porous media. DAM uses six monomials for basis and nine noded support domains to evaluate first and second derivatives of physical fields required to solve the problem at hand, namely the Navier Stokes equation. To determine the approximation coefficients $\boldsymbol{\alpha}$, a norm

[TABLE]

is minimized, where $\boldsymbol{W}$ is a diagonal matrix with elements $W_{jj}=\sqrt{w(\vec{p}_{j})}$ with

[TABLE]

where $\sigma$ stands for weight shape parameter, $\vec{p}_{0}$ for centre of support domain and $p_{\text{min}}$ for the distance to the nearest support domain node. There are different computational approaches to minimizing (4). The most intuitive and also computationally effective is to simply compute the gradient of $R^{2}$ with respect to $\boldsymbol{\alpha}$ and setting it to zero, resulting in a positive definite system

[TABLE]

The problem of this approach is bad conditioning, as the condition number of $\boldsymbol{B}^{\mathsf{T}}\boldsymbol{W}^{2}\boldsymbol{B}$ is the square of the condition number of $\boldsymbol{WB}$, unnecessarily increasing numerical instability. A more stable and more expensive approach is QR decomposition. An even more stable approach is SVD decomposition, which is of course even more expensive. Nevertheless, the solution of equation (6) can be written generally in matrix form as

[TABLE]

where $\boldsymbol{A}^{+}$ stands for a Moore–Penrose pseudo inverse of matrix $\boldsymbol{A}$. By explicitly inserting equation (7) for $\boldsymbol{\alpha}$ into (1), the equation

[TABLE]

is obtained, where $\boldsymbol{\chi}=\boldsymbol{b}(\vec{p})^{\mathsf{T}}(\boldsymbol{W}\boldsymbol{B})^{+}\boldsymbol{W}$ is called a shape function. Now, we can apply a partial differential operator $\mathcal{L}$ to the trial function, and get

[TABLE]

In this paper we deal with a Cauchy-Navier equation and therefore following shape functions are needed, expressed explicitly as

[TABLE]

The shape functions depend only on the numerical setup, namely nodal distribution, shape parameter, basis and support selection, and can as such be precomputed for a specific computation.

3 Governing problem

The goal in this paper is to numerically determine the stress and displacement distributions in a solid body subjected to the applied external force. To obtain a displacement vector field $\vec{u}$ throughout the domain, a Cauchy-Navier equation is solved, which can expressed concisely in vector form as

[TABLE]

where $\mu$ and $\lambda$ stand for Lamé constants. In two dimensions we express $\vec{u}=(u,v)$ and the equation reads

[TABLE]

Two types of boundary conditions are commonly used when solving these types of problems, namely essential boundary conditions and traction (also called natural) boundary conditions. Essential boundary conditions specify displacements on some portion of the boundary of the domain, i.e. $\vec{u}=\vec{u}_{0}$, while traction boundary conditions specify surface traction $\sigma\vec{n}=\vec{t}_{0}$, where $\vec{n}$ is an outside unit normal to the boundary of the domain and

[TABLE]

is the stress tensor. In terms of displacement vector $\vec{u}$ the traction boundary conditions read

[TABLE]

where ${t_{0}}_{i}$ and $n_{i}$ denote the Cartesian components of $\vec{t}_{0}$ and $\vec{n}$.

4 Solution procedure

4.1 Discretization of the problem

The elliptic boundary value problem at hand is discretized into a linear system of $2N$ algebraic equations by approximating the differential operations using MLSM, as described in section 2. A block system of linear equations for two vectors $\boldsymbol{u}$ and $\boldsymbol{v}$ of unknowns representing values $u(\vec{p_{i}})$ and $v(\vec{p_{i}})$, respectively, is constructed. This system is a discrete analogy of PDE (15) and can symbolically be represented as

[TABLE]

where $\boldsymbol{u}$ and $\boldsymbol{v}$ stand for unknown displacements, $\boldsymbol{b_{1}}$ and $\boldsymbol{b_{2}}$ for values of boundary conditions and blocks $U1,V1,U2,V2$ contain precomputed shape functions (10–14). With $\mathcal{N}(i)$ standing for a list of indices of the chosen $n$ neighbours of a point $\vec{p}_{i}$, as introduced in the beginning of section 2, we can, for all indices $i$ of internal nodes, express

[TABLE]

[TABLE]

for each $j=1,\ldots,n$. Note that equation (22) represents direct discrete analogue of (16) and, likewise, (23) of (17).

Similarly, for all indices $i$ of boundary nodes with traction boundary conditions we express

[TABLE]

[TABLE]

for each $j=1,\ldots,n$, where $n_{i}$ are the Cartesian components of the outside unit normal to the boundary in node $\vec{p}_{i}$. Again, equation (24) is a direct analogue of (19) and (25) of (20). And finally, for indices $i$ of nodes with essential boundary condition, we express

[TABLE]

System (21) is sparse with nonzero ratio of less then $2n/N$. An example of the matrix of this system for the cantilever beam problem described in section 5.1 is shown in Figure 1, where the block structure and different patterns for boundary and internal nodes are clearly visible.

4.2 Positioning of nodes in a complex domain

Meshless methods are advertised as the methods that do not require any topological relations among nodes. That implies that even randomly distributed nodes could be used [19]. However, it is well-known that with regularly distributed nodes one achieves much better results in terms of accuracy and stability [20]. This has also been recently reported for MLSM in a solution of a Navier-Stokes problem [7]. The reason behind the sensitivity regarding the distribution of nodes lies in the generation of shape functions. To construct a stable method well balanced support domains are needed, i.e. the nodes in support domain need to be distributed evenly enough [7]. This condition is obviously fulfilled in regular nodal distributions, but when working with more interesting geometries, the positioning of nodes requires additional treatment. In literature one can find several algorithms for distributing the nodes within the domain of different shapes [21, 22]. In this paper we will use an extremely simple algorithm, introduced in [7] to minimize the variations in distances between nodes in the support domain. The basic idea is to “relax” the nodes based on a potential between them. Since a Gaussian function is a suitable potential and already used as weight in the shape functions, the nodes are translated simply as

[TABLE]

where $\delta\vec{p}$, $\vec{p}_{i}$, $\sigma_{k}$ and $N_{S}$ stand for the translation step of the node, position of $i$-th support node, relaxation parameter and number of support nodes, respectively (Figure 2(a)). After offsets in all nodes are computed, the nodes are repositioned as

[TABLE]

Presented iterative process procedure begins by positioning the boundary nodes, which are considered as the definition of the domain and are kept fixed throughout the process.

4.3 h-refinement

Besides flexibility regarding the shape of the domain, nodal refinement is often mandatory to achieve desired accuracy in cases with pronounced differences in stress within the domain. A typical example of such situation is a contact problem [17]. To mitigate the error in areas with high stresses the h-refinement scheme is used., which has already been introduced into different meshless solutions [23, 24]. In context of local RBF approximation the h-refinement has been used in the solution of the Burger’s equation [15], where a quad-tree based algorithm has been used to add and remove child nodes symmetrically around the parent node in transient solution of Burgers’ equation. However, the algorithm presented in  [15] supported only regular nodal distribution. In this paper we generalize it also to irregular nodal distribution.

In each node to be refined, new nodes are added on the half distances between the node itself and its support nodes

[TABLE]

where index $j$ indicates $j$-th support node. When adding new nodes, checks are performed if the newly added node is too close to any of the existing nodes; in that case the node is not added. Moreover, if the refined node and support node are both boundary nodes, newly added node is positioned on the boundary (Figure 2(b)). This procedure can be repeated several times if an even more refined domain is desired. These subsequent refinements will be called levels of refinement and will be denoted as level $i$ for the refinement that resulted from $i$ applications of the described algorithm.

The described algorithm follows the concept of meshless methods and as such does not require any special topological relations between nodes to refine a certain part of the computation domain. It is also flexible regarding the dimensionality of the domain, i.e. there is no difference in implementation of 2D or 3D variant of the algorithm.

An example on a non-trivial refinement is demonstrated in Fig. 3, where a domain with a hole is considered. The vicinity of the hole is four times refined and then, to mitigate possible irregularities during refinement, relaxed.

4.4 Asymptotic complexity of MLSM

The asymptotic complexity analysis begins with an assumption that evaluations of basis functions, weights, linear operators and boundary conditions take $O(1)$ time. For simple domain discretization, such as the uniform grid in a rectangle or random positioning, $O(N)$ time is required, where $N$ stands for number of computational nodes. To find the neighbours of each point, a tree based data structure such as kd-tree [6], taking $O(N\log N)$ time to construct and $O(n\log N)$ time to query $n$ closest nodes, is used. The relaxation of nodal positions (see section 4.2) with $I$ iterations costs additional $O(InN\log^{2}N)$ time. Re-finding the support nodes by rebuilding the tree and querying for support nodes once again, requires another $O((N+n)\log N)$ time. Calculation of the shape functions requires $N$ SVD decompositions, each taking $O(nm^{2})$ time, as well as some matrix and vector multiplication of lower complexity. Assembling the matrix takes $O(nN)$ time and assembling the right hand side takes $O(N)$ of time. Then, the system is solved using BiCGSTAB iterative algorithm. The final time complexity is thus $O(InN\log^{2}N+(N+n)\log N+m^{2}nN)+T$, where $T$ stands for the time spent by BiCGSTAB.

For comparison, the complexity of a well-known weak form Element Free Galerking method (EFG) [25] differs from MLSM in construction of the shape functions, whose computation requires $O(Nn_{q}m^{2}n)$ time using EFG method, with $n_{q}$ standing for the number of Gauss integration points per node. Additionally, the number of nonzero elements in the final system of EFG is of order $n_{q}$ times higher than that of MLSM, again increasing the complexity of EFG.

5 Numerical examples

5.1 Cantilever beam

First, the standard cantilever beam test is solved to assess accuracy and stability of the method. Consider an ideal thin cantilever beam of length $L$ and height $D$ covering the area $[0,L]\times[-D/2,D/2]$. Timoshenko beam theory offers a closed form solution for displacements and stresses in such a beam under plane stress conditions and a parabolic load on the left side. The solution is widely known and derived in e.g. [16], giving stresses in the beam as

[TABLE]

and displacements as

[TABLE]

where $I=\frac{1}{12}D^{3}$ is the moment of inertia around the horizontal axis, $E$ is Young’s modulus, $\nu$ is the Poisson’s ratio and $P$ is the total load force.

In the numerical solution, traction free boundary conditions are used on the top and bottom of the domain, essential boundary conditions given by (31) are used on the right and traction boundary conditions given by (30) on the left

[TABLE]

[TABLE]

The problem is solved using MLSM method with $n=9$ or $n=13$ support nodes and Gaussian weight with $\sigma=1$ (see (5)). Two sets of basis functions are considered, 9 monomials

[TABLE]

and 9 Gaussian basis functions (see (5) for definition) centred in support nodes. In the following discussions these two choices of basis functions will be referred to as M9 and G9, respectively.

System (21) is solved with BiCGSTAB iterative algorithm [26] with ILUT preconditioning [27]. Values of $L=30\,\mathrm{m}$, $D=5\,\mathrm{m}$, $E=72.1\,\mathrm{GPa}$, $\nu=0.33$ and $P=1000\,\mathrm{N/m}$ were chosen as physical parameters of the problem.

The acquired numerical solution of the cantilever beam problem is shown in Figure 4.

The error of the numerical approximation of stresses and displacements is measured in relative discrete $L_{\infty}$ norm, using

[TABLE]

as error indicators, with $X$ representing the set of all nodes. Convergence with respect to the number of computational nodes is shown in Figure 5. The numerical approximations converge towards the correct solution in stress $(e_{\infty}(\sigma))$ norm as well, with approximately the same convergence rate.

It can be seen that monomials converge very regularly with order 1 as expected, while Gaussian functions exhibit slightly worse convergence. Such behaviour has already been reported in solution of diffusion equation, where MLSM with Gaussian basis failed to obtain accurate solution with a high number of computational nodes. More details about the phenomenon and further reading can be found in [28].

The method was compared to the standard Element Free Galerkin (EFG) method [29]. The EFG method used circular domains of influence with radius $d_{I}$ equal to 3.5 times internodal distance, a cubic spline

[TABLE]

for a weight function, $n_{q}=4$ Gaussian points for approximation of line integrals and $n_{q}=16$ points for approximating area integrals. Lagrange multipliers were used to impose essential boundary conditions.

The performance of EFG with respect to the number of nodes is much better than MLSM. However, a more fair comparison would take into account also a higher complexity of the EFG. This can be achieved by comparing error with respect to the number of MLS evaluations, which is the most time consuming part of the solution procedure. In Figure 5 it is demonstrated that although EFG provides much better results in comparison to MLSM at a given number of nodes, its accuracy becomes comparable to MLSM, when compared per number of MLS evaluations.

To asses the stability of the method regarding the nodal distribution, the following analysis was performed. A regular distribution of points as used in the solution in Figure 4 was distorted by adding a random perturbation to each internal node. Its position is altered by

[TABLE]

where $\delta$ is the distance to the closest node, and measuring the accuracy of the solution with respect to $\sigma$, representing magnitude of the perturbation. An example of original and perturbed node distributions are shown in Figure 6.

Accuracy of the solution with respect to the perturbation magnitude is presented in Figure 7. It is demonstrated that using monomials as a basis with 9 support nodes results in an unstable setup. On the other hand monomials with 13 support nodes are much more stable and equally accurate, while using Gaussian basis with high shape parameter is the most unstable setup. To mitigate the stability issue, a lower shape parameter can be chosen, however, at the cost of accuracy. Regardless of the setup one can expect the solution to be stable at least up to $\sigma\approx 0.1$. Note that using more nodes in support domain can also increase stability. Refer to [7] for more details.

Time spent on each part of the solution procedure is shown in Figure 8. All measurements were performed on a laptop computer with an Intel(R) Core(TM) i7-4700MQ @2.40GHz CPU and with 16 GiB of DDR3 RAM. MLSM is implemented in C++ [18] and compiled using g++ 7.1.7 for Linux with -std=c++14 -O3 -DNDEBUG flags.

It can be seen that solving the system (21) makes up for more than 50 % of total time spent. Around 70 % of that time is spent on computing the preconditioner. The only other significant factor is computing the shape functions taking approximately 40 % of total time. Domain construction and matrix assembly take negligible amounts of time, matching the predictions made by complexity analysis in section 4.4.

To emphasize the generality of MLSM method, a “drilled” domain is considered in the next step. Arbitrarily positioned holes are added to the rectangular domain. The positioning algorithm described in section 4.2 and h-refinement algorithm described in section 4.3 are used to distribute the nodes inside the domain and refine the areas around the holes. The boundary conditions in this example are $\vec{u}=0$ on the right, traction free on the inside of the holes and on top and bottom and uniform load of $P/D$ on the left. The computed solution is shown in Figure 9 along with the ordinary cantilever beam example.

Both solutions are coloured using von Mises stress $\sigma_{v}$, computed for the plane stress case as

[TABLE]

To further illustrate the generality of the method, an even more deformed domain is considered (Figure 10).

5.2 Hertzian contact

Another more interesting case arises from basic theory of contact mechanics, called Hertzian contact theory [30]. Consider two cylinders with radii $R_{1}$ and $R_{2}$ and parallel axes pressed together by a force per unit length of magnitude $P$. The theory predicts they form a small contact surface of width $2b$, where

[TABLE]

and

[TABLE]

Elastic modulus and Poisson’s ratio for the first material are denoted with $E_{1}$ and $\nu_{1}$, and with $E_{2}$ and $\nu_{2}$ for the second material. The pressure distribution between the bodies along the contact surface is semi-elliptical, i.e. of the form

[TABLE]

A problem can be reduced to two dimensions using plane stress assumption. A special case of this problem is when $E_{1}=E_{2}$, $\nu_{1}=\nu_{2}$ and $R_{2}\to\infty$, describing a contact of a cylinder and a half plane. This is the second numerical example tackled in this paper. The setup is ideal for testing the refinement, since a pronounced difference in behaviour of numerical solution near the contact in comparison to the rest of the domain is expected.

A displacement field $\vec{u}$ satisfying (15) on $(-\infty,\infty)\times(-\infty,0)$ with boundary conditions

[TABLE]

is sought. Vector $\vec{t}$ represents traction force on the surface and $\vec{\jmath}=(0,1)$ the upwards direction. Analytical solution for internal stresses in the plane in general point $(x,y)$ is calculated using the method of complex potentials [31] and the stresses are given in terms of $m$ and $n$, defined as

[TABLE]

where $m=\sqrt{m^{2}}$ in $n=\operatorname{sgn}(x)\sqrt{n^{2}}$. The stresses are then expressed as

[TABLE]

Numerically the problem is solved by truncating the infinite domain to a rectangle $[-H,H]\times[-H,0]$ for large enough $H$ and setting the essential boundary conditions $\vec{u}=0$ everywhere but on the top boundary. The top boundary has a traction boundary condition with normal traction given by $p(x)$ and no tangential traction. An illustration of the problem domain along with the boundary conditions is given in Figure 11.

The described contact problem is solved numerically and the error is measured between calculated and given stresses in relative $L_{\infty}$ norm as before, using

[TABLE]

as an error indicator. Values $P=543\,\mathrm{N/m}$, $E_{1}=E_{2}=72.1\,\mathrm{GPa}$, $\nu_{1}=\nu_{2}=0.33$, $R_{1}=R=1\,\mathrm{m}$ were chosen for the physical parameters of the problem. These values yield contact half-width $b=0.13\,\mathrm{mm}$ and peak pressure $p_{0}=2.6\,\mathrm{MPa}$. A value of $H=10\,\mathrm{mm}$ for domain height is chosen, approximately 38 times greater than width of the contact surface. Convergence of the method is shown in Figure 12.

It is clear that the convergence of the method is very irregular and slow. This is to be expected, as $N=10^{6}$ means only approximately 30 nodes positioned within the contact surface, and that naturally leads to large changes as a change of a single node bears a relatively high influence. Another problem is that the boundary conditions are only continuous, and exhibit no higher regularity, not even Lipschitz continuity. The accuracy of the approximation may seem bad, but is in fact comparable to the cantilever beam case. Using the comparable value of $N=30\cdot 15=450$ nodes in the contact area $[-b,b]\times[-b,0]$ it can be seen from Figure 5 that the approximation using this number of nodes in cantilever beam case achieved very similar results.

The total error of the approximation is composed of two main parts, the truncation error due to the non-exact boundary conditions and the discretization error, due to solving a discrete problem instead of the continuous one. First, we analyse the total error in terms of domain height $H$. A graph showing the total error with respect to domain height $H$ is shown in Figure 13.

The total error decreases as domain height increases, regardless of the discretization density used. However, as soon as truncation error becomes lower than discretization error, increasing the height further yields little to no gain in total error. The higher the discretization density is, the later this happens. When convergence of a method stops or significantly decreases in order, an error limit imposed by the truncation error was reached.

It soon becomes impossible to uniformly increase discretization density due to limited resources, and the immediate solution is to refine the discretization in the contact area with the h-refinement algorithm introduced in section 4.3. A domain of height $H=1\,\mathrm{m}\approx 75\,000b$ is chosen. Primary refinement is done in rectangle areas of the form

[TABLE]

and secondary refinement around points $\pm b$ on the surface is done in rectangle areas

[TABLE]

The refined domain as described above is shown in Figure 14. This domain was used to solve the considered contact problem. Different levels of secondary refinement were tested to prove that refinement helps with accuracy. Convergence of the method on the refined domain is shown in Figure 15.

Comparing Figure 15 to Figure 12, it can be seen that refinement greatly improves the accuracy of the method. Using $N=10^{6}$ nodes without refinement yields worse results than $N=10^{4}$ nodes with only primary refinement. Each additional level of secondary refinement helps to decrease the error even further while keeping the same order of convergence. A solution of the problem on the final mesh is shown in Figure 16.

6 Conclusions

A MLSM solution of a linear elasticity problem on regular and irregular domains with a refined nodal distribution of two different numerical examples is presented in this paper. The method is analysed in terms of accuracy by comparison against available closed form solutions and by comparison against weak form EFG method. The convergence of the method is evaluated with respect to the number of computational nodes, selection of different basis functions, different refinement strategies and different boundary conditions. MLSM is also analysed from complexity point of view, first, theoretically, and then also experimentally by timing the computer execution time of all main blocks of the method. It is clearly demonstrated that the method is accurate and stable. Furthermore, it is demonstrated that nodal adaptivity is mandatory when solving contact problems in order to obtain accurate results and that the proposed MLSM method can handle extensive refinement with the smallest internodal distance being $2^{17}$ times smaller than the initial one. It is also demonstrated that proposed MLSM configuration can handle computations in complex domains.

In our opinion the presented meshless setup can be used, not only to solve academic cases with the sole goal to show excellent convergences, but also in more complex engineering problems. The C++ implementation of presented MLSM is freely available at [18].

In future work we will continue to develop a meshless solution of a contact problem with a final goal to simulate a crack propagation due to the fretting fatigue [17] in a general 3D domain with added p-adaptivity to treat singularities near the crack tip.

Acknowledgement

The authors would like to acknowledge the financial support of the Research Foundation Flanders (FWO), The Luxembourg National Research Fund (FNR) and Slovenian Research Agency (ARRS) in the framework of the FWO Lead Agency project: G018916N Multi-analysis of fretting fatigue using physical and virtual experiments and the ARRS research core funding No. P2-0095.