A Review on the Modified Finite Point Method

Abstract

The objective of this paper is to make a review on recent advancements of the modified finite point method, named MFPM hereafter. This MFPM method is developed for solving general partial differential equations. Benchmark examples of employing this method to solve Laplace, Poisson, convection-diffusion, Helmholtz, mild-slope, and extended mild-slope equations are verified and then illustrated in fluid flow problems. Application of MFPM to numerical generation of orthogonal grids, which is governed by Laplace equation, is also demonstrated.

1. Introduction

In recent decades, a new type of numerical methods for solving partial differential equations (PDEs) without meshes has been developed. These methods are called meshless or meshfree methods, which refer to the idea that information of connectivity of each mesh or element is no longer required as the input data in this new development. It saves a lot of manpower in the tedious preprocessing stage in preparing numerical computations. However, the information of the discretized nodes, such as the position vectors, is still a necessary requirement.

The idea of meshless methods can be traced back to the smoothed particle hydrodynamics (SPH) for modeling the astrophysical phenomenon [1]. The basic concept of meshless methods is to approximate the solution of a PDE as a linear combination of basis functions or shape functions. The process for solving the PDE is to seek the undetermined factors which are mostly just coefficients of the basis functions and in some other cases could be the values of the solution itself at discretized nodes.

Meshless methods are categorized mainly into two parts [2], namely, methods based on weak-form formulations and methods based on collocation techniques, also known as the strong-form methods. Though strong-form methods are much simpler and more straightforward, it was reported in [2] that there has been less research devoted to strong-form methods. This is partly because they are less robust and partly because researches were concentrated on the finite element method (FEM) which used weak-forms. It was then a natural step to favor the use of weak-form methods.

In the 1990s, a general collocation method to solve linear partial differential equations was proposed by adopting the radial basis functions (RBFs) [3]. There are two kinds of RBF collocation methods, domain-type and boundary-type. The application of domain-type RBF collocation methods is broader, including tissue engineering [4, 5], heat transfer [6], nonlinear Burgers’ equations [7], shallow water equations [8–10], Stokes flows [11], electromagnetic problems [12, 13], and finance [14, 15]. However, it was reported that domain-type RBF collocation methods usually malfunction in the seeking of partial derivatives around the boundaries [16, 17]. On the contrary, boundary-type RBF collocation methods, in which governing equations are satisfied automatically in the entire domain, are more applicable to problems in which accurate partial derivatives of the solutions around the boundaries are desirable. Typical examples are the application of MFS (method of fundamental solutions) to water wave problems [18–20]. Nevertheless, boundary-type RBF collocation methods are limited to problems governed by some specific equations, such as Laplace or Helmholtz equations. In [21], a modified domain-type RBF collocation method was proposed. With additional satisfaction to the governing equation on boundary collocation points, the gradient of the velocity potential at any free surface node, which represents the velocity vector at that specific node, was accurately estimated. Thus this method was successfully applied to the simulation of water waves in a swaying tank.

Apart from the combination of RBFs, one could use a polynomial to approximate the solution of a PDE and seek the coefficients of all the monomials by applying the collocation technique. However, polynomials could not be used in a global solution form such as what is usually did in RBF collocation methods. This is because the values of high degree terms could become extremely large at points far from the origin. In most polynomial collocation methods, polynomials are just applied to approximate the solution piece-wisely around discretized nodes. The moving least square (MLS) or weighted least square (WLS) approaches are always accompanied by the polynomial approximation to transform and combine the locally approximated solutions into a global solution form. In each local solution, basis functions are the monomials and their factors are just the coefficients. The least square (LS) approach with Taylor series expansion was developed for the finite difference with arbitrary meshes [22] while the WLS and MLS are introduced to the polynomial collocation methods in [23, 24] which stated that in most cases local approximation with just LS is not stable. The WLS and MLS are identical to each other when they are used for the local approximation at discretized nodes. The difference is on how to calculate the partial derivatives of the approximation. In [23] one could find that both the performances of MLS and WLS are acceptable.

In the global solution form, the basis functions are called the shape functions and their factors are the values of the solution itself. This localizing approach is somewhat like localized RBF collocation methods [25]. Representatives of this family are the hp-meshless cloud method [26] and the finite point method (FPM) [23, 24].

Like most strong-form methods, collocation methods with local polynomial approximation are reported to be unstable [27] and usually can work only for problems with simple and regular geometry and boundary conditions. This coincides with what was reported in [2]. Contrarily, weak-form methods with MLS or WLS approaches [28–35], which need background cells for local integration, yield much higher stability and are applied more widely in many numerical implementations, such as viscous flow problems [31], nonlinear water wave problems [32, 33], and wave breaking problems [34, 35]. In some studies efforts are made to generalize the FPM to weak-form formulation by applying alternative weights [36] or alternative derivative approximation [37, 38] in the FPM. A new meshless method, based on subdomain collocation approach, is also present and named “the generalized finite point method” (GFPM) [38].

Though sometimes unstable, the simplicity renders the strong-form methods worthy of further investigations. Studies on improving the collocation techniques and the construction of local clouds, which means the selection of neighboring nodes involved in the local approximation, are important issues for increasing the numerical stability [39–41].

Besides methods mentioned above, one could also identify a few literatures about meshless methods employing the reproducing kernel approximation [42–46]. With complicated formulation, these methods accurately approximate PDE’s solutions as linear combinations of shape functions in a global way. In [44] it was reported that “satisfying the boundary conditions and the governing equations simultaneously at the boundary can provide a more accurate solution to a given problem.” This coincides with what was reported in [21]. Improvements of the approximation for making the governing equation satisfied at the boundaries were also proposed in [44]. Though in polynomial collocation methods with MLS approach the global solution forms are also the linear combinations of shape functions, the basic concept differs from the reproducing kernel approximation so that the way for collocation is not the same. In methods with reproducing kernel approximation, such as finite cloud methods (FCM), the partial derivatives of the solution are obtained from the differentiation of the global shape functions while in collocation methods with local polynomial approximation they are acquired by differentiating the local basis functions, which could turn out to be just the coefficients of higher monomials as part of the numerical solutions when the WLS approach is employed.

Adopting the advantage of localization, a robust local polynomial collocation method, which is similar to some finite point methods in which polynomials are localized by placing their origins at collocation points [39–41], was proposed in [47]. Inspired by [21, 44], the satisfaction of the governing equation is additionally required on boundary collocation points. This makes the method proposed in [47] not just as simple and straightforward as other similar methods, but also more robust and stable. Because this method mainly originates from the finite point method but improved, we hereby nominate it as the modified finite point method (MFPM). This method has been employed to several engineering and science applications [21, 48–53]. This paper is to review the modified finite point method, particularly emphasizing on its improvements and recent findings. Verifications of MFPM are illustrated in benchmark problems, and its recent applications to engineering and science problems are also demonstrated.

2. Formulations of Modified Finite Point Method

The local polynomial collocation method proposed in [47], which originates from the finite point method, is a mesh-free method for solving general partial differential equations (PDEs). A mesh-free method does not require the connectivity of nodal points of a mesh or element. Instead, the position vectors of computational nodes are the only information required to be provided as the input in the preprocessing stage of computations. This advantage of mesh-free methods renders them increasingly popularity in selection of numerical methods.

Taking a general 2-D linear second order PDE as an example, the descriptions of this method and the conventional method are given as follows.

2.1. The Local Approximation

