In this chapter we will describe about numerical simulation with meshfree methods. We know; phenomena in nature, whether physical, geological, mechanical, electrical, or biological, can often be describe by means of algebraic, differential, or integral equations. Obtaining exact solutions for these equations is ideal. Unfortunately, we can only obtain exact ones for limited practical problems because most of these problems are complex. Therefore using numerical procedure to obtain approximate solutions is inevitable. One of the most important tools in the field of numerical methods that has been developed newly is meshfree or meshless methods.
A meshfree method is a method used to establish system algebraic equations for the whole domain of problem without using a predefined mesh for the domain discretization. This infant method uses a set of scattered nodes, called field nodes, to establish the problem domain and boundaries, which do notrequire any priori information on the relationship between the nodes for the interpolation or approximation of the unknown functions of field variables. Inthe FEM, a continuum with a complicated shape is divided into elements,finite elements. The individual elements are connected together by atopological map called a mesh.
Meshfree methods have been proposed and achieved remarkable progress over the past few years. According to the formulation procedure, meshfree methods fall into three categories: meshfree weak form methods (like: EFG, MLPG, LRPIM,…), meshfree strong form methods (like: SPH, Collocation method,…) and meshfree weak-strong form methods based on the combination of both weak form and strong form (like: MWS method).These three categories and their limitations, applications, advantages and other descriptions will be introduced.In seeking for an approximate solution to the problem governed by PDEs and boundary conditions, one first needs to approximate the unknown field function using shape (trial or base) functions before any formulation procedure can be applied to establish the discretized system equations. In this chapter definition of base and shape functions and various techniques for meshfree shape function constructions are discussed. These shape functions are locally supported, because only a set of field nodes in a small local domain are used in the construction and the shape function is not used or regarded as zero outside the local domain. Such a local domain is termed the support domain or influence domain. The concept and kinds of support domain and determination of the dimension of the support domain will be described.
After introducing the concept of support domain, the point interpolation method (PIM) in detail will be discussed. Point interpolation method is one of the series representation methods for the function approximation, and useful for creating meshfree shape functions. A scalar function defined in the problem domain that is represented by a set of scattered nodes will be shown. There are two types of PIM shape functions have been developed so far using different forms of basis functions Polynomial basis functions and radial basis functions (RBF) have often been used in meshfree methods. These two types of PIMs will be discussed in the following chapter.
For satisfying the boundary conditions, penalty method, direct method, lagrange multiplier method and direct interpolation method can be used to enforce essential boundary conditions. One of these methods due to using meshfree method can be elected and they will be explained and compared.
To simulate some problems, the partial differential equations and boundary conditions for two dimensional solid mechanics and fluid mechanics problem and heat transfer problem especially thermodynamics of plates and shells will be given in sub-sections. These problems are solved with meshfree methods.
2. Meshfree methods categories
According to the formulation procedure, meshfree methods fall into three categories: meshfree weak form methods (like: EFG, MLPG, LRPIM,…), meshfree strong form methods (like: Collocation method, SPH,…) and meshfree weak-strong form methods based on the combinations of both weak forms and strong forms (like: MWS method).
2.1. Strong form methods
Regarding to formulate the governing equations, the direct approximate solution from the differential equations is used. It means the strong form of governing equations for boundary conditions are directly discretized at the field nodes to obtain a set of discretized system equations. If Taylor series is used and the differentiations are replaced, the method is the strong form method. The strong form method does not need the numerical integration. Thus
The background mesh even locally is not needed for the strong form methods.
Meshfree strong form methods have some attractive advantages: a simplealgorithm, computational efficiency, and truly meshfree. However, Meshfreestrong form methods are often unstable, not robust, and inaccurate,especially for problems with derivative boundary conditions. Severalstrategies may be used to impose the derivative (Neumann) boundary conditions in thestrong form methods, such as the use of fictitious nodes, the use of theHermite-type Meshfree shape functions, the use of a regular grid on the derivative boundary.
2.2. Weak form methods
In Meshfree weakformmethods, the governing partial differential equations (PDEs) withderivative boundary conditions are first transformed to a set of so called weakformintegral equations using different techniques. Theweak forms are then used to derive a set of algebraic system equations througha numerical integration process using sets of background cells that may beconstructed globally or locally in the problem domain.Meshfreeweak form methods were relatively under developed before 1990,but there has been a substantial increase in research effort since then.
There are now many different versions of Meshfreeweak form methods. Meshfreeweak form methods based on the global weak forms are called Meshfree global weak form methods, and those based on local weak forms are called Meshfree local weak form methods. Meshfree global weak form methods are based on the global Galerkinweak form for equations of problems and the Meshfree shape functions. Two typical Meshfree global weak form methods: the element freeGalerkin (EFG) method (Belytschko et al., 1994a) and the radial point interpolation method (RPIM) (GR Liu and Gu, 2001c; Wang and GR Liu, 2000; 2002a).
Another typical Meshfree global weak form method isthe reproducing kernel particle method (RKPM) proposed by Liu and coworkersin 1995 (Liu et al., 1995). The main idea of RKPM is to improvethe SPH approximation to satisfy consistency requirements using acorrection function. RKPM has been used in nonlinear and largedeformation problems (Chen et al., 1996; Chen et al., 1998; Liu and Jun,1998), inelastic structures (Chen et al., 1997), structural acoustics (Uras et al.,1997), fluid dynamics (Liu and Jun et al., 1997), et cetera.Meshfree local weak form methods were developed by Atluri andcoworkers based on the local Petrov-Galerkinweak form, and the Meshfreeshape functions.Some other Meshfreeweak form methods have also been developed, suchas the hp-cloud method (Armando and Oden, 1995), the partition of unityfinite element method (PUFEM) (Melenk and Babuska, 1996; Babuska andMelenk, 1997), the finite spheres method (De and Bathe, 2000), the freemesh method (Yagawa and Yamada, 1996), et cetera.
2.3. Weak-strong form methods
These Meshfree methods are called Meshfree weak-strong (MWS) formmethods in thisbook because are based on the combination of weak and strong form methods. The MWS method was developed by GR Liu and Gu(2002d, 2003b). The key idea of the MWS method is that in establishing thediscretized system equations, both the strong form and the local weak form are used for the same problem, but for different groups of nodes that carriesdifferent types of equations/conditions. The local weak form is used for allthe nodes that are on or near boundaries with derivative (Neumann) boundary conditions.The strong form is used for all the other nodes. The MWS method uses least background cells forthe integration, and it is currently the almost ideal Meshfree method that canprovide stable and accurate solutions for mechanics problems.
There are also Meshfree methods based on the integral representationmethod for function approximations, such as the Smooth ParticleHydrodynamics (SPH) methods (Lucy, 1977; Gingold and Monaghan, 1977;GR Liu and Liu, 2003, etc.). In the standard SPH method, the functionapproximation is performed in a weak (integral) form, but strong formequations are directly discretized at the particles.
2.4. Comparisons between three meshfree categories
Each meshfree method has features with advantages and defects. With these properties, the appropriate method can be selected to solve the problem. The features of methods are presented in sub-sections.
2.4.1. Advantages and disadvantages
Convergence rate and highest accuracy are important properties in numerical methods. When the problems include Dirichlet boundary conditions, the strong form methods are the best but in case of Neumann boundary conditions, weak form methods are optimum and when the both of Dirichletand Neumann boundaries are used in problems, the weak-strong form methods are useful.
The strong form methods are with good convergence rate andthey are truly meshless.The procedure is straightforward, and the algorithms and coding are simple. They are computationally efficient, and the solution is accurate when there are only Dirichlet boundary conditions.
However, Meshfreestrong form methods have disadvantages: they are often unstable and less accurate, especially for problemsgoverned by PDEs - with derivative boundary conditions.Derivative boundary conditions (DBCs) involve a set of separatedifferential equations defined on the boundary; these are different fromthe governing equations defined in the problem domain. These DBCsrequire special treatments.Unlike integration, which is a smoothing operator, differentiation is aroughening operator; it magnifies errors in an approximation. Thismagnified error is partially responsible for the instability of the solution of PDEs. Hence, Meshfreestrong formmethods are often unstable. Special treatments are employed to implement the derivative boundaryconditions in Meshfreestrong form methods. However, such treatmentscannot always control the error. Atechnique suitable for one problem may not work for another, even oneof the same types. A set of parameters tuned for one problem may notwork for another.
Thecommon feature of Meshfreeweak form methods is that the PDE of a problem is first replaced by or converted into an integral equation(global or local) based on a principle (weighted residual methods, energyprinciple). Weak form system equations can then be derived byintegration by parts.A set of system equations of Meshfreeweak form methods can be obtainedfrom the discretization of the weak form using meshfree interpolationtechniques.There are four features of the local weak form.The integral operation can smear the errorover the integral domainand, therefore improve the accuracy in the solution. It acts like somekind of regularization to stabilize the solution.The requirement of the continuity for the trial function is reduced orweakened, due to the order reduction of the differential operationresulting from the integration by parts.The force (derivative) boundary conditions can be naturallyimplemented using the boundary integral term resulting from theintegration by parts.The system equations in the domain and the derivative boundaryconditions are conveniently combined into one single equation.
These features give Meshfreeweak form methods the following advantages.They exhibit good stability and excellent accuracy for manyproblems. The derivative (Neumann) boundary conditions can be naturally andconveniently incorporated into the same weak form equation. Noadditional equations or treatments are needed and no errors areintroduced in the enforcement of traction boundary conditions.A method developed properly using a weak form formulation isapplicable to many other problems. A set of parameters tuned forone method for a problem can be used for a wide range of problems.This robustness of the weak form methods have been demonstratedthrough many practical problems. It is this robustness that makesthe weak form methods applicable to many practical engineeringproblems.
However, Meshfree global weak form methods are meshfree only in terms of the interpolation of the field variables. Background cells have to be used to integrate a weak form over the global domain. The numerical integration makes them computationally expensive, and the background mesh for the integration means that the method is not truly meshless.
In the Meshfree local weak form methods, the local integral domain in theinterior of the problem domain is usually of a regular shape. It can be assimple as possible and can be automatically constructed in the process ofcomputation. The Meshfree local weak form methods have obtainedsatisfactory results in solid mechanics and fluid mechanics (Atluri and Shen,2002; GR Liu, 2002).
Although the Meshfree local weak form methods made a significant step indeveloping ideal meshfree methods, the numerical integration is stillburdensome, especially for nodes on or near boundaries with complex shape.The local integration can still be computationally expensive for somepractical problems. It is therefore desirable to minimize the need fornumerical integrations.
The Mesh Weak-Strong method is designed to combine the advantages of strong formand weak form methods and to avoid their shortcomings. This can beperformed only after a thorough examination of the features of both types ofmethods, presented in the above sentences. AnMeshfree weak-strong (MWS) form method was proposed recently byGR Liu and Gu (2002d); it aimed to remove the background mesh forintegration as much as possible, and yet to obtain stable and accuratesolutions even for PDEs with derivative boundary conditions. The MWSmethod has been successfully developed and used in solid mechanics (Guand GR Liu, 2005; GR Liu and Gu, 2003b) and fluid mechanics (GR Liuand Wu et al., 2004; GR Liu and Gu et al., 2003c).
The convergence of the MWS method is studied numerically by comparison with other methods. The weak form method treats the Neumann boundary condition naturally and easily. In addition, the accuracy achieved by meshfree methods based on the weak form equations are generally much better than those based on strong form equations. However, the efficiency is a big problem for the weak form methods because of the need for weak form integration.
The MWS method proposed by Liu and Gu was based on both collocation and local radial point interpolation formulation. In the present MWS method, the strong form of meshfree collocation method is applied to the internal nodes and the nodes on the essential boundaries, while the local radial point interpolation weak form is applied to the nodes on the natural boundaries. The advantages of this MWS method are:
The Neumann boundary condition can be imposed straightforwardly and accurately with arbitrary nodal distributions.
Stable and accurate solution can be obtained with high efficiency.
2.4.2. Applications of each category
Strong form methods are suitable for Dirichlet boundary conditions problems and weak form methods are used more with problems that have Neumann boundary conditions. Weak-strong form methods are appropriate for problems with both of Dirichlet and Neumann boundary conditions.
3. Shape functions
In seeking for an approximate solution to a problem governed by PDEs and boundary conditions, one first needs to approximate the equation variables using shape functions, befor any formulation procedure can be applied to establish the descretized system equations.
This section discusses various techniques for MFree shapefunction constructions. These shape functions are locally supported, because only a set of field nodes in a small local domain are used in the constructionand the shape function is not used or regarded as zero outside the localdomain. Such a local domain is termed the support domain or influencedomain or smoothing domain.
3.1. Point interpolation methods shape functions
The point interpolation method (PIM) is one of the series representation methods for the function approximation, and is useful for creating Meshfree shape functions. Consider a scalar function T(x) defined in the problem domain Ω that is represented by a set of scattered nodes. The PIM approximates T(x) at a point of interest x in the form of
where the Bi(x) are the basis function defined in the space Cartesian coordinatesXT=[x, y],m is the number of basis functions, and the aiare the coefficients.
For function approximation, a local support domain is first formed for the point of interest at xwhich includes a total of n field nodes. For the conventional point interpolation method (PIM), n=m is used that results in the conventional PIM shape functions that pass through the function values at methods. The RPIM interpolation augmented with each scattered node within the defined support domain.
For the weighted least square (WLS) approximation or the moving least squares (MLS) approximation, n is always larger than m. There are two types of PIM shape functions have been developed so far using different forms of basis functions. Polynomial basis functions (GR Liu and Gu, 1999; 2001a) and radial basis functions (RBF) (Wang and GR Liu, 2000; GR Liu, 2002) have often been used in Meshfree methods.
3.1.1. Conventional polynomial PIM
Using polynomials as the basis functions in the interpolation is one of the earliest interpolationschemes. It has been widely used in establishing numerical methods, such as the FEM.Consider a continuous function u(x) defined in a domain, which is represented by a set of field nodes. The u(x) at a point of interestxis approximated in the form of
PT(x)=(1 x y x2 xy y2) for m=6, p=2 (2-D)
wherepi(x) is a given monomial in the polynomial basis function in the space coordinates xT=[x,y], m is the number of monomials, and ai is the coefficient for pi(x) which is yet to be determined. The pi(x) in Equation is built using Pascal's triangles, and a complete basis is usually (but not always)
3.1.2. Radial point interpolation shape functions
In order to avoid the singularity problem in the polynomial PIM, the radial basis function (RBF) is used to develop the radial point interpolation method (RPIM) shape functions for Meshfreeweak form methods (GR Liu and Gu, 2001c; Wang and Liu, 2000; 2002a,c). The RPIM shape functions will be used for both Meshfreeweak form and strong form polynomials can be written as
Where Rj(x) is a radial basis function (RBF), n is the number of RBFs, pi(x) is monomial in the space coordinates xT=[x, y], and m is the number of polynomial basis functions. When m=0, pure RBFs are used. Otherwise, the RBF is augmented with m polynomial basis functions. Coefficients ai and bj are constants. r is the distance between the point of interest (x,y) and a node (xi,yi) atr=
There are a number of types of radial basis functions (RBF), and the characteristics of RBFs have been widely investigated (Kansa,1990; Sharan et al.,1997;Franke and Schaback, 1997; etc). Four often used RBFs, the multi-quadrics (MQ) function, the Gaussian (Exp) function, the thin plate spline (TPS)function, and the Logarithmic radial basis function, are listed in Table.1.
Note: dcis a characteristic length that relates to the nodal spacing in the local support domain of the point of interest x, and it is usually the average nodal spacing for all the nodes in the local support domain.
the polynomial moment matrix is
The above equations are brought to show the procedure of shape function produce. The shape functions Φ are obtained and then the discretized derivatives can be used to governing equations and the parameters are shown with the equation
The derivatives of u(x) are easily obtained as
wherel denotes either the coordinates x or y.
The accuracy of interpolation for the point of interest depends on the nodes in the support domain as shown in Fig.2. Therefore, a suitable support domain should be chosen to ensure an efficient and accurate approximation. For a point of interest at Xα, the dimension of the support domain is determined by
α is the dimensionless size of the support domain, and dcis the nodal spacing near the point at. If the nodes are uniformly distributed, dcis simply the distance between two neighboring nodes.When nodes are non uniform and where α is a constant of shape parameter, dc can be defined as an average nodal spacing in the support domain of Xa. The exponential function of the support domain αs controls the actual dimension of the support domain.
Rectangular support domains (and r: dimensions of the support domain in x and y directions). The support domain is centred by Xq.
The actual number of nodes, n, can be determined by counting all the nodes included in the support domain. Generally, an α=2.0~3.0 leads to good results for many problems that we have studied. Note that the support domain is usually centered by a point of interest at xq.
4. Satisfying boundary conditions
For the Dirichlet boundary condition, the essential boundary conditions for ucan be simply given as follows: (when node is on the boundary)u=
The essential boundary condition can be directly imposed using the direct interpolation method. another method is the Penalty method has been used to enforce essential boundary conditions in the MLPG and LRPIM Methods. Since RPIM shape functions possess the Kronecker delta function property, the essential boundary conditions can be easily enforced as in the FEM (see, e.g., GR Liu and Quek, 2003).
The natural boundary conditions can besatisfy automatically when we use weak-strong form method and no additional equation or treatment is needed.
4.1. Direct method
The ith component is prescribed by setting
Such an essential boundary condition can then be enforced directly into the system Equation through the following modifications to the global matrix and the global right vector. The global matrix, K, is changed to
The components in the global right vector are changed to
The direct method can exactly enforce essential boundary conditions, but changing matrices and vectors needs additional computational operations. In addition, the algorithm of the direct method is also complicated.
The penalty method is a convenient alternative for enforcing the essential boundary conditions, in which the diagonal entry in the stiffness matrix, is changed to
whereis the penalty coefficient that is the much larger number than the components of the global matrix K. In the global right vector F, only the component Fiis changed as follows
The penalty method has some advantages: there are only two changes of matrices, and the algorithm is very simple. However, the penalty method can only approximately satisfy the essential boundary conditions. In addition, the accuracy is affected by selection of the penalty coefficient.
the global matrix, K, is then changed to
5. Examplesfor numerical simulations
In this section, some problems are brought to show the abilitiesof meshfree methods for solving the heat, solid and fluid mechanics problems.
5.1. Heat conduction
Meshfre methods are used to solve the heat transfer problem. For example, heat conduction in the plate is solved.
5.1.1.Formulation of heat transfer in the plate
Formulation of heat transfer in Cartesian coordinate is
Ifk(conductivity coefficient) is constant and for steady state and without any energy generation we have:
T is the Temperature and is the rate of energy generation and is the density and Cp is the specific heat in the formula.
5.1.2. Numerical results and discussion Domain representation for heat transfer
First, the temperature distribution in square plate is obtained. In problem 1, the bottom wall is in temperature T0 and other walls are in temperature 0 and in problem 2, the up wall has Neumann boundary condition. To check the validity of the method, three different problems are considered. Fig.3 shows the domain representation for problems 1 and 2 by the scattered nodes. The essential and natural boundary conditions should be satisfied on the boundary nodes.
5.1.3.Problem1 with essential boundary conditions
The constant temperature lines in this figure are shown by solid lines. Table.2 compares the results obtained by collocation method and LRPIM method with those obtained by the analytical method.
the analytical solution of the problem can be written as
(x,y) are the coordinate of points in the plate. T is the temperature. We showed the difference between three meshfree methods and the difference between using different number of nodes to give the better results.
We used the error norm
and Total Error defined as
5.1.4. Problem 2 with naturalboundary condition
The MWS method is used to solve the same problem with both essential and natural boundary conditions by 256(16x16) and 1156(34x34) nodes. The temperature distributions for those problems are given in Fig.7. It should be noted that the essential boundary conditions are satisfied exactly whereas the natural (Neumann) boundary conditions are satisfied in the weak form formulation.
5.2.Lid driven cavity problem
In this section, the lid driven cavity problem is solved by meshfree method. In this problem, the cavity is full of fluid and the upper plate in the cavity drive horizontally. It is shown that the moving boundary conditions in the top wall are easily applied and natural boundary condition can be satisfied.
5.2.1.Formulation and boundary conditions of driven cavity problem
The application of Navier-Stokes equation in solving fluid flow has evolved in the past few decades with meshfree method as one of the most adopted techniques. In this section the Navier-Stockes equation is solved:
The boundary conditions are shown in fig.8 and it is shown that three walls are without motion and the upper wall move with the fix speed.
5.2.2. Numerical results
The results from the solution of driven cavity problem when Reynolds number is 100 are shown below:
The vortex on the corner is created and is related to Reynolds number and the corner vortex will grow if Reynolds number increases. The vorticity and velocity contours are shown:
In Table 6 the results are compared with the minimum stream function and the location of large vortex. The results show the ability of meshfree methods to simulate the fluid mechanics problems.
5.3. Cantilever beam problem
Numerical studies are conducted for a cantilever beam that is often used for benchmarking numerical methods because the analytic solution for this problem is known. This problem is a sample of solid mechanics.
5.3.1. Formulation of cantilever beam problem
The equilibrium equation is used with the formula:
σ is the stress vector and is the body force vector components. The strain-displacement relations are another formula that are brought in two directions:
is the strain component and u is the displacement in the horizontal direction.
is the strain component and v is the displacement in the vertical direction.
The last equation is Hook’s law:
De is the matrix of elastic constant.
5.3.2. Numerical results
The analytical solution is obtained for displacement of points of beam:
u is the displacement of points in horizontal direction and Pis the force at the end of the beam. E is the elasticity modulus and υ is the poisson ratio and moment of inertia is I, D is the height and L is the length
The energy norm is defined to compare the results.
Ωis the problem domain and is the strain vector with numerical and analytical solutions.
These results show the capacity of MWS and LRPIM methods to solve the solid mechanics problems. The Table 7 shows the errors are minimums with comparison with the analytic solution.
This research has been mainly financed by Shiraz university. The author appreciate the support from Shiraz university.
The meshfree methods are numerical methods that can be used to solve the many different and complicated problems. The heat transfer problems, solid and fluid mechanics problems have been solved with meshfree methods.
Three categories are used to solve the problems. Strong form methods, weak form methods and weak-strong form methods(MWS) are meshfree categories. They can be used to solve the problems with Dirichlet and Neumann boundary conditions. For examples the heat conduction problem and lid driven cavity and cantilever beam are solved that they have different type of boundary conditions.Solutions are related to many parameters: the selected meshfree method, number of nodes, shape function parameters et cetera.
Nowadays, many changes are employed to different types of meshfree methods. The advantages are improved and the high convergence rate and high accuracy are accessible.
- Partial Differential Problems