Notice: Undefined offset: 2690 in /var/www/scholarpedia.org/mediawiki/includes/parser/Parser.php on line 5961
Meshless methods for PDEs - Scholarpedia

# Meshless methods for PDEs

Post-publication activity

Curator: Darrell Pepper

Although many numerical and analytical schemes exist for solving engineering problems, the meshless method is a particularly attractive method that is receiving attention in the engineering and scientific modeling communities. Finite difference (FDM), finite volume (FVM), and finite element (FEM) methods have been historically used to model a wide variety of engineering problems in complex geometries that may require extensive meshing. The meshless method is simple, accurate, and requires no meshing.

The need to accurately simulate various physical processes in complex geometries is important, and has perplexed modelers utilizing conventional numerical schemes for many years. Today, advances in numerical schemes and enhanced hardware have lead to many commercial codes that can employ Herculean efforts to solve complex stress-strain, heat transfer, fluid flow, and other nearly intractable problems [1,2]. Recently, advances in the development and application of meshless techniques show they can be strong competitors to the more classical finite difference/volume and finite element approaches. Textbooks by Liu  and Fasshauer  discuss meshfree methods, implementation, algorithms, and coding issues for stress-strain problems; Liu  includes Mfree2D, an adaptive stress analysis software package available for free from the web, and Fasshauer  include MATLAB modules. Atluri and Shen  also produced a textbook that describes the meshless method in more detail, including much in-depth mathematical basis.

## History of Meshless Methods

Meshless methods are uniquely simple, yet provide solution accuracies for certain classes of equations that rival those of finite elements and boundary elements without requiring the need for mesh connectivity. Ease in programming, no domain or surface discretization, no numerical integration, and similar formulations for 2-D and 3-D make these methods very attractive.

There exist several types of meshless methods. The more common techniques include kernel methods, moving least square method, meshless Petrov- Galerkin, partition of unity methods, smooth-particle hydrodynamics, and radial basis functions. Each technique has particular traits and advantages for specific classes of problems. Generally the simplest and easiest to implement is the radial basis function approach. The examples illustrated here are based on radial basis functions (RBFs) and Kansa’s approach .

Over the last 10 years, development in using RBFs as a meshless method approach for approximating partial differential equations has accelerated. Kansa’s method, which is a domain-type meshless method, was developed by Kansa in 1990  by directly collocating RBFs, especially multiquadric approximations (MQ). The use of MQ was first developed by Hardy  in 1971 as an interpolation method for modeling the earth’s gravitational field. Franke  published a review paper on 2-D interpolation methods and ranked MQ as the best based on its accuracy, speed, storage requirements, and ease of implementation. It wasn’t long before researchers began to employ Kansa’s method to solve a number of engineering problems including the 1-D nonlinear Burgers equation, heat transfer, and free boundary problems (see Li et al, ).

More recently, the use of the Method of Fundamental Solutions (MFS) with the dual reciprocity method (DRM) and RBFs has appeared . The MFS method is also meshless, fast, and can be extended to multiple dimensions – and is especially attractive for problems involving complex geometry. The MFS is equivalent to the more widely known Boundary Element Method (BEM) (see [11-13]). Previous use of the MFS limited its application to homogeneous equations; however, by coupling the MFS with the DRM and RBFs has lead to the creation of an effective boundary-type meshless method for solving a wide range of partial differential equations. The MFS-DRM is slightly more computationally expensive than Kansa’s method, more difficult to implement, and requires the existence of fundamental solutions. However, it is more mathematically robust.

Meshless methods utilizing RBFs create mesh-free algorithms that are significantly simpler to employ than FDM/FVM, FEM, and BEM approaches, and truly eliminate the need for meshes requiring connectivity and optimization.

RBFs are the natural generalization of univariate polynomial splines to a multivariate setting. The main advantage of this type of approximation is that it works for arbitrary geometry with high dimensions and it does not require a mesh at all. A RBF is a function whose value depends only on the distance from some center point. Using distance functions, RBFs can be easily implemented to reconstruct a plane or surface using scattered data in 2-D, 3-D or higher dimensional spaces. Due to the uses of the distance functions, the RBFs can be easily implemented to reconstruct the surface using scattered data in 2D, 3D or higher dimensional spaces.