When solving a general 2-D linear, second order PDE is ℒ𝜙=𝑐₁𝜙+𝑐₂𝜕𝜙𝜕𝑥+𝑐₃𝜕𝜙𝜕𝑦+𝑐₄𝜕²𝜙𝜕𝑥²+𝑐₅𝜕²𝜙𝜕𝑦²+𝑐₆𝜕²𝜙𝜕𝑥𝜕𝑦=𝑠 (1) subjected to specified boundary conditions ℬ𝜙=𝑞₁𝜙+𝑞₂𝜕𝜙𝜕𝑥+𝑞₃𝜕𝜙𝜕𝑦=𝑓⇀−𝑥∈Γ₁,𝜙=𝜙_𝑏,⇀−𝑥∈Γ₂, (2)(3) where 𝑐₁,𝑐₂,…,𝑐₆, 𝑞₁, 𝑞₂, 𝑞₃, 𝑓, and 𝑠 are all functions of 𝑥, 𝑦 and may change in time in transient problems. The boundary Γ₁ could be nonsmooth and then at a corner there could be more than one condition like (2). Therefore, 𝑞₁, 𝑞₂, 𝑞₃, and 𝑓could be multivalued in this case. Boundary condition expressed as (2) is just for conciseness. It will be explained later on how boundary conditions will be treated at a point where more than one boundary condition exists. In seeking the numerical solutions, the entire domain is distributed with 𝑁 nodes as needed. At each node, 𝑥=𝑥_𝑗, 𝜙 is approximated as 𝜙⇀−𝑥_{𝑥−⇀𝑥≈−⇀𝑗}≈𝜙_𝑗⇀−𝑥=^𝑚_𝑖=1𝛼_𝑗𝑖𝑝_𝑖𝑋⇀− (4) in which 𝑥⇀−𝑋=⇀−𝑥−⇀−_𝑗 is the relative position vector, 𝑝_𝑖(⇀−𝑋) is the 𝑖th monomial of the polynomial, and 𝛼_𝑗𝑖 are coefficients to be determined. The subscript 𝑗 indicates that this approximation is valid only in the vicinity of ⇀−𝑥_𝑗. Once a new ⇀−𝑥_𝑗 is chosen, there will be a new set of 𝛼_𝑗𝑖. For a 2-D problem governed by a second order partial differential equation, the monomials are 𝑝_𝑖𝑋=,𝑖=1∼𝑚1𝑋𝑌𝑋²𝑌²𝑋𝑌𝑋³𝑋²𝑌𝑋𝑌²𝑌³⋯ (5) in which 𝑗⇀−𝑋=𝑋⇀−𝑖+𝑌⇀−. The value of 𝑚 which represents that the chosen degree of the polynomial should be at least the same as the degree of the polynomial. It has been suggested [47] that the degree of the polynomial higher than the order of the differential equation by one might be better. The approximations of the solution and its partial derivatives are just formulated simply by the local approximation. In 2-D problems for example, they are 𝜙_𝑥^−⇀_𝑥=^−⇀_𝑗≈𝜙_{𝑥−⇀𝑥=−⇀𝑗}=𝛼_𝑗1,𝜕𝜙𝜕𝑥_𝑥^−⇀_𝑥=^−⇀_𝑗≈𝜕𝜙𝜕𝑥_{𝑥−⇀𝑥=−⇀𝑗}=𝛼_𝑗2,𝜕𝜙𝜕𝑦_𝑥^−⇀_𝑥=^−⇀_𝑗≈𝜕𝜙𝜕𝑦_{𝑥−⇀𝑥=−⇀𝑗}=𝛼_𝑗3,𝜕²𝜙𝜕𝑥²_{𝑥−⇀𝑥=−⇀𝑗}≈𝜕²𝜙𝜕𝑥²_{𝑥−⇀𝑥=−⇀𝑗}=2𝛼_𝑗4,𝜕²𝜙𝜕𝑦²_{𝑥−⇀𝑥=−⇀𝑗}≈𝜕²𝜙𝜕𝑦²_{𝑥−⇀𝑥=−⇀𝑗}=2𝛼_𝑗5,𝜕²𝜙𝜕𝑥𝜕𝑦_{𝑥−⇀𝑥=−⇀𝑗}≈𝜕²𝜙𝜕𝑥𝜕𝑦_{𝑥−⇀𝑥=−⇀𝑗}=𝛼_𝑗6,…. (6) Here the error residual of the local approximation around ⇀−𝑥=⇀−𝑥_𝑗 is defined as 𝐸_𝑗=^𝑁_𝑙=1𝑊_𝑗𝑙𝜙⇀−𝑥_𝑙−𝜙_𝑗⇀−𝑥_𝑙², (7) where 𝑊_𝑗𝑙 is a weighting factor determined by the distance between ⇀−𝑥_𝑗 and ⇀−𝑥_𝑙. Usually, a normalized Gaussian function is selected for determining the weighting factor 𝑊_𝑗𝑙=exp−𝜀𝑟_𝑗𝑙/𝜌_𝑗²−exp(−𝜀)1−exp(−𝜀),𝑟_𝑗𝑙<𝜌_𝑗0,𝑟_𝑗𝑙≥𝜌_𝑗, (8) where 𝑟_𝑗𝑙 is the distance between ⇀−𝑥_𝑗 and ⇀−𝑥_𝑙 (i.e., 𝑟_𝑗𝑙=|⇀−𝑥_𝑙−⇀−𝑥_𝑗|), 𝜀 is the shape parameter, and 𝜌_𝑗 is the supporting range measured from the point of  ⇀−𝑥_𝑗. The weight, 𝑊_𝑗𝑙, is intended to keep this parameter roughly within a range of optimal choice, regardless of the physical size of computational domain. Considering only the nonzero terms, (7) can be rewritten as 𝐸_𝑗=^𝑛_𝑘=1𝑊_𝑗𝑘𝜙⇀−𝑥_𝑘−𝜙_𝑗⇀−𝑥_𝑘², (9) where 𝑘 is the local index number of ⇀−𝑥_𝑙 in the 𝑗th subdomain and 𝑛 is number of nodes inside the subdomain for local collocation. The nodal index relationship between local and global domains is illustrated in Figure 1. Other alternative functions for considering the value of the weight can be found in references of the meshless methods using the MLS or WLS approach, both in weak-form and strong-form formulations [22–24, 26–41]. The weights are treated as functions in the MLS approach while they are just factors in seeking the approximations of the derivatives when WLS approach is employed.

Figure 1 

Illustration of the nodal index relationship between local and global domains.

2.2. The Conventional Collocation Method

Since the modified finite point method is originated from the conventional finite point method, the collocation method in the conventional finite point method will be presented briefly for completeness in this section.

The coefficients 𝛼_𝑗𝑖 corresponding to the minimized 𝐸_𝑗 in (9) can be expressed explicitly by the least square approach as 𝛼_𝑗=[Λ]𝛽 (10) in which 𝛼_𝑗𝑚×1=𝛼_𝑗1⋯𝛼_𝑗𝑖⋯𝛼_𝑗𝑚^𝑇[Λ]_𝑚×𝑛=[𝐴]^𝑇[𝐴]⁻¹[𝐴]^𝑇[𝐴]_𝑛×𝑚=𝑎₁₁𝑎₁₂⋯⋯𝑎_1𝑚𝑎₂₁⋱𝑎_2𝑚⋮𝑎_𝑘𝑖⋮𝑎⋮⋱⋮_𝑛1𝑎_𝑛2⋯⋯𝑎_𝑛𝑚𝛽_𝑛×1=𝛽₁⋯𝛽_𝑘⋯𝛽_𝑛^𝑇=𝑤₁𝜙₁⋯𝑤_𝑘𝜙_𝑘⋯𝑤_𝑛𝜙_𝑛^𝑇, (11) where 𝑤_𝑘=𝑊_𝑗𝑘,  𝜙_𝑘=𝜙(⇀−𝑥_𝑘), and 𝑎_𝑘𝑖=𝑤_𝑘𝑝_𝑖(⇀−𝑥_𝑘−⇀−𝑥_𝑗).

It should be noted that this approximation is only valid in the vicinity of ⇀−𝑥_𝑗. Once a new ⇀−𝑥_𝑗 is chosen, the entries of matrix [Λ] and components of vector [𝛽] should all be renewed.

When employing a conventional collocation method to solve a PDE, what one needs to do is just use the approximations of the solution and its partial derivatives to formulate a global matrix system and solve it. In the process of collocation, because it is anticipated to find local approximations of the solution and its partial derivatives very close to the relevant exact values which are not known a priori, the approximations listed in (6) are deemed as the exact values. This is similar to what is done in finite difference methods (FDM). Therefore, one gets a global matrix equation when assembling the local approximations into the global system [𝐾]_𝑁×𝑁𝜙_𝑁×1=[𝑏]_𝑁×1. (12) In the case that ⇀−𝑥_𝑗 is on the boundary and the value of 𝜙 is known, the entities in (12) are 𝜅_𝑗𝑙=1,if𝑗=𝑙0,otherwise𝑏_𝑗=𝜙_𝑗. (13)(14) When the collocation point is inside the domain, the entities in (12) are 𝜅_𝑗𝑙=𝑤_𝑘𝑐₁𝜆_1𝑘+𝑐₂𝜆_2𝑘+𝑐₃𝜆_3𝑘+2𝑐₄𝜆_4𝑘+2𝑐₅𝜆_5𝑘+𝑐₆𝜆_{6𝑘𝑥−⇀𝑥=−⇀𝑗},if⇀−𝑥_𝑗−⇀−𝑥_𝑙<𝜌_𝑗𝑏0,otherwise_𝑗=𝑠|_𝑥^−⇀_𝑥=^−⇀_𝑗, (15) where 𝜆_𝑖𝑘 indicates the 𝑖th entry in 𝑘th column in matrix [Λ]. When the collocation point is on the boundary where a Neumann condition or a Robin condition is given, the entities in (12) are 𝜅_𝑗𝑙=𝑤_𝑘𝑞₁𝜆_1𝑘+𝑞₂𝜆_2𝑘+𝑞₃𝜆_{3𝑘𝑥−⇀𝑥=−⇀𝑗},if⇀−𝑥_𝑗−⇀−𝑥_𝑙<𝜌_𝑗𝑏0,otherwise_𝑗=𝑓_𝑥^−⇀_𝑥=^−⇀_𝑗. (16)

2.3. The Modified Collocation Method

