Meshless methods for PDEs
Darrell Pepper (2010), Scholarpedia, 5(5):9838. | doi:10.4249/scholarpedia.9838 | revision #91485 [link to/cite this article] |
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 [3] and Fasshauer [4] discuss meshfree methods, implementation, algorithms, and coding issues for stress-strain problems; Liu [3] includes Mfree2D, an adaptive stress analysis software package available for free from the web, and Fasshauer [4] include MATLAB modules. Atluri and Shen [5] also produced a textbook that describes the meshless method in more detail, including much in-depth mathematical basis.
Contents |
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 [6].
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 [6] by directly collocating RBFs, especially multiquadric approximations (MQ). The use of MQ was first developed by Hardy [7] in 1971 as an interpolation method for modeling the earth’s gravitational field. Franke [8] 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, [9]).
More recently, the use of the Method of Fundamental Solutions (MFS) with the dual reciprocity method (DRM) and RBFs has appeared [10]. 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.
Radial Basis Functions
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 [14]; Fasshauer [4] 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 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 [15], 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}\) [6]. 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 [2], 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 [2]. 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 [18].
References
- Lewis, R. W., Morgan, K., Thomas, H. R., and Seetharamu, K. N. (1996): The Finite Element Method in Heat Transfer Analysis, J. Wiley & Sons: Chichester, UK.
- Huang, H-C. and Usmani, A. S. (1994): Finite Element Analysis for Heat Transfer, Springer-Verlag: London, UK.
- Liu, G. R. (2002): Mesh Free Methods: Moving Beyond the Finite Element Method, CRC Press, Boca Raton, FL.
- Fasshauer, G. E. (2007), Meshfree Approximation Methods with MATLAB, World Scientific Pub. Co., Hackensack, NJ.
- Atluri, S. N. and Shen, S. (2002): The Meshless Local Petrov Galerkin (MLPG) Method, Tech Science Press, Encino, CA.
- Kansa, E. J. (1990): Multiquadrics – A scattered data approximation scheme with applications to computational fluid dynamics II, Computers Math. Appl., 19, 8/9, pp. 147-161.
- Hardy, R. L. (1971), Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res., 76, pp. 1905-1915.
- Franke, R. (1982), Scattered Data Interpolation: Tests of Some Methods, Math. Comput., 48, pp. 181-200.
- Li, J., Hon, Y. C., and Chen, C. S. (2002): Numerical comparisons of two meshless methods using radial basis functions, Engineering Analysis with Boundary Elements, Elsevier Science Ltd., 26, pg. 205-225.
- Alves, C. J. S. and Chen, C. S. (2005), A new method of fundamental solutions applied to nonhomogeneous elliptic problems, Adv. Comput. Math., 23, 1-2, pp. 125-142.
- Chen, C. S., Kuhn, G. Li., J., and Mishuris, G. (2003), Radial basis functions for solving near singular Poisson problems, Comm. Num. Meth. Eng., 19, 5, pp. 333-347.
- Brebbia, C. A. and Dominguez, J. J. (1989), Boundary Elements: An Introductory Course, Computational Mechanics Pub., Southampton, UK, and McGraw-Hill, NY.
- Kassab, A. J., Wrobel, L. C., Bialecki, R. A., and Divo, E. A. (2006), Chapter 4: Boundary-Element Method, in Handbook of Numerical Heat Transfer, 2nd Ed., W. J. Minkowycz, E. M. Sparrow, and J. Y. Murthy (eds.), Wiley & Sons, NJ, pp. 125-165.
- Powell, M. J. D. (1992): The theory of radial basis function approximation in 1990, in Advances in Numerical Analysis, Vol. II, W. Light (Ed.), Oxford Sci. Pub., Oxford, UK, pp. 105-210.
- Incropera, F. P. and DeWitt, D. P. (2002): Fundamentals of Heat and Mass Transfer, 5th Ed., J. Wiley & Sons: New York.
- Pepper, D. W., Carrington, D. B., and Gewali, L. (2000): A Web-based, Adaptive Finite Element Scheme for Heat Transfer and Fluid Flow. ISHMT/ASME 4th Conf. on Heat and Mass Transfer, Jan. 12-14, Pune, India.
- Pepper, D. W. and Chen, C. S. (2002): A Meshless Method for Modeling Heat Transfer, S. N. Atluri and D. W. Pepper (Eds.), Proceedings of the 22nd ICES 2002 Conference, July 31-August 2, 2002, Reno, NV.
- Pepper, D. W., Kassab, A. J., and Divo, E. A. (to appear 2010): Finite Element, Boundary Element, and Meshless Methods, ASME Press, NY.
Internal references
- Lawrence M. Ward (2008) Attention. Scholarpedia, 3(10):1538.
- Rob Schreiber (2007) MATLAB. Scholarpedia, 2(7):2929.
- Andrei D. Polyanin, William E. Schiesser, Alexei I. Zhurov (2008) Partial differential equation. Scholarpedia, 3(10):4605.
- Roberto Benzi and Uriel Frisch (2010) Turbulence. Scholarpedia, 5(3):3439.