## The Kansa Approach

In Kansa’s method, it is assumed that a variable (assume temperature) can be expressed as an approximation in the form

$\tag{1} T(\textbf x)=\sum_{j=1}^{N}\phi_{j} (\textbf x)T_{j}$

where $${T_{j}}$$ are the unknown temperature values to be determined, and $$\phi_{j}(x) = \phi ||x-x_{j}||$$ is some form of RBF where $${\textbf x_i}, i=1,\ldots,N_I,$$ are interior points and $${\textbf x_i}, i = N_I + 1,\ldots,N,$$ are boundary points. Some popular choices of RBFs include linear $$(r)\ ,$$ cubic $$(r^3)\ ,$$ multiquadrics (MQ) $$((r^2+c^2)^{1/2})\ ,$$ polyharmonic splines ($$r^{2n+1}\log r$$ in 2-D, $$r^{2n+1}$$ in 3-D), Gaussian ($$\exp(-cr^2)$$). The theory of RBFs interpolation is discussed in Powell ; Fasshauer  lists many RBFs and their derivatives.

To illustrate the application of the meshless method with RBF, consider the 2-D Poisson’s Equation

$\tag{2} \begin{array}{l} {\nabla^{2} T=f(\textbf x),\quad \textbf x\in \Omega } \\ {T=g(\textbf x),\quad \quad \textbf x\in \Gamma } \end{array}$

where $$\textbf x=(x,y)\ .$$ Now approximate $$T(\textbf x)$$ assuming

$\tag{3} T(\textbf x)=\sum_{j=1}^{N}\phi (r_{j} )T_{j}^{}$

where $$r$$ is defined as

$\tag{4} r_{j} =\sqrt {(x-x_{j} )^{2} +(y-y_{j} )^{2} }$

To be more specific, we choose MQ as the basis function

$\tag{5} \phi (r_{j} )=\sqrt {r_{j}^{2} +c^{2} } =\sqrt {(x-x_{j} )^{2} +(y-y_{j} )^{2} +c^{2} }$

where $$c$$ is a predetermined shape parameter. Likewise, the derivatives can be expressed as

$\tag{6} \begin{array}{l} \frac{\partial \phi }{\partial x} =\frac{x-x_{j} }{\sqrt {r_{j}^{2} +c^{2} } } ,\quad \frac{\partial \phi }{\partial y} =\frac{y-y_{j} }{\sqrt {r_{j}^{2} +c^{2} } } \\ \frac{\partial^{2} \phi }{\partial x^{2} } =\frac{(y-y_{j} )^{2} +c^{2} }{(r_{j}^{2} +c^{2} )^{3/2}} ,\quad \frac{\partial^{2} \phi }{\partial y^{2} } =\frac{(x-x_{j} )^{2} +c^{2} }{(r_{j}^{2} +c^{2} )^{3/2}} \end{array}$

Substituting into the original equation set, one obtains

$\tag{7} \begin{array}{l} {\sum_{j=1}^{N}\nabla^{2} \phi (r_{j} ) T_{j} =f(\textbf x),\quad i=1,2,\cdots ,N_{I} } \\ {\sum_{j=1}^{N}\phi (r_{j} )T_{j} =g(\textbf x),\quad i=N_{I} +1,N_{I} +2,\cdots ,N} \end{array}$

For a parabolic equation, an implicit scheme can be used, i.e.,

$\tag{8} \frac{T^{n+1} -T^{n} }{\Delta t} -\alpha \left(\frac{\partial^{2} T^{n+1} }{\partial x^{2} } +\frac{\partial^{2} T^{n+1} }{\partial y^{2} } \right)=f(x,y,T^{n} ,\frac{\partial T^{n} }{\partial x} ,\frac{\partial T^{n} }{\partial y} )$

where $$\Delta$$t denotes the time step and superscript $$n+1$$ is the unknown (or new) value to be solved. The approximate solution can be expressed as

$T(x,y,t^{n+1} )=\sum_{j=1}^{N}\phi_{j} (x,y)T_{j}^{n+1}$

Substituting into Eq. (8), one obtains