In the conventional method, when the collocation point is on the boundary, the satisfactory of the governing equation is not used. When the collocation point is on the boundary where two or more conditions are given, only one of them is chosen for the conventional collocation. The partial derivatives are thus inaccurately estimated as a result of lacking the complete satisfaction of all equations. Inspired by [21, 44], an alternative local approximation for direct collocation was proposed in [47]. The coefficients of the local polynomial corresponding to the minimal error residual at the node ⇀−𝑥_𝑗 (9) under the condition that ℒ𝜙−𝑠²+ℬ₁𝜙−𝑓₁²+⋯+ℬ_𝑛nd𝜙−𝑓_𝑛nd²⟶0, (17) where 𝑛_nd is the number of non-Dirichlet boundary conditions at the node ⇀−𝑥=⇀−𝑥_𝑗, can be expressed as 𝛼_𝑗=Λ^′𝛽𝛽^′ (18) in which Λ^′_{𝑚×(𝑛+𝑛nd+1)}=𝐴𝐴^′𝑇𝐴𝐴^′−1𝐴𝐴^′𝑇𝐴^′_{(𝑛nd+1)×𝑚}=𝑤^′𝑐₁⋯⋯⋯𝑤^′𝑐_𝑝𝑤0⋯0_′𝑞_1,1⋯𝑤^′𝑞_1,3𝑤0⋯⋯⋯0⋮⋮⋮⋮_′𝑞_𝑛nd,1⋯𝑤^′𝑞_𝑛nd,3𝛽0⋯⋯⋯0_{′(𝑛nd+1)×1}=𝑤^′𝑠𝑤^′𝑓₁⋯𝑤^′𝑓_𝑛nd^𝑇, (19) where the definitions of 𝑤_𝑘, 𝜙_𝑘, 𝑎_𝑘𝑖 are all the same as in the conventional collocation method, while 𝑤^′=𝑊^′ in which 𝑊^′ represents a penalty weighting factor whose value is much greater than 1. In the case of 𝑛_nd in (17) is greater than 1, it obviously indicates that the collocation point rests on an edge or at a corner. At an interior node, there is only one term in (13) (i.e., 𝑛_nd=0).

Assembling the local approximations into a global matrix system, one also gets a global matrix equation as (12), but the entries in this matrix equation are different from those described in (13)–(16), except those specified by the Dirichlet boundary condition. At all other points no matter inside the domain or on the boundaries, where non-Dirichlet conditions are given, the entries in (12) are 𝜅_𝑗𝑙=𝑤_𝑘𝜆^′_1𝑘−𝛿_𝑗(𝑙),if⇀−𝑥_𝑗−⇀−𝑥_𝑙<𝜌_𝑗𝛿0,otherwise_𝑗(𝑙)=𝑏1,if𝑗=𝑙0,otherwise_𝑗=−^𝑛nd+1_𝑘=1𝜆^′_{()1,𝑛+𝑘}𝛽^′_𝑘, (20)

where 𝜆^′_1𝑘 indicates the first entry in kth column in matrix [Λ^′].

2.4. The Parameters

Both in the conventional finite point method and in the modified one, there are several parameters to be chosen. They are the shape parameter 𝜀 in the normalized Gaussian function as the weight, the subdomain size, 𝜌_𝑗, and the degree of polynomial for the local approximation. The penalty weighting factor 𝑊^′ needs to be additionally given in MFPM. In [47] the shape parameter 𝜀 in the normalized Gaussian function for the modified method has been shown to be much more insensitive than the one needed in the conventional method. It was suggested that 𝜀=22, 𝜌_𝑗 as 1.05 times of the distance from ⇀−𝑥_𝑗 to its 25th nearest neighboring node, and 𝑊^′=10⁴ be chosen for numerical computations in different problems in this paper. The third degree polynomial gives better performance in both the conventional and the modified method for second order partial differential equations.

3. Verifications of Benchmark Problems

In this section, the following problems of different form of partial differential equations are chosen as benchmark problems to verify the accuracy of MFPM.

3.1. Laplace Equation

A potential uniform flow passing around a single cylinder in an infinite 2-D domain is chosen as the benchmark. The radius of the cylinder is 𝑎. Viscous effects which could be significant in real fluids are not considered here. The velocity at 𝑥=±∞ or 𝑦=±∞ is ⇀−𝑣=𝑈𝑗⇀−𝑖+0⇀−. The velocity potential 𝜙 is governed by the Laplace equation: 𝜕²𝜙𝜕𝑥²+𝜕²𝜙𝜕𝑦²=0. (21) The velocity field is determined as ⇀−𝑣=𝑢𝑗⇀−𝑖+𝑣⇀−: 𝑢=𝜕𝜙𝜕𝑥𝑣=𝜕𝜙.𝜕𝑦 (22) There is an exact solution expressed as 𝜙_𝑒𝑎=𝑈𝑥1+²𝑥²+𝑦². (23) The subscript “𝑒” denotes the “exact solution.”

Due to symmetry, only the upper left quadrant is modeled. The values of 𝑈 and 𝑎 are chosen as 1.0 and 0.5, respectively. These quantities are dimensionless. The computational domain is within the range of −2≤𝑥≤0 and 0≤𝑦≤2. The specification of the boundary conditions and the arrangement of the collocation points are shown in Figure 2. Two nodal arrangements shown in Figure 3 are chosen for the verifications. The velocity potential at the upper boundary in Figure 2 is specified by the analytical solutions along this boundary. No-flux boundary condition is applied to the lower boundary, due to symmetry of the problem on the straight line and physical solid boundary on the quad-circle curve, respectively.

Figure 2 

Boundary conditions for a problem of uniform flow passing a circular cylinder.

(a)

(b)

Figure 3 

Two types of nodal arrangement for a potential uniform flow passing a circular cylinder.

The profile of computed flow speed (𝑉=√𝑢²+𝑣²) along the cylinder surface is shown in Figure 4. It is obvious that the modified method performs much better than the conventional one. In the comparison between these two panels in Figure 4, one could find that nodal distribution is much less sensitive in the modified method than in the conventional one. In the right panel of Figure 4, the result by using the conventional collocation method is actually acceptable. This indicates that the conventional collocation method functions properly when the nodes are carefully arranged. This observation coincides with what was reported in [29]: polynomial collocation methods with MLS approach usually can work only for problems with simple and regular geometry and boundary conditions. It should also be noted that, at 𝜃=𝜋, which is a stagnation point, the result of the MFPM solution is exact while the result of the conventional collocation method is incorrect. This is because in the conventional collocation method only one Neumann boundary condition is employed in the collocation but in the MFPM the Neumann boundary conditions of both directions and the satisfactory of the governing equation are all used. This comparison has been done in [47]. Here we replot the figures and put them in the corresponding positions for the convenience of comparison.

(a)

(b)

Figure 4 

Computed flow speed profile along the body surface in the problem of a uniform flow passing a circular cylinder by two different nodal arrangements.

3.2. Poisson Equation

Consider the following 2-D Poison equation: 𝜕²𝜙𝜕𝑥²+𝜕²𝜙𝜕𝑦²=sin(𝜋𝑥)sin𝜋𝑦 (24) within the domain of 0≤𝑥≤1, 0≤𝑦≤1. The boundary conditions are 𝜕𝜙=𝜕𝑦−1()2𝜋sin𝜋𝑥cos𝜋𝑦,along𝑦=0,𝑦=1𝜙=0,along𝑥=0,𝑥=1. (25)

The exact solution of this Poisson equation can be found as 𝜙_𝑒=−sin(𝜋𝑥)sin(𝜋𝑦)/(2𝜋²).

The verification of the MFPM to this problem was demonstrated in [47]. Here we redo this problem by increasing the nodal density and plot the contours of the numerical result together with the exact solution for the convenience of comparison. A uniform nodal arrangement is chosen for numerical computations of this problem. The nodal spacing is 0.05. Therefore, totally 441 nodes are used in this test. The comparison of the numerical results with the analytical solutions is shown in Figure 5.

(a)

(b)

Figure 5 

Comparisons of the numerical results (a) with the exact solutions (b) for the benchmark problem of Poisson equation. From top to bottom, they are the value of the solution and partial derivatives in 𝑥- and 𝑦-direction, respectively.

3.3. Convection-Diffusion Equation

Consider the following steady-state convection-diffusion equation [54]: 𝜕²𝜙𝜕𝑥²+𝜕²𝜙𝜕𝑦²−𝑈𝜕𝜙𝜕𝑥=0intheregionof0≤𝑥≤1,0≤𝑦≤0.6, (26) with boundary conditions: 𝜙_𝑥=0=1,𝜙_𝑥=1=2,𝜕𝜙𝜕𝑦_𝑦=0=0,𝜕𝜙𝜕𝑦_𝑦=0.6=0. (27) This problem has an exact solution as 𝜙_𝑒[𝑈(=2−1−exp𝑥−1)]().1−exp−𝑈 (28)

A strong convective case of 𝑈=100 is chosen. The nodal arrangement for this problem is shown in Figure 6. The nodal spacing is mainly Δ𝑥=Δ𝑦=0.02 while near the right boundary the nodal density is increased to Δ𝑥=Δ𝑦=0.01, in order to catch the behavior of exponential growth near the boundary of 𝑥=1. Totally 1,825 nodes are used in this test. Only 1.48% of the entries in the global matrix [𝐾] are nonzero. Profiles of 𝜙 and 𝜕𝜙/𝜕𝑥 along 𝑦=0 and 𝑦=0.3 are shown in Figure 7. Very good agreement between present numerical results and the exact solutions of function values and its partial derivatives is obtained. The verification of the MFPM to this problem was demonstrated in [47]. Here we replot the figures of this problem for the convenience of comparison.

Figure 6 

Nodal arrangement for the problem of the steady-state convection-diffusion equation.

(a)

(b)

Figure 7 

Results of problem of the steady-state convection-diffusion equation.

3.4. Helmholtz Equation

In previous three benchmark problems, unknowns of real-valued governing equations are dealt with. There are physical problems governed by equation with complex-valued unknowns. It is intended to choose free-surface water waves problems to demonstrate applicability of MFPM to such problems.

In the constant water depth, the motion of progressive, monochromatic free-surface waves can be reduced to a 2-D problem governed by Helmholtz equation 𝜕²𝜙𝜕𝑥²+𝜕²𝜙𝜕𝑦²+𝑘²𝜙=0, (29) where 𝜙, a complex valued variable, is the pseudo velocity potential of the fluid on the free surface. The symbols used in physical problems hereafter should be differentiated from those used in the formulation of the MFPM. Wave number 𝑘 is defined as 2𝜋/𝐿, in which 𝐿 is the wavelength of the monochromatic waves. The real free-surface velocity potential of monochromatic waves is defined as Φ=Re𝜙𝑒^{−𝑖𝜔𝑡} (30) in which 𝑖=√−1 is the pure imaginary number. The symbol 𝜔 indicating the angular frequency is defined as 𝜔=2𝜋/𝑇 and 𝑇 is the period of the monochromatic wave motion.

In the problem of monochromatic waves diffracted by a vertical cylinder in a constant water depth, shown as in Figure 8, the outer boundary has to be located far away from the cylinder and is given as the Sommerfeld [55] radiation condition of the scattered waves, 𝜙^𝑆: lim_{√𝑟→∞}𝜕𝜙^𝑆𝜕𝑟−𝑖𝑘𝜙^𝑆=0. (31)

Figure 8 

A sketch of progressive monochromatic waves scattered by a vertical cylinder in a constant water depth.

For linear progressive, monochromatic waves, the total velocity potential of wave motion equals the linear combination of incident wave and scattered waves, that is, 𝜙=𝜙^𝐼+𝜙^𝑆, and (31) can be implemented at an artificial open boundary far away from the scatter as 𝜕𝜙𝜕𝑟−𝑖𝑘𝜙=𝜕𝜙^𝐼𝜕𝑟−𝑖𝑘𝜙^𝐼,at𝑟=𝑟₀≫0, (32) where 𝑟=𝑥²+𝑦² and 𝜙^𝐼=𝑎_𝐼𝑒^𝑖𝑆 is the pseudo free surface velocity potential of incident wave in which 𝑆 denotes its phase function, and 𝑎_𝐼 denotes wave amplitude. For waves propagating toward the +𝑥 direction, the phase 𝑆 of the incoming wave is thus 𝑖𝑘𝑥. The no-flux boundary condition given at the cylinder surface is 𝜕𝜙𝜕𝑟=0,at𝑟=𝑎, (33) where 𝑎 is the radius of cylinder.

The analytical solution [56] of the problem governed by (29) subjected with the boundary condition given as (31)-(32) is 𝜙=𝑎_𝐼^∞_𝑛=0𝜀_𝑛𝑖^𝑛𝐽_𝑛(𝑘𝑟)𝐻_𝑛()×𝐽𝑘𝑟^′𝑛(𝑘𝑎)𝐻^′_𝑛(𝑘𝑎)cos𝑛𝜃−𝜃_𝐼, (34) where 𝜀=1,if𝑛=02,if𝑛>0 (35) and 𝜃_𝐼 is the incident angle of waves, and 𝐽_𝑛 and 𝐻_𝑛 are the Bessel and Hankel functions of the first kind of order 𝑛. The superscript of the symbol of prime indicates their derivatives, respectively. The radial and tangential velocities are obtained by the differentiation of the analytical solution with respect to 𝑟 and 𝜃, respectively: 𝑈_𝑟=𝜕𝜙=𝜕𝑟𝑔𝑎_𝐼𝜔×^∞_𝑛=0𝜀_𝑛𝑖^𝑛−1𝑘𝐽^′_𝑛(𝑘𝑟)−𝑘𝐻^′_𝑛(𝑘𝑟)𝑘𝐽^′_𝑛(𝑘𝑎)𝑘𝐻^′(𝑘𝑎)×cos𝑛𝜃−𝜃_𝐼𝑈_𝜃=1𝑟𝜕𝜙=𝜕𝜃𝑔𝑎_𝐼×𝜔𝑟_∞𝑛=0𝜀_𝑛𝑖^𝑛−1−𝐽_𝑛(𝑘𝑟)+𝑘𝐻^′_𝑛(𝑘𝑟)𝑘𝐽^′_𝑛(𝑘𝑎)𝐻^′_𝑛(𝑘𝑎)×𝑛sin𝑛𝜃−𝜃_𝐼. (36)

A case of water depth of 10 m and incident waves with period 𝑇=20 s, which corresponds to wave number 𝑘=0.0308 m⁻¹, is chosen for demonstration. The radius of the cylinder is chosen as 30 m while the open outer boundary is located at 𝑟₀=295 m. The nodal arrangement for the computation is shown in Figure 9. Around the first circle (i.e., the cylinder boundary), 30 nodes are placed. The number of nodes in each outer circle is increased by 6 while the radius of the circle is increased by 6.28 m. Totally 43 circles (including cylinder boundary) are arranged so there are 6,708 nodes. Contours of numerical results of amplitude, phase, and magnitude of velocity within part of the computational domain are compared with analytical solutions as shown in Figures 10, 11, and 12. Numerical results of amplitude, phase, and magnitude of velocity around the cylinder are compared with the analytical solutions in Figure 13. Very good agreement in these comparisons implies that the MFPM is applicable not only to partial differential equation of real variables, but also to those of complex variables.

Figure 9 

Nodal arrangement for MFPM computations of wave scattering.

(a)

(b)

Figure 10 

Contours of equal-amplitude of waves scattered by a cylinder in constant water depth, 𝑇=20 s. (a) MFPM results; (b) analytical solutions.

(a)

(b)

Figure 11 

Contours of equal-phase of waves scattered by a cylinder on constant water depth, 𝑇=20 s. (a) MFPM results; (b) analytical solutions.

(a)

(b)

Figure 12 

Contours of magnitude of velocity of waves scattered by a cylinder on constant water depth, 𝑇=20 s. (a) MFPM results; (b) analytical solutions.

(a)

(b)

(c)

Figure 13 

Comparisons of MFPM results and analytical solutions of amplitude, phase, and magnitude of velocity around at 𝑟=𝑎.

4. Further Applications

After the correctness and applicability of MFPM are verified in the previous section, applications of MFPM to several problems are further demonstrated in this section.

4.1. Numerical Generation of Orthogonal Grids

Though meshless numerical methods are becoming more and more popular in recent decades, grid-based numerical methods, which have been fully developed, are still indispensable. Transformation of coordinates from Cartesian system (𝑥-𝑦 coordinate) to a curvilinear system (𝜉-𝜂 coordinate) plays an important role in the numerical simulations using grid-based methods. The most desirable transformation may be the one that could generate orthogonal grid lines on the boundaries as well as inside the domain, so the discretization of the governing equation is simplified and thus the discretization error in solving the governing equation could be minimized [57]. The orthogonal coordinate transformation can be accomplished by solving two pairs of Laplace equations and subjected to the Cauchy-Riemann boundary conditions [58].

One application of MFPM is to employ it to solve the Laplace equations of numerically generate 2-D orthogonal grids. Based on Cauchy-Riemann conditions, the orthogonal transformation from the physical domain (𝑥-𝑦 coordinate) to the computational domain (𝜉-𝜂 coordinate) can be reduced to a problem governed by two Laplace equations:𝜕²𝜉𝜕𝑥²+𝜕²𝜉𝜕𝑦²=0;𝜕²𝜂𝜕𝑥²+𝜕²𝜂𝜕𝑦²=0 (37a)(37b)

with the boundary conditions⇀−𝑛⋅∇𝜂=0on𝜉=0and𝜉=𝜉₀⇀−𝑛⋅∇𝜉=0on𝜂=0and𝜂=𝜂₀, (38a)(38b)

where ⇀−𝑛 is the outward normal vector of the boundaries. After (37a) and (37b) are solved separately, the obtained 𝑥-𝑦 to 𝜉-𝜂 one-to-one corresponding relationship can be used as the boundary conditions for seeking the inverse transformation which is also governed by two Laplace equations:𝜕²𝑥𝜕𝜉²+𝜕²𝑥𝜕𝜂²=0;𝜕²𝑦𝜕𝜉²+𝜕²𝑦𝜕𝜂²=0. (39a)(39b)

The relations between the partial derivatives are listed below:𝜕𝜉=1𝜕𝑥𝐽𝜕𝑦;𝜕𝜂𝜕𝜉=𝜕𝑦−1𝐽𝜕𝑥𝜕𝜂𝜕𝜂=𝜕𝑥−1𝐽𝜕𝑦;𝜕𝜉𝜕𝜂=1𝜕𝑦𝐽𝜕𝑥,𝜕𝜉 (40a)(40b)(40c)(40d)

where 𝐽 is the Jacobian and is defined as 𝐽=𝜕𝑥𝜕𝜉𝜕𝑦−𝜕𝜂𝜕𝑥𝜕𝜂𝜕𝑦.𝜕𝜉 (41)

A sketch of the transformation is illustrated in Figure 14. A detailed derivation of these equations could be found in [58].

(a)

(b)

Figure 14 

Sketch of a transformation between two orthogonal coordinates from the physical domain to the transformed domain.