$\begin{array}{l} {\sum_{j=1}^{N}\left(\frac{\phi_{j} }{\Delta t} -\alpha \left(\frac{\partial^{2} \phi }{\partial x^{2} } +\frac{\partial^{2} \phi }{\partial y^{2} } \right)\right)(x_{i} ,y_{i} )T_{j}^{n+1} =} \\ {\frac{\phi_{j}}{\Delta t} T^{n} (x_{i} ,y_{i} )+f(x_{i} ,y_{i} ,t^{n} ,T^{n} (x_{i} ,y_{i} ),T_{x}^{n} (x_{i} ,y_{i} ),T_{y}^{n} (x_{i} ,y_{i} ))\quad i=1,2,...,N_{I} } \\ {\sum_{j=1}^{N}\phi (x_{i} ,y_{i} )T_{j}^{n+1} =g(x_{i} ,y_{i} ,t^{n+1} )\quad i=N_{I} +1,...,N} \end{array}$

which produces an $$NxN$$ linear system of equations for the unknown $$T^{n+1}\ .$$ Note that

$T^{n} (x_{i} ,y_{i} )=\sum_{j=1}^{N}\phi_{j} (x_{i} ,y_{i} )T_{j}^{n} ,\quad T_{x}^{n} (x_{i} ,y_{i} )=\sum_{j=1}^{N}\frac{\partial \phi_{j}(x_{i} ,y_{i} )}{\partial x}{T_{j}^{n}} ,\quad T_{y}^{n} (x_{i} ,y_{i} )=\sum_{j=1}^{N}\frac{\partial \phi_{j} (x_{i} ,y_{i} )}{\partial y}{T_{j}^{n}}$ Figure 1: Irregular domain discretized using (a) 3-noded triangular finite elements, b) boundary element, and (c) arbitrary interior and boundary points using a meshless method.

Figure 1 shows an arbitrary domain discretized using 3-noded triangular elements, boundary elements, and a meshless method. An internal mesh is required in the FEM (Fig. 1a) and linear elements are needed along the boundary in the BEM (Fig. 1b). Both methods require the use of efficient matrix solvers to obtain values at the prescribed nodes, which can become resource limiting and time consuming. The meshless and boundary element methods, with arbitrarily distributed interior and boundary points, require no mesh as illustrated in Fig. 1 (b,c).

## Examples

To illustrate the use of meshless methods, let us begin with a simple heat transfer problem. The governing equation for energy can be written as

$\tag{9} \frac{\partial T}{\partial t} + \textbf V \cdot \nabla T = \alpha \nabla^2 T + Q$

$\tag{10} q + k \nabla T - h (T - T_\infty) - \varepsilon \sigma (T^4 - T^{4}_{\infty}) = 0$

$\tag{11} T(\textbf x,0) = T_0$

where $$\textbf V$$ is the vector velocity, $$\textbf x$$ is vector space, $$T(\textbf x,t)$$ is temperature, $$T_\infty$$ is ambient temperature, To is initial temperature, $$\alpha$$ is thermal diffusivity ($$\kappa /\rho cp$$), $$\epsilon$$ is emissivity, $$\sigma$$ is the Stefan-Boltzmann constant, $$h$$ is the convective film coefficient, $$q$$ is heat flux, and $$Q$$ is heat source/sink. Velocities are assumed to be known and typically obtained from solution of the equations of motion (a separate program is generally used for fluid flow).

For the 2-D transport equation for heat transfer, the relation becomes

$\tag{12} \begin{array}{l} {\sum_{j=1}^{N}\left(\frac{\phi (r_{j} )}{\Delta t} -\alpha \nabla^{2} \phi (r_{j} )\right) T_{j}^{n+1} =\frac{T^{n} (r_{j} )}{\Delta t} -V\bullet \nabla T^{n} (r_{i} ),\; \; i=1,2,\cdots ,N_{I} } \\ {\sum_{j=1}^{N}\phi (r_{j} )T_{j}^{n+1} =g(x,t),\; \; i=N_{I+1},N_{I+2},\cdots ,N} \end{array}$

from which one can solve the N X N linear system above for the unknowns $${T_{j}}$$ and obtain the approximate solution at any point in the domain, $$\Omega \ .$$

In this example problem, a two-dimensional plate is subjected to prescribed temperatures applied along each boundary , as shown in Fig. 2. The temperature at the mid-point $$(1,0.5)$$ is used to compare the numerical solutions with the analytical solution. The analytical solution is given as