Conventionally the four corners of the irregular domain have to be right angles and the irregular domain is denominated as a “hyper-rectangle” [59–62], which consists of four right corners and four smooth curves. A method of conformal mapping in the complex-value domain is integrated to transform any irregular geometry into a hyper-rectangle has been developed in [59–62], except the case which contains interior corners of zero-degree angle. The grid generation in a region shown in Figure 15 was a bottleneck in previous works [59–62], because of a local complex transformation failed at points with zero-degree internal angle. The orthogonal grid generation in this irregular region with zero-degree internal corners is presented here for completeness. The values of 𝜉₀ and 𝜂₀ are all set to be 1 because of the symmetry. For forward transformation from 𝑥-𝑦 to 𝜉-𝜂, the collocation points inside the domain are regularly distributed with a 45-degree inclination, as shown in Figure 16. The spacing of collocation points on the boundary is 1/80 of the arc length. Collocation points are selected in the way that the nearest distance between any two collocation points is greater than 0.007. So some collocation points around the zero-degree internal corners are skipped.

Figure 15 

The domain of an irregular area with zero-degree internal corners.

Figure 16 

Arrangement of the collocation points in the up-right quarter of the irregular area with corners with zero-degree angle.

There is no exact solution for this problem, but one can still verify the correctness of the transformation by checking the orthogonality. Because the orthogonality is guaranteed mathematically, the following relative error criteria are used to show the correctness of the results:ARE_𝑂=∑^𝑁_𝑗=1∇𝜉⋅∇𝜂²_{𝑥−⇀𝑥=−⇀𝑗}∑^𝑁_𝑗=1∇𝜉²∇𝜂²_{𝑥−⇀𝑥=−⇀𝑗}MRE_𝑂=max∇𝜉⋅∇𝜂²_{𝑥−⇀𝑥=−⇀𝑗},𝑗=1∼𝑁1𝑁∑^𝑁_𝑗=1∇𝜉²∇𝜂²_{𝑥−⇀𝑥=−⇀𝑗}. (42a)(42b)

The subscript “𝑂” denotes “orthogonal.” The smaller the values of these two relative error criteria are, the higher the correctness of the numerical results will be. MFPM numerical result shows that the averaged relative error and maximum relative error calculated by the collocation points are 7.60×10⁻⁵ and 7.53×10⁻⁴, respectively.

Choosing Δ𝜉=Δ𝜂=0.05, generated grid lines are plotted in Figure 17. Each line of one set is observed very well perpendicular to all the lines of the other set. To illustrate more details of the grid lines around a zero-degree internal corner, grid lines with Δ𝜉=Δ𝜂=0.01 are plotted in Figure 18. The density of the grid nodes is still very coarse near the corner. The point (𝜉,𝜂)=(0.99,0.99) is far away from the corner, compared with the distances from other neighboring 𝜉-𝜂 grid nodes to it. This implies that the scaling factors 𝜕𝜉/𝜕𝑥, 𝜕𝜉/𝜕𝑦, 𝜕𝜂/𝜕𝑥, and 𝜕𝜂/𝜕𝑦 approach to zero at the corners. One could find more details of this application and other applications of MFPM to such kind of problems in [47].

Figure 17 

Orthogonal grid lines in an irregular area with internal corners of zero-degree angle, Δ𝜉=𝜉0/20,  Δ𝜂=𝜂0/20,  𝜉0=1, and 𝜂0=1.

Figure 18 

Orthogonal grid lines around a zero-degree internal corner, Δ𝜉=𝜉0/100, Δ𝜂=𝜂0/100,  𝜉0=1, and 𝜂0=1.

4.2. Frequency-Domain Water Wave Problems

4.2.1. Combined Refraction and Diffraction of Water Waves

Various numerical methods have been developed to study the combined water-wave refraction and diffraction problems in the last decade. Notably, Tsay and Liu [63] based on the linear long-wave theory developed a hybrid finite element model, which can be used for studying harbor oscillations and calculating wave forces and moments on offshore structures. Following the same idea of using a hybrid method, Yue et al. [64] developed a three-dimensional finite element model which is applicable for a wider range of wave lengths. However, the three-dimensional model requires huge storage and computational time when the incident wave length becomes short.

Alternatively, the mild-slope equation derived by Berkhoff [65] and Smith and Sprinks [66] has gained popularity as a base for studying combined refraction and diffraction of water waves, since it reduces a three-dimensional problem to a two-dimensional one. Furthermore, this equation is also applicable to a wide range of wave frequencies. Of course, the drawbacks for using the mild-slope equation are that water depth must be assumed to be slowly varying in space and the evanescent modes are ignored. Nevertheless, Bettess and Zienkiewicz [67], Houston [68, 69], and Tsay and Liu [63] developed by hybrid methods based on the mild-slope equation. Excellent results were obtained for problems concerning wave scattering by islands as well as harbor oscillations. In Bettess and Zienkiewicz’s [67] work, infinite elements were used in the region away from the scatters. Houston [68, 69] showed that his model can generate accurate solutions by using very fine elements. Neither model can, however, calculate wave forces on a semi-submerged offshore structure. Tsay and Liu [63] used a hybrid finite element model originally developed by Chen and Mei [70] for shallow water wave problems. The physical domain is divided into two regions as a finite inner region of variable depth and infinite outer region with constant or deep water depth. The infinite outer region is represented by analytical solutions with undetermined coefficients which satisfy the governing equation and the radiation condition on the outer boundary. Conventional discretization of finite element is used in the variable depth region.

Following the procedures that Berkhoff [65] used, the mild-slope equation for water waves could be briefly shown as ∇⋅𝐶𝐶_𝑔∇𝜂+𝑘²𝐶𝐶_𝑔𝜂=0, (43) where 𝐶=𝜔/𝑘 and 𝐶_𝑔=(1/2)(1+2𝑘ℎ/sinh2𝑘ℎ)𝐶 represent wave celerity and group celerity, respectively; 𝑘 represents the wave number. In case the water depth is a constant, (43) can be reduced to (29).

Generally, the treatment of the boundary far from the cylinder is the same as in (32). When there are some man-made structures such as breakwaters and wharfs in harbors, a no-fluxing boundary at such solid boundaries is required and can be expressed as 𝜕𝜂𝜕𝑛=0. (44)

In this section, the MFPM is implemented for the solution of the mild-slope equation for combined refraction and diffraction when monochromatic waves encountered with scatters, due to change of water depth or artificial structures. In present numerical model, the velocity of fluid particles at the grid points can be found easily as part of the coefficients resulted from the coefficients of collocation at each computational node. To verify present numerical model, analytical solutions of a benchmark problem in shallow water waves scattered by a circular island on a paraboloidal shoal by Homma [71] have been applied. It should be noted that the unknown variables of the numerical and analytic solutions are complex numbers.

As shown in Figure 19, the depth between the shoreline and outer boundary of the shoal is described as ℎ=ℎ₀𝑟𝑟₀²,𝑎≤𝑟≤𝑟₀, (45) where 𝑟 is the distance measured from the center of the island, 𝑎 is radius of the circular island, 𝑟₀ is the radius of paraboloidal shoal, and ℎ₀ is the constant depth outside the paraboloidal shoal.

Figure 19 

Definition of a circular island on a paraboloidal shoal.

Analytical solutions of propagating shallow water wave responses around this circular island for 𝑎≤𝑟≤𝑟₀, derived by Homma [71], are presented as the following for completeness: 𝜂=𝑎₀^∞_𝑛=0𝑘𝜀_𝑛𝑖^2𝑛𝐽^′_𝑛𝑘𝑟₀𝐻_𝑛𝑘𝑟₀−𝐽_𝑛𝑘𝑟₀𝐻^′_𝑛𝑙𝑟₀×𝛼_𝑛+1𝛼_𝑛𝑎−1^−2𝛼𝑛𝑟^𝛼𝑛−1+𝑟^{−𝛼𝑛−1}+𝛼cos𝑛𝜃_𝑛+1𝑎^−2𝛼𝑛𝑟^𝛼𝑛₀⁻²−𝑟^−𝛼𝑛₀⁻²𝐻_𝑛𝑘𝑟₀−𝑘𝐻^′_𝑛𝑘𝑟₀𝛼_𝑛+1𝛼_𝑛𝑎−1^−2𝛼𝑛𝑟^𝛼𝑛₀⁻¹+𝑟^−𝛼𝑛₀⁻¹. (46)𝜃_𝐼 is the incident angle of propagating waves and 𝐽_𝑛 and 𝐻_𝑛 are the Bessel and Hankel functions of the first kind of order 𝑛, as in the fourth benchmark problem. The superscript of the symbol of prime indicates their derivatives, respectively.

Differentiating (47) with respect to radial (𝑟) and angular (𝜃) variables, respectively, the analytical solutions of radial and tangential velocities in the inner region of 𝑎≤𝑟≤𝑟₀ can be obtained analytically with any difficulty.

A case for a circular island with radius, 𝑎=10 km, as shown in Figure 19, is considered. The water depth is constant, ℎ₀=4 km, at distance greater than 𝑟₀=30 km. The amplitude and period of incident wave are 1 m and 720 s, respectively. The radius of the local radiation boundary is 275 km which is about 2 times of wave length from the paraboloidal shoal boundary. 120 nodes are arranged at the circumference of the circular island. For the next circle, the spaces of radial and tangential directions of nodes increased by 5%. After 55 circles of nodes, there are 4 circles increased by 6 more nodes for each one. In the case there are totally 9,910 nodes and 59 circles arranged. The numerical results of contours of amplitude, phase, and magnitude of velocity are compared with analytical solutions as shown in Figures 20, 21, and 22. Comparisons of numerical of results of amplitude, phase, and magnitude of velocity at 𝑡=0 with the analytical solutions around the margin of cylinder are shown in Figure 23. Good agreements in these comparisons imply the applicability of MFPM to complex variable problems and simulate water-wave scattering problems. It should be noted that the velocity components in 𝑥- and 𝑦-direction, respectively, are part of the solutions in the coefficients of local polynomial collocations. Therefore, present MFPM is more efficient than other numerical methods using linear shape functions.

(a)

(b)

Figure 20 

Contours of wave amplitude lines around a circular island on a paraboloidal shoal, 𝑇=720 s. (a) Numerical result; (b) analytic results.

(a)

(b)

Figure 21 

Contours of wave phase lines around a circular island on a paraboloidal shoal, 𝑇=720 s, at 𝑡=0. (a) Numerical results; (b) analytic results.

(a)

(b)

Figure 22 

Contours of magnitude of velocity around a circular island on a paraboloidal shoal, 𝑇=720 s, 𝑡=0. (a) Numerical results; (b) analytic results.

(a)

(b)

(c)

Figure 23 

Numerical and analytic results of amplitude, phase, and magnitude of velocity around at 𝑟=𝑎.

4.2.2. Harbor Resonance

Although present meshless methods have advantages over element- or mesh-based numerical methods to avoid tedious preprocessing of the mesh connectivity information in most the numerical computations, there are problems in engineering and science consisting of different regions in its physical behaviors divided by a very thin barrier or crack, such as a breakwater in water-wave diffraction, a crack in thin plates, and a cutoff in groundwater seepage. In present MFPM, the local collocation is performed by searching some closest nodes near the base point. When physical behaviors in separated regions by a thin barrier, it requires to identify the nodes searched by relative distance are not from the different regions with different physical behaviors. Therefore, it is required to divide the whole domain into subdomains of different physical behaviors. It is referred to as regional connectivity to complement the searching procedures of the local nodes for appropriate collocations [49]. In this section, applications of present MFPM will be demonstrated in coastal engineering problem for harbor resonance with breakwaters.

To study resonance phenomena in a harbor basin, Lee [72] first developed deep basin resonance theory (arbitrary-shape harbor theory) and used separation of variables method to obtain the basin resonance for a circular harbor and a narrow channel. Compared with experimental measurements, the theoretical results of the narrow channel showed good agreement. Many numerical methods have been developed based on different linear, small-amplitude wave theory and obtained very good results [63, 65–67, 70, 73].

However, there exist two drawbacks of these kinds of numerical models. First, a tedious preprocessing procedure is required to provide the connectivity information of the computational nodes and meshes in FEM. Secondly, numerical models with linear shape functions fall short of providing velocity fields of the water-wave motions.

For water-wave combined refraction and diffraction problems with abrupt change of geometry in computational domains, such as existence of breakwater(s), a concept of subdomains of regional connectivity has been established to guarantee application of local collocation that is performed within appropriate subdomains in applications of MFPM. The governing equation and boundary conditions can be referred to as the same of (43)–(45).

Protruding breakwaters or sharp geometrical corners exist in most harbors and may behave like a branch-cut line to divide the water area into subdomains with different characters such as up-wave and down-wave sides, respectively. In the collocation method of MFPM, the approximate polynomial function is established by identifying mostly and closely nearby computational grids in relative distances to a base point. When a protruding breakwater exists, special care must be taken to avoid identifying grid points from different subdomains in this searching process to evaluate coefficients of local polynomials in collocation approximation. Therefore, the computational domain is divided into subdomains and regional connectivity is specified for further local polynomial approximation at each node. In diffraction of water-wave problems, a computational domain can be divided into three subdomains which are the region of ocean, the entrance region of the harbor, and the basin region, shown as in Figure 24 for circular and rectangular harbor, respectively. In searching mostly closed nodes, only nodes in its own and adjacent subdomains indicated by regional connectivity information are included. For example, in Figure 24, when collocation approximation is applied to subdomain I, only nodes in subdomains I and II are considered. For subdomain III, only nodes in subdomains II and III are taken into account to search some closest modes. For subdomain II, all of nodes in its own and adjacent subdomains, namely, II, I, and III, are considered.

Figure 24 

Regional connectivity of subdomains for closest nodal searching.

A circular harbor with protruding breakwaters into the ocean side is considered, as shown in Figure 25. In this benchmark problem, when monochromatic incident waves propagating into a circular harbor with protruding breakwaters, Mei and Petroni [74] had obtained analytical solutions for wave responses within the region of harbor, as shown in (47)–(49) as follows: 𝜂_𝐻=𝑎₀𝑁𝐷^∞_𝑛=0𝜀_𝑛𝐽2𝜋⁰(𝑛𝛿)𝑘𝑎𝐽^′_𝑛()𝐽𝑘𝑎^𝑛(𝑘𝑟)cos𝑛𝜃−𝜃_𝐼, (47) where 𝜂_𝐻 is water free-surface displacement within the harbor, 𝑎₀ is incident wave amplitude, 𝛿 is half of the harbor entrance opening angle, 𝜃₀ is the angle from the positive 𝑥-coordinate to the centerline of entrance opening, 𝜃_𝐼 is direction of incident waves with respect to positive 𝑥-axis, 𝑎 is radius of the circular harbor, 𝑘 is incident wave number, and 𝑁 is weighted average of the forcing pressure: 𝑁=^∞_𝑛=02𝑖𝜋𝑘𝑎2𝜀_𝑛𝐽₀(𝑛𝛿)𝐻^′_𝑛()×𝑘𝑎cos𝑛𝜋2cos𝑛𝜃_𝐼cos𝜃₀+𝑖sin𝑛𝜋2sin𝑛𝜃_𝐼sin𝑛𝜃₀, (48)

Figure 25 

Definition of a circular harbor with protruding breakwaters.

and 𝐷 is term of the coefficient: 𝐷=^∞_𝑛=0𝜀_𝑛𝐽2𝜋^𝑛(𝑘𝑎)𝐽²₀(𝑛𝛿)𝑘𝑎𝐽^′_𝑛()−𝑘𝑎_∞𝑛=0𝜀_𝑛𝐻2𝜋^𝑛(𝑘𝑎)𝐽²₀(𝑛𝛿)𝑘𝑎𝐻^′_𝑛(𝑘𝑎)1+cos𝑛𝜋cos2𝑛𝜃₀. (49)

In this problem, verifications of numerical calculations will be presented by comparisons with analytical results of Mei & Petroni [74]. Two protruding breakwaters into the ocean side clearly separate wave field into two distinct regions. A bridging subdomain at the harbor entrance connects numerical computational information between ocean and harbor sides in present local collocation approach. There are total 6,336 computational nodes used for ocean side (subdomain I), 10 nodes at the harbor mouth (subdomain II), and 449 nodes for the harbor region (subdomain III). The total computational nodes are 6,795. The computational time consumes about 5 minutes of CPU time for one frequency case with computer facility of the personal computer with the Intel(R) Core(TM)2 duo CPU E7300 (2.6 GHz 2.6 GHz) and Microsoft Windows 32 bits system.

A case of this application with incident wave direction 𝜃_𝐼 of 225° is presented. Present numerical results for an idealized circular harbor of radius 𝑎=0.75 ft with harbor mouth 60° opening (2𝛿).

Response curves at 𝑟=0.70 ft, 𝜃=0^∘, are shown in Figure 26. It is observed that there is very good agreement between present MFPM results and analytical solutions [74]. In the case of these obliquely incident waves, one major resonant wavenumber happens at 𝑘𝑎=0.5 and the amplification factor, 𝑅, is equal to 1.46. At the resonant wavenumber, the corresponding wavelength and period are equal to 2.87 m and 1.37 s, respectively. Contours of wave amplitude, Figure 27, clearly show different features of wave fields in the up- and down-wave sides, respectively, due to the sheltering effects of the breakwaters. The maximum wave amplitude occurs at the most inner end of the circular harbor. Velocity distributions of the harbor at 𝑘𝑎=0.5 indicate the flow directions of the obliquely incident waves at two different time instants, Figure 28. The maximum velocity happens again at the harbor entrance.

Figure 26 

Response curve at 𝑟=0.7 ft, 𝜃=0∘ of a circular harbor of radius, 𝑎=0.75 ft, 60∘ opening for 𝜃𝐼=225∘.

Figure 27 

Contours of wave amplitude of a circular harbor of radius, 𝑎=0.75 ft, 60∘ opening, with protruding breakwaters at 𝑘𝑎=0.5 for 𝜃𝐼=225∘.

(a)

(b)

Figure 28 

Velocity vectors of a circular harbor of radius, 𝑎=0.75 ft, 60° opening, with protruding breakwaters at 𝑘𝑎=0.5. (a) 𝑡=0, (b) 𝑡=𝑇/2 for 𝜃𝐼=225∘.

4.2.3. Kelvin Waves

For waves with very long periods, such as tidal waves, the Coriolis effects cannot be ignored. Based on the mild-slope equation, an extended mild-slope equation has been formulated to take the effects of Coriolis effects into account, when wave period becomes very long or the domain of concerns is large [73]: ∇⋅𝑃𝐺∇𝜂+𝜆²𝑓𝑃𝐺1−²𝜔²𝜂=0, (50) where 𝑓 represents the Coriolis parameter and is rotating in a counter-clockwise direction in the northern hemisphere; 𝜔 represents frequency of incident wave; 𝜆 is the eigenvalue of the characteristic equation, 𝜔²=𝑔𝜆tanh𝜆ℎ, which may be real or imaginary values; 𝑃 and 𝐺 represent, respectively, the wave and group celerity as 𝜔𝑃=𝜆1𝐺=21+2𝜆ℎsinh2𝜆ℎ𝑃, (51) and 𝑘 is the wave number which could be expressed as 𝑘²=𝜆²𝑓1−²𝜔². (52)