$\theta (x,y)\equiv \frac{T-T_{1} }{T_{2} -T_{1} } =\frac{2}{\pi } \sum_{n=1}^{\infty }\frac{(-1)^{n+1} +1}{n} \sin \left(\frac{n\pi x}{L} \right)\frac{\sinh (n\pi y/L)}{\sinh (n\pi W/L)} :$

which yields $$\theta (1,0.5) = 0.445\ ,$$ or $$T(1,0.5) = 94.5\; \,^{\circ}{\rm C}$$ . Table 1 lists the final temperatures at the mid-point using a finite element method, a boundary element method, and a meshless method.

Table 1. Comparison of results for problem 1 from Exact, FEM, BEM, and Meshless Methods

 Method mid-point ($$\,^{\circ}{\rm C}$$) Elements Nodes Exact 94.512 0 0 FEM 94.605 256 289 BEM 94.471 64 65 Meshless 94.514 0 325

As a second example problem, a two-dimensional domain is prescribed with

Dirichlet and Neumann boundary conditions applied along the boundaries, as shown in Fig. 3. This problem, described in Huang and Usmani , was used to assess an h-adaptive FEM technique for accuracy. A fixed temperature of $$100\; \,^{\circ}{\rm C}$$ is set along side AB; a surface convection of $$0\; \,^{\circ}{\rm C}$$ acts along edge BC and DC with $$h = 750\; W/m^2\,^{\circ}{\rm C}$$ and $$k = 52\; W/m\,^{\circ}{\rm C}\ .$$ The temperature at point $$E$$ is used for comparative purposes. The severe discontinuity in boundary conditions at point $$B$$ creates a steep temperature gradient between points$$B$$ and $$E\ .$$ Figures 3(b,c) show the initial and final FEM meshes after two adaptations using bilinear triangles . The analytical solution for the temperature at point $$B$$ is $$T = 18.2535\; \,^{\circ}{\rm C}\ .$$

Table 2 lists the results for the three methods compared with the exact solution. The initial 3-noded triangular mesh began with 25 elements and 19 nodes.

Table 2. Comparison of results for problem 2 from Exact, FEM, BEM, and Meshless Methods

 Method Point E ($$\,^{\circ}{\rm C}$$) Elements Nodes Exact 18.2535 0 0 FEM 18.1141 256 155 BEM 18.2335 32 32 Meshless 18.2531 0 83

A simple irregular domain is used for the last example problem and results compared from the three methods. Results from a fine mesh FEM technique (without adaptation) are used as a reference benchmark [16,17]. The discretized domain and accompanying boundary conditions set along each surface are shown in Fig. 4. The FEM results are displayed as contour intervals. Figure 5(a,b) shows meshless results (using FEM fine mesh for contouring) versus FEM solutions using adapted quadrilateral elements. Heat conduction occurs as a result of constant temperatures set on the top and bottom surfaces, adiabatic faces in the upper right cutout and lower cutout portions, and convective heating along the right and left vertical walls. Adaptive meshing occurs in the corners as a result of steep temperature gradients; this is not evident when using meshless methods.

The FEM, BEM, and meshless mid-point values at $$(0.5,0.5)$$ are listed in Table 3.

Table 3. Comparison of results for problem 3 from FEM, BEM, and Meshless Methods

 Method mid-point ($$\,^{\circ}{\rm C}$$) Elements Nodes FEM 75.899 138 178 BEM 75.885 36 37 Meshless 75.893 0 96

Results obtained using a meshless method based on the Kansa approach compare very well to solutions obtained using an h-adapting finite element and boundary elements for heat transfer in two-dimensions. All three techniques provide accurate results for the three example cases. The meshless method was clearly the fastest, simplest, and least storage demanding method to employ. Advances now being made in meshless methods will eventually enable the scheme to compete directly with the FEM and BEM on a much broader range of problems. Details regarding the development and use of meshless methods can be obtained by accessing the web site http://www.ncacm.unlv.edu. Recent work using a meshless method for incompressible fluid flow with turbulence and heat transfer, including compressible flows with shock capture as well as bioengineering applications, are discussed in Pepper et al .