For real value of 𝜆, 𝑘 is a real value when 𝑓<𝜔, there exists a propagating mode of waves, namely, Kelvin waves. For waves with shorter periods, it reduces the traditional mild-slope equation [65, 66] when Coriolis parameter is negligible.

Pedlosky [75] had studied the case of Kelvin waves by using the linearized momentum equations expressed as 𝜕𝑢𝜕𝑡−𝑓𝑣=−𝑔𝜕𝜍𝜕𝑥𝜕𝑣𝜕𝑡+𝑓𝑢=−𝑔𝜕𝜍𝜕𝑦𝜕𝑤𝜕𝑡=−𝑔𝜕𝜍.𝜕𝑧 (53)

In this section, Coriolis effects on long waves propagating in an infinite channel are investigated. The channel has two parallel side vertical walls of width 𝐿, which is shown in Figure 29. Water depth is 𝐻₀ and waves propagating from the left to the right.

Figure 29 

An infinite channel of width 𝐿, with Coriolis parameter 𝑓.

According to Pedlosky [75], the analytic solution in Kelvin wave could be expressed as 𝜂=Re𝜂_𝑜𝑒^{−𝑓𝑦/𝐶𝑜}𝑒^𝑖𝜆𝑥𝜂𝑢=⁰𝐻₀𝐶₀𝑒^{−𝑓𝑦/𝐶0}cos𝑘𝑥±𝜔𝑡+𝜑, (54) where 𝐶₀=(𝑔𝐻₀)^1/2; 𝐻₀ is constant depth of the channel; 𝑔 is the gravity; 𝜂₀ is a constant amplitude of incident wave; 𝜆 is the wave number without effect of Coriolis force.

For numerical computations of MFPM, both radiation condition at the far-field open boundary and no-flux condition on the solid boundaries need to be modified accordingly. Based on the wave modulation of a vertical distribution with 𝐹(𝑧)=cosh𝜆(𝑧+ℎ)/cosh𝜆ℎ, the fluid displacement, 𝜍, can be expressed as 𝜍𝑥,𝑦,𝑧;𝑡=𝜂𝑥,𝑦𝐹(𝑧)𝑒^{−𝑖𝜔𝑡}𝜂𝑥,𝑦=𝐴𝑥,𝑦𝑒^{𝑖∫𝑘𝑛𝑥⋅𝑑−⇀}=𝐴𝑥,𝑦𝑒^{𝑖∫(𝑘𝑥𝑑𝑥+𝑘𝑦𝑑𝑦)}, (55) where 𝐴 is amplitude of free surface displacement, 𝜂; 𝑘 is wave number; 𝑘_𝑥 and 𝑘_𝑦 are the components of 𝑘 in 𝑥- and 𝑦-direction, respectively.

Assuming that the flow velocities associated with wave action, (𝑢,𝑣), also have the same vertical distribution and with (𝑈,𝑉) as the components at the free surface: (𝑢,𝑣)=(𝑈,𝑉)𝐹(𝑧). (56)

From (53) and (55)-(56), the velocity components formula with Coriolis effects can be derived and expressed as 𝑈=𝑔𝑓𝜔²−𝑓²𝜕𝜂−𝜕𝑦𝑖𝑔𝜔𝜔²−𝑓²𝜕𝜂𝜕𝑥𝑉=−𝑔𝑓𝜔²−𝑓²𝜕𝜂−𝜕𝑥𝑖𝑔𝜔𝜔²−𝑓²𝜕𝜂.𝜕𝑦 (57)

Then from the no-flux condition, 𝑈⋅𝑛=0, the solid boundary condition could be expressed as 𝜕𝜂𝜕𝑛=±𝑖𝑓𝜔𝜕𝜂,𝜕𝑠 (58) where 𝑛 and 𝑠 represent normal and tangent direction on the solid boundary, respectively. It is noted that the modified no-flux condition, (59), reduces to a traditional one when Coriolis parameter is equal to zero.

For practical purposes, outgoing propagation of scattered waves satisfy the Sommerfeld [55] radiation condition as (31) can be approximated as a local one as 𝜕𝜂^𝑆𝜕𝑛−𝑖𝜆𝜂^𝑆=0. (59)

In order to illustrate the Coriolis effects of waves propagating in a channel with width 𝐿=60 km and water depth 𝐻₀=3 m have been used. A wave period 𝑇=46800 sec, amplitude = 1.0 m, and Coriolis parameter 𝑓=0.0000935638 sec⁻¹ (equivalently at latitude = 40°) have been taken. The wave in the channel is propagating from the left to right. The result of analytical solution is shown as in Figure 30. The distribution of wave amplitude decreases exponentially from right to left for an observer looking in the same direction of wave propagation. The rate of decrease of the amplitude from right to left is proportional to 𝑐/𝑓, the Rossby radius of deformation. The distribution of the wave velocity is parallel to the walls of the channel and also proportional to 𝑐/𝑓.

(a)

(b)

Figure 30 

The results of analytical solutions: (a) distribution of wave amplitude; (b) distribution of wave velocity.

For numerical computations, the incident wave amplitude is chosen to be 𝜂_𝑜𝑒^{−𝑓𝑦/𝐶𝑜} and propagates from left to right. The length of the channel of 1,000 km was used for numerical computations. The local radiation boundary condition as (59) was used on the other end of the channel for the incident wave and scatter wave. By applying this numerical MFPM model, the velocity at each node can be also obtained simultaneously with the free surface displacement. Present numerical results of amplitude and velocity, Figure 31, show very good agreement with analytical solutions. The amplitude and velocity of 𝑥-direction between the numerical and analytical results at cross section 𝑥 = 10,000 m are further compared as shown in Figure 32. Very good agreement between numerical and analytical results is observed. The velocity of 𝑦-direction is equal to 0 with a negligible error. Very good agreement between numerical results and analytical solutions is observed.

(a)

(b)

Figure 31 

The results of numerical model: (a) distribution of wave amplitude; (b) distribution of wave velocity.

(a)

(b)

(c)

Figure 32 

Comparisons of amplitude and velocity at 𝑥-diretion between numerical and analytical results at the cross section 𝑥 = 10,000 m.

4.3. Time-Domain Water Wave Problems

For several decades, free-surface water wave problems have been studied as potential-flow problems governed by the Laplace equation with nonlinear free surface boundary conditions. If the motions of the free-surface are periodic, the problems can be treated as frequency-domain problems as illustrated above. But if one wants to obtain more detail in the transient process of the wave motion, time-domain water wave problems are recommended to be considered.

Most of these studies have been carried out with boundary-element methods (BEM), subjected to a mixed Eulerian-Lagrangian (MEL) time marching approach [76–81]. The governing equation is the Laplace equation ∇²𝜙=0, (60) where 𝜙(⇀−𝑥,𝑡) is the velocity potential with ⇀−𝑥 being the position vector and 𝑡 being the time. The velocity vector is defined as ⇀−𝑣=∇𝜙 ((𝜕𝜙/𝜕𝑥)𝑗⇀−𝑖+(𝜕𝜙/𝜕𝑦)⇀− for 2-D problems and (𝜕𝜙/𝜕𝑥)𝑘⇀−𝑖+(𝜕𝜙/𝜕𝑦)⇀−𝑗+(𝜕𝜙/𝜕𝑧)⇀− for 3-D problems). The kinematic free-surface boundary condition (KFSBC) and dynamic free-surface boundary condition (DFSBC) can be expressed as 𝑑𝑥∇𝜙=⇀−,𝑑𝑡⇀−𝑥∈theentiredomain,includingthefreesurface,𝑑𝜙1𝑑𝑡=−𝑔𝑧+2∇𝜙⋅∇𝜙,⇀−𝑥∈thefreesurface, (61) where 𝑔 is the gravitational acceleration, 𝑧 is the vertical coordinate, and 𝐶(𝑡) is Bernoulli’s constant, which can be set to zero for a quiescent ambiance. Both the KFSBC and the DFSBC are stated in the Lagrangian frame. The boundary condition at the water-body interface is the no-flux boundary condition, which can be expressed as ⇀−𝑛⋅∇𝜙=⇀−𝑛⋅⇀−𝑣_𝑏, (62) where ⇀−𝑛 is the unit normal vector outward from the body surface, and ⇀−𝑣_𝑏 is the velocity of the body surface. The right-hand side of this equation equals zero for a fixed bottom in a finite depth and will be given by the motion of the solid boundary.

For solving this time dependent problem, the time-domain firstly has to be discretized. At each time step, the Laplace equation needs to be solved once to obtain the velocity potential for the entire domain thus to get the velocity. Boundary conditions are updated by the motion of the solid boundary and prediction from the free-surface boundary, which involves nonlinear terms. One could use the a mixed Eulerian-Lagrangian (MEL) time marching approach in which the 4th order Runge-Kutta time integration is used, or choose a simpler but satisfactory time marching scheme proposed in [21]. Adopting the time marching scheme proposed in [21], the DFSBC (66) can be formulated as 𝜙_𝑥^−⇀_𝑥=^−⇀_𝑗^(𝑛)=𝜙_𝑥^−⇀_𝑥=^−⇀_𝑗^(𝑛−2)1+2Δ𝑡−𝑔𝑧+2∇𝜙⋅∇𝜙_{𝑥−⇀𝑥=−⇀𝑗}^(𝑛−1), (63) where ⇀−𝑥_𝑗 denotes the position of the 𝑗th node and this equation is only valid in case the node is on the free surface. In this formulation, the required data on the right-hand side for seeking the velocity potential in the entire domain at the 𝑛th time step, including the nonlinear terms, are already known. Generally speaking, this explicit time marching scheme may be regarded as a leap-frog approach. What one needs to do first is just to determine the position of each traced “particle,” ⇀−𝑥_𝑗^(𝑛). It was proposed to use the second-order finite difference scheme in the time-domain, also a leap-frog approach, to determine the position of ⇀−𝑥_𝑗^(𝑛): ⇀−𝑥_𝑗^(𝑛)=⇀−𝑥_𝑗^(𝑛−2)+2Δ𝑡∇𝜙_𝑥^−⇀_𝑥=^−⇀_𝑗^(𝑛−1). (64)

In the case of strong nonlinearity, the free-surface particle trajectory ⇀−𝑥_𝑗^(𝑛) should be modified by applying the Crank-Nicolson formula for a better numerical stability: ⇀−𝑥_𝑗^(𝑛)=⇀−𝑥_𝑗^(𝑛−1)+Δ𝑡2∇𝜙_𝑥^−⇀_𝑥=^−⇀_𝑗^(𝑛)+∇𝜙_𝑥^−⇀_𝑥=^−⇀_𝑗^(𝑛−1). (65)

The above equation is applicable after the velocity potential for the entire domain is solved. Note that there is no need to solve the Laplace equation again because there is barely difference between the free-surface velocity potential at ⇀−𝑥_𝑗^(𝑛) predicted by using (64) and that predicted by using (65). The accuracy of this time marching scheme has been verified in [21, 50–53].

4.3.1. Solitary Wave Generation by a Piston-Type Wave Maker

The first application of the MFPM to transient water wave problem is the simulation of the solitary wave generation by a piston-type wave maker [52]. In [52] a modification of Goring’s solitary wave generation method [82] was proposed. The motion of the solid boundary in the model (i.e., (66)) is described as 𝑑𝜉=𝑑𝑡𝑢(𝜉,𝑡)=𝐶𝜂_𝑥=𝜉ℎ+𝜂_𝑥=𝜉, (66) where 𝜉(𝑡) is the position of the wave paddle at time 𝑡, while 𝑡 is elapsed time since the start of the motion, 𝐶 is the wave celerity, 𝜂 is the free surface displacement, and ℎ is the still water depth. In [82] it was proposed to apply Boussinesq’s solitary wave solution [83] into (66). That is, 𝜂=𝐻𝑆², (67) where 𝐻 is the wave height, and 𝑆=cosh⁻¹[𝐾𝑋] (68) in which 𝐾 is the boundary outskirt decay coefficient, and 𝑋=𝜉−𝐶𝑡−𝑥₀, (69) where 𝑥₀ is the initial position of the wave crest. Because the solitary wave is generated in a quiescent water flume, 𝑥₀ must be chosen as a negative value, which means the wave crest is initially out of the domain. Its value is chosen by considering the length of the wave. The boundary outskirt decay coefficient of Boussinesq’s solitary wave solution is determined as 𝐾=3𝐻4ℎ³, (70) where the wave celerity is determined as 𝐶=𝑔(ℎ+𝐻). (71)

Figure 33 shows the snapshots of fluid particles during the wave generating process. In the case, the water depth is 20 cm and Goring’s solitary wave generation method is used. The input in (71) is 7.4 cm. Details about the setup of this case and the verification of the model could be found in [52]. How the height of the wave crest decays during the wave propagation could also be found in [52].

(a)

(b)

(c)

(d)

Figure 33 

Snapshots of fluid particles during wave generating process, for the case of water depth with 20 cm and input wave height of 7.4 cm.

In [52] it was suggested that one should apply the 9th order solitary wave solution [84] into (70) so the generated solitary wave could propagate for a longer distance without the wave crest decay. In [84], the free surface displacement can be described as 𝐻𝜂=ℎℎ𝜂₁+𝐻ℎ²𝜂₂𝐻+⋯+ℎ⁹𝜂₉. (72)

The formulae to determine the values of 𝜂₁ to 𝜂₉ could be found in [84]. In an analogous way, the boundary outskirt decay coefficient and the wave celerity are determined as 𝐾=3𝐻4ℎ³𝐻1+ℎ𝐾₁+𝐻ℎ²𝐾₂𝐻+⋯ℎ⁹𝐾₉,𝐶²𝐻=𝑔ℎ1+ℎ𝐶₁+𝐻ℎ²𝐶₂𝐻+⋯+ℎ⁹𝐶₉, (73) where the formulae to determine 𝐾₁,…,𝐾₉ and 𝐶₁,…,𝐶₉ are all listed in [84].

The numerical results show that by applying the 9th order solitary wave solution [84] into (70) the wave generated can propagate without decay for a distance hundred times longer than the water depth, as shown in Figure 34.

Figure 34 

Evolution of the free surface wave for the case of 𝐻0=0.14 m and ℎ=0.4 m generated by the new solitary wave generation method. The subscript “0” denotes the original wave height setup in the wave generation control.

4.3.2. Liquid Sloshing in a Swaying Tank

A numerical model for simulating the free surface waves of the liquid sloshing in a swaying tank was developed in [53]. The model was verified by the comparison with the experimental data in [85]. The layout of the experiment is shown in Figure 35. This experiment was carried out to verify the numerical model developed in [85] which is a VOF (volume of fluid) model for the simulation of real fluid motion governed by the Reynolds-averaged Navier-Stokes equations. Among several cases carried out in [85], the nonbreaking case with the strongest nonlinear effect in the experiment is chosen as the verification of present model. The period of the oscillation is 1.0372 sec. The amplitude of the oscillation is 0.5 cm.

Figure 35 

Layout of the sloshing experiment of Liu and Lin (2008).

Discretizing a wave length with at least 20 segments, the initial nodal spacing on the free surface is chosen as 5.18 cm. The collocation points are initially distributed as a hexagonal close packing array so that the most compact nodal arrangement can be achieved. Therefore, the vertical nodal spacing on the side walls is 3 cm. Totally, there are 127 collocation points. The time step chosen in the simulation is 1/80 of the swaying period. The initial nodal distribution is shown in Figure 36.

Figure 36 

Initial nodal distribution of the numerical model of the liquid sloshing in a swaying tank.

Figure 37 shows the comparison of the numerical results with the experimental data. Very good agreement is found. It is also found that the higher the peak grows, the flatter the trough becomes. This indicates an increase of nonlinear effect as the tank oscillation continues. In the simulation the side walls are set to be infinitely high but in the experiment the tops of the two side walls are just 15 cm high from the still water level. What happened when the free surface elevation goes higher than the height of side walls had not been mentioned in [85]. In present simulation, free surface elevation goes over the top of the right side wall at 𝑡=6.184 sec. This might be the reason why the simulated free surface displacement, 𝜂, in the last wave period is slightly higher than observed. Figure 38 shows the positions of the traced fluid particles in the time interval of 𝑡=6.0𝑇∼6.5𝑇. This figure shows that, at the end of simulation, the run-up becomes much higher than the initial water depth. It is an indication of very strong nonlinearity on the free surface. Only 9 seconds is needed to simulate this case by using just one processing unit of Intel(R) Core(TM) i7-3370 CPU. The nodal spacing of present model is an order larger than the grid size in the model of [75]. Present model is much more efficient in this case: 𝑡=6.0𝑇∼6.5𝑇. (74)

(a)

(b)

(c)

Figure 37 

Comparison of the numerical results with the experimental data.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 38 

Snapshots of traced fluid particles in the time interval of 𝑡=6.0𝑇∼6.5𝑇.

5. Conclusions

Recent advancement on the modified finite point method is reviewed in this paper. A normalized parameter in the Gaussian distribution of the weighting factor is presented with recommendation of relevant parameters for better numerical results, regardless of the size of the computational domain. A concept of regional connectivity is introduced and demonstrated for computational domains consisting of different physical behaviors in subdomains. This concept will be needed for similar problems such as a crack in a thin plate and a cutoff wall in groundwater seepage.

Benchmark examples of employing this method to solve Laplace, Poisson, convection-diffusion, Helmholtz, mild-slope, and extended mild-slope equations are verified. This method could be employed to solve governing equations with real-valued unknowns as well as complex ones. Application of MFPM to numerical generation of orthogonal grids is illustrated to overcome difficulty of other methods. Problems of free-surface water waves represented in frequency- and time-domain are demonstrated for the applicability of MFPM. The derivatives of the potentials can be easily obtained as part of the solutions after the local collocations are achieved. This advantage of MFPM will provide great potentials in numerical simulations for physical problems when accurate derivatives of the solution functions are of particular interest.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors would like to thank the National Science Council of Taiwan for the financial support through Projects nos. 102-2221-E-002-155 and NSC 103-2911-I-006-302 in these years for developing this MFPM method.