# A Review on the Modified Finite Point Method

Authors: Nan-Jing Wu; Boe Shiun Chen; Ting-Kuei Tsay
Journal: Mathematical Problems in Engineering
 (2014)

---

## 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.


---

## Body

## 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 thehp-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)L{ϕ}=c1ϕ+c2∂ϕ∂x+c3∂ϕ∂y+c4∂2ϕ∂x2+c5∂2ϕ∂y2+c6∂2ϕ∂x∂y=s
subjected to specified boundary conditions
(2)B{ϕ}=q1ϕ+q2∂ϕ∂x+q3∂ϕ∂y=fiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiix⃑∈Γ1,(3)ϕ=ϕb,x⃑∈Γ2,
where c1,c2,…,c6, q1, q2, q3, f, and s are all functions of x, y and may change in time in transient problems. The boundary Γ1 could be nonsmooth and then at a corner there could be more than one condition like (2). Therefore, q1, q2, q3, and fcould 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 N nodes as needed. At each node, x→=x→j, ϕ is approximated as
(4)ϕ(x⃑)|x⃑≈x⃑j≈ϕ^j(x⃑)=∑i=1mαjipi(X⃑)
in which X⃑=x⃑-x⃑j is the relative position vector, pi(X⃑) is the ith monomial of the polynomial, and αji are coefficients to be determined. The subscript j indicates that this approximation is valid only in the vicinity of x⃑j. Once a new x⃑j is chosen, there will be a new set of αji. For a 2-D problem governed by a second order partial differential equation, the monomials are
(5){pi(X→),i=1~m}={1XYX2Y2XYX3X2YXY2Y3⋯}
in which X⃑=Xi⃑+Yj⃑. The value of m 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
(6)ϕ|x⃑=x⃑j≈ϕ^|x⃑=x⃑j=αj1,∂ϕ∂x|x⃑=x⃑j≈∂ϕ^∂x|x⃑=x⃑j=αj2,∂ϕ∂y|x⃑=x⃑j≈∂ϕ^∂y|x⃑=x⃑j=αj3,∂2ϕ∂x2|x⃑=x⃑j≈∂2ϕ^∂x2|x⃑=x⃑j=2αj4,∂2ϕ∂y2|x⃑=x⃑j≈∂2ϕ^∂y2|x⃑=x⃑j=2αj5,∂2ϕ∂x∂y|x⃑=x⃑j≈∂2ϕ^∂x∂y|x⃑=x⃑j=αj6,….
Here the error residual of the local approximation around x⃑=x⃑j is defined as
(7)Ej=∑l=1N(Wjl(ϕ(x⃑l)-ϕ^j(x⃑l))2),
where Wjl is a weighting factor determined by the distance between x⃑j and x⃑l. Usually, a normalized Gaussian function is selected for determining the weighting factor
(8)Wjl={exp⁡(-ε(rjl/ρj)2)-exp⁡(-ε)1-exp⁡(-ε),rjl<ρj0,rjl≥ρj,
where rjl is the distance between x⃑j and x⃑l (i.e., rjl=|x⃑l-x⃑j|), ε is the shape parameter, and ρj is the supporting range measured from the point of  x⃑j. The weight, Wjl, 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
(9)Ej=∑k=1n(Wjk(ϕ(x⃑k)-ϕ^j(x⃑k))2),
where k is the local index number of x⃑l in the jth subdomain and n 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αji corresponding to the minimized Ej in (9) can be expressed explicitly by the least square approach as
(10)[αj]=[Λ][β]
in which
(11)[αj]m×1=[αj1⋯αji⋯αjm]T[Λ]m×n=([A]T[A])-1[A]T[A]n×m=[a11a12⋯⋯a1ma21⋱a2m⋮aki⋮⋮⋱⋮an1an2⋯⋯anm][β]n×1=[β1⋯βk⋯βn]T=[w1ϕ1⋯wkϕk⋯wnϕn]T,
where wk=Wjk,  ϕk=ϕ(x⃑k), and aki=wkpi(x⃑k-x⃑j).It should be noted that this approximation is only valid in the vicinity ofx⃑j. Once a new x⃑j 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
(12)[K]N×N[ϕ]N×1=[b]N×1.
In the case that x⃑j is on the boundary and the value of ϕ is known, the entities in (12) are
(13)κjl={1,ifj=l0,otherwise(14)bj=ϕj.
When the collocation point is inside the domain, the entities in (12) are
(15)κjl={wk(c1λ1k+c2λ2k+c3λ3k+2c4λ4k+2c5λ5k+c6λ6k)x⃑=x⃑j,if|x⃑j-x⃑l|<ρj0,otherwisebj=s|x⃑=x⃑j,
where λik indicates the ith entry in kth 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
(16)κjl={wk(q1λ1k+q2λ2k+q3λ3k)x⃑=x⃑j,if|x⃑j-x⃑l|<ρj0,otherwisebj=f|x⃑=x⃑j.


### 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 x⃑j (9) under the condition that
(17)(L{ϕ}-s)2+(B1{ϕ}-f1)2+⋯+(Bnnd{ϕ}-fnnd)2⟶0,
where nnd is the number of non-Dirichlet boundary conditions at the node x⃑=x⃑j, can be expressed as
(18)[αj]=[Λ′][ββ′]
in which
(19)[Λ′]m×(n+nnd+1)=([AA′]T[AA′])-1[AA′]T[A′](nnd+1)×m=[w′c1⋯⋯⋯w′cp0⋯0w′q1,1⋯w′q1,30⋯⋯⋯0⋮⋮⋮⋮w′qnnd,1⋯w′qnnd,30⋯⋯⋯0][β′](nnd+1)×1=[w′sw′f1⋯w′fnnd]T,
where the definitions of wk, ϕk, aki are all the same as in the conventional collocation method, while w′=W′ in which W′ represents a penalty weighting factor whose value is much greater than 1. In the case of nnd 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., nnd=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
(20)κjl={wkλ1k′-δj(l),if|x⃑j-x⃑l|<ρj0,otherwiseδj(l)={1,ifj=l0,otherwisebj=-∑k=1nnd+1λ1,(n+k)′βk′,whereλ1k′ 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, ρj, and the degree of polynomial for the local approximation. The penalty weighting factor W′ 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, ρj as 1.05 times of the distance from x⃑j to its 25th nearest neighboring node, and W′=104 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.


## 2.1. The Local Approximation

When solving a general 2-D linear, second order PDE is(1)L{ϕ}=c1ϕ+c2∂ϕ∂x+c3∂ϕ∂y+c4∂2ϕ∂x2+c5∂2ϕ∂y2+c6∂2ϕ∂x∂y=s
subjected to specified boundary conditions
(2)B{ϕ}=q1ϕ+q2∂ϕ∂x+q3∂ϕ∂y=fiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiix⃑∈Γ1,(3)ϕ=ϕb,x⃑∈Γ2,
where c1,c2,…,c6, q1, q2, q3, f, and s are all functions of x, y and may change in time in transient problems. The boundary Γ1 could be nonsmooth and then at a corner there could be more than one condition like (2). Therefore, q1, q2, q3, and fcould 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 N nodes as needed. At each node, x→=x→j, ϕ is approximated as
(4)ϕ(x⃑)|x⃑≈x⃑j≈ϕ^j(x⃑)=∑i=1mαjipi(X⃑)
in which X⃑=x⃑-x⃑j is the relative position vector, pi(X⃑) is the ith monomial of the polynomial, and αji are coefficients to be determined. The subscript j indicates that this approximation is valid only in the vicinity of x⃑j. Once a new x⃑j is chosen, there will be a new set of αji. For a 2-D problem governed by a second order partial differential equation, the monomials are
(5){pi(X→),i=1~m}={1XYX2Y2XYX3X2YXY2Y3⋯}
in which X⃑=Xi⃑+Yj⃑. The value of m 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
(6)ϕ|x⃑=x⃑j≈ϕ^|x⃑=x⃑j=αj1,∂ϕ∂x|x⃑=x⃑j≈∂ϕ^∂x|x⃑=x⃑j=αj2,∂ϕ∂y|x⃑=x⃑j≈∂ϕ^∂y|x⃑=x⃑j=αj3,∂2ϕ∂x2|x⃑=x⃑j≈∂2ϕ^∂x2|x⃑=x⃑j=2αj4,∂2ϕ∂y2|x⃑=x⃑j≈∂2ϕ^∂y2|x⃑=x⃑j=2αj5,∂2ϕ∂x∂y|x⃑=x⃑j≈∂2ϕ^∂x∂y|x⃑=x⃑j=αj6,….
Here the error residual of the local approximation around x⃑=x⃑j is defined as
(7)Ej=∑l=1N(Wjl(ϕ(x⃑l)-ϕ^j(x⃑l))2),
where Wjl is a weighting factor determined by the distance between x⃑j and x⃑l. Usually, a normalized Gaussian function is selected for determining the weighting factor
(8)Wjl={exp⁡(-ε(rjl/ρj)2)-exp⁡(-ε)1-exp⁡(-ε),rjl<ρj0,rjl≥ρj,
where rjl is the distance between x⃑j and x⃑l (i.e., rjl=|x⃑l-x⃑j|), ε is the shape parameter, and ρj is the supporting range measured from the point of  x⃑j. The weight, Wjl, 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
(9)Ej=∑k=1n(Wjk(ϕ(x⃑k)-ϕ^j(x⃑k))2),
where k is the local index number of x⃑l in the jth subdomain and n 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αji corresponding to the minimized Ej in (9) can be expressed explicitly by the least square approach as
(10)[αj]=[Λ][β]
in which
(11)[αj]m×1=[αj1⋯αji⋯αjm]T[Λ]m×n=([A]T[A])-1[A]T[A]n×m=[a11a12⋯⋯a1ma21⋱a2m⋮aki⋮⋮⋱⋮an1an2⋯⋯anm][β]n×1=[β1⋯βk⋯βn]T=[w1ϕ1⋯wkϕk⋯wnϕn]T,
where wk=Wjk,  ϕk=ϕ(x⃑k), and aki=wkpi(x⃑k-x⃑j).It should be noted that this approximation is only valid in the vicinity ofx⃑j. Once a new x⃑j 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
(12)[K]N×N[ϕ]N×1=[b]N×1.
In the case that x⃑j is on the boundary and the value of ϕ is known, the entities in (12) are
(13)κjl={1,ifj=l0,otherwise(14)bj=ϕj.
When the collocation point is inside the domain, the entities in (12) are
(15)κjl={wk(c1λ1k+c2λ2k+c3λ3k+2c4λ4k+2c5λ5k+c6λ6k)x⃑=x⃑j,if|x⃑j-x⃑l|<ρj0,otherwisebj=s|x⃑=x⃑j,
where λik indicates the ith entry in kth 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
(16)κjl={wk(q1λ1k+q2λ2k+q3λ3k)x⃑=x⃑j,if|x⃑j-x⃑l|<ρj0,otherwisebj=f|x⃑=x⃑j.


## 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 x⃑j (9) under the condition that
(17)(L{ϕ}-s)2+(B1{ϕ}-f1)2+⋯+(Bnnd{ϕ}-fnnd)2⟶0,
where nnd is the number of non-Dirichlet boundary conditions at the node x⃑=x⃑j, can be expressed as
(18)[αj]=[Λ′][ββ′]
in which
(19)[Λ′]m×(n+nnd+1)=([AA′]T[AA′])-1[AA′]T[A′](nnd+1)×m=[w′c1⋯⋯⋯w′cp0⋯0w′q1,1⋯w′q1,30⋯⋯⋯0⋮⋮⋮⋮w′qnnd,1⋯w′qnnd,30⋯⋯⋯0][β′](nnd+1)×1=[w′sw′f1⋯w′fnnd]T,
where the definitions of wk, ϕk, aki are all the same as in the conventional collocation method, while w′=W′ in which W′ represents a penalty weighting factor whose value is much greater than 1. In the case of nnd 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., nnd=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
(20)κjl={wkλ1k′-δj(l),if|x⃑j-x⃑l|<ρj0,otherwiseδj(l)={1,ifj=l0,otherwisebj=-∑k=1nnd+1λ1,(n+k)′βk′,whereλ1k′ 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, ρj, and the degree of polynomial for the local approximation. The penalty weighting factor W′ 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, ρj as 1.05 times of the distance from x⃑j to its 25th nearest neighboring node, and W′=104 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 isa. Viscous effects which could be significant in real fluids are not considered here. The velocity at x=±∞ or y=±∞ is v⃑=Ui⃑+0j⃑. The velocity potential ϕ is governed by the Laplace equation:
(21)∂2ϕ∂x2+∂2ϕ∂y2=0.
The velocity field is determined as v⃑=ui⃑+vj⃑:
(22)u=∂ϕ∂xv=∂ϕ∂y.
There is an exact solution expressed as
(23)ϕe=Ux(1+a2x2+y2).
The subscript “e” denotes the “exact solution.”Due to symmetry, only the upper left quadrant is modeled. The values ofU and a are chosen as 1.0 and 0.5, respectively. These quantities are dimensionless. The computational domain is within the range of -2≤x≤0 and 0≤y≤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.Two types of nodal arrangement for a potential uniform flow passing a circular cylinder.
(a)
(b)The profile of computed flow speed (V=u2+v2) 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.Computed flow speed profile along the body surface in the problem of a uniform flow passing a circular cylinder by two different nodal arrangements.
(a)
(b)


### 3.2. Poisson Equation

Consider the following 2-D Poison equation:(24)∂2ϕ∂x2+∂2ϕ∂y2=sin(πx)sin(πy)
within the domain of 0≤x≤1, 0≤y≤1. The boundary conditions are
(25)∂ϕ∂y=-12πsin(πx)cos⁡(πy),alongy=0,y=1ϕ=0,alongx=0,x=1.The exact solution of this Poisson equation can be found asϕe=-sin(πx)sin(πy)/(2π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.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 inx- and y-direction, respectively.
(a)
(b)


### 3.3. Convection-Diffusion Equation

Consider the following steady-state convection-diffusion equation [54]:
(26)iiiiiiiiiiiiiiiiiiiiiiiiii∂2ϕ∂x2+∂2ϕ∂y2-U∂ϕ∂x=0intheregionof0≤x≤1,0≤y≤0.6,
with boundary conditions:
(27)ϕ|x=0=1,ϕ|x=1=2,∂ϕ∂y|y=0=0,∂ϕ∂y|y=0.6=0.
This problem has an exact solution as
(28)ϕe=2-1-exp⁡[U(x-1)]1-exp⁡(-U).A strong convective case ofU=100 is chosen. The nodal arrangement for this problem is shown in Figure 6. The nodal spacing is mainly Δx=Δy=0.02 while near the right boundary the nodal density is increased to Δx=Δy=0.01, in order to catch the behavior of exponential growth near the boundary of x=1. Totally 1,825 nodes are used in this test. Only 1.48% of the entries in the global matrix [K] are nonzero. Profiles of ϕ and ∂ϕ/∂x along y=0 and y=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.Results of problem of the steady-state convection-diffusion equation.
(a)
(b)


### 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(29)∂2ϕ∂x2+∂2ϕ∂y2+k2ϕ=0,
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 k is defined as 2π/L, in which L is the wavelength of the monochromatic waves. The real free-surface velocity potential of monochromatic waves is defined as
(30)Φ=Re{ϕe-iωt}
in which i=-1 is the pure imaginary number. The symbol ω indicating the angular frequency is defined as ω=2π/T and T 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 Figure8, 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, ϕS:
(31)lim⁡r→∞(∂ϕS∂r-ikϕS)=0.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,ϕ=ϕI+ϕS, and (31) can be implemented at an artificial open boundary far away from the scatter as
(32)∂ϕ∂r-ikϕ=∂ϕI∂r-ikϕI,atr=r0≫0,
where r=x2+y2 and ϕI=aIeiS is the pseudo free surface velocity potential of incident wave in which S denotes its phase function, and aI denotes wave amplitude. For waves propagating toward the +x direction, the phase S of the incoming wave is thus ikx. The no-flux boundary condition given at the cylinder surface is
(33)∂ϕ∂r=0,atr=a,
where a 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
(34)ϕ=aI∑n=0∞εnin{Jn(kr)Hn(kr)×Jn′(ka)Hn′(ka)}cos⁡n(θ-θI),
where
(35)ε={1,ifn=02,ifn>0
and θI is the incident angle of waves, and Jn and Hn are the Bessel and Hankel functions of the first kind of order n. 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 r and θ, respectively:
(36)Ur=∂ϕ∂r=gaIω×∑n=0∞εnin-1{kJn′(kr)-kHn′(kr)kJn′(ka)kH′(ka)}×cos⁡n(θ-θI)Uθ=1r∂ϕ∂θ=gaIωr×∑n=0∞εnin-1{-Jn(kr)+kHn′(kr)kJn′(ka)Hn′(ka)}×nsinn(θ-θI).A case of water depth of 10 m and incident waves with periodT=20 s, which corresponds to wave number k=0.0308 m−1, is chosen for demonstration. The radius of the cylinder is chosen as 30 m while the open outer boundary is located at r0=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.Contours of equal-amplitude of waves scattered by a cylinder in constant water depth,T=20 s. (a) MFPM results; (b) analytical solutions.
(a)
(b)Contours of equal-phase of waves scattered by a cylinder on constant water depth,T=20 s. (a) MFPM results; (b) analytical solutions.
(a)
(b)Contours of magnitude of velocity of waves scattered by a cylinder on constant water depth,T=20 s. (a) MFPM results; (b) analytical solutions.
(a)
(b)Comparisons of MFPM results and analytical solutions of amplitude, phase, and magnitude of velocity around atr=a.
(a)
(b)
(c)


## 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 isa. Viscous effects which could be significant in real fluids are not considered here. The velocity at x=±∞ or y=±∞ is v⃑=Ui⃑+0j⃑. The velocity potential ϕ is governed by the Laplace equation:
(21)∂2ϕ∂x2+∂2ϕ∂y2=0.
The velocity field is determined as v⃑=ui⃑+vj⃑:
(22)u=∂ϕ∂xv=∂ϕ∂y.
There is an exact solution expressed as
(23)ϕe=Ux(1+a2x2+y2).
The subscript “e” denotes the “exact solution.”Due to symmetry, only the upper left quadrant is modeled. The values ofU and a are chosen as 1.0 and 0.5, respectively. These quantities are dimensionless. The computational domain is within the range of -2≤x≤0 and 0≤y≤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.Two types of nodal arrangement for a potential uniform flow passing a circular cylinder.
(a)
(b)The profile of computed flow speed (V=u2+v2) 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.Computed flow speed profile along the body surface in the problem of a uniform flow passing a circular cylinder by two different nodal arrangements.
(a)
(b)


## 3.2. Poisson Equation

Consider the following 2-D Poison equation:(24)∂2ϕ∂x2+∂2ϕ∂y2=sin(πx)sin(πy)
within the domain of 0≤x≤1, 0≤y≤1. The boundary conditions are
(25)∂ϕ∂y=-12πsin(πx)cos⁡(πy),alongy=0,y=1ϕ=0,alongx=0,x=1.The exact solution of this Poisson equation can be found asϕe=-sin(πx)sin(πy)/(2π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.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 inx- and y-direction, respectively.
(a)
(b)


## 3.3. Convection-Diffusion Equation

Consider the following steady-state convection-diffusion equation [54]:
(26)iiiiiiiiiiiiiiiiiiiiiiiiii∂2ϕ∂x2+∂2ϕ∂y2-U∂ϕ∂x=0intheregionof0≤x≤1,0≤y≤0.6,
with boundary conditions:
(27)ϕ|x=0=1,ϕ|x=1=2,∂ϕ∂y|y=0=0,∂ϕ∂y|y=0.6=0.
This problem has an exact solution as
(28)ϕe=2-1-exp⁡[U(x-1)]1-exp⁡(-U).A strong convective case ofU=100 is chosen. The nodal arrangement for this problem is shown in Figure 6. The nodal spacing is mainly Δx=Δy=0.02 while near the right boundary the nodal density is increased to Δx=Δy=0.01, in order to catch the behavior of exponential growth near the boundary of x=1. Totally 1,825 nodes are used in this test. Only 1.48% of the entries in the global matrix [K] are nonzero. Profiles of ϕ and ∂ϕ/∂x along y=0 and y=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.Results of problem of the steady-state convection-diffusion equation.
(a)
(b)


## 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(29)∂2ϕ∂x2+∂2ϕ∂y2+k2ϕ=0,
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 k is defined as 2π/L, in which L is the wavelength of the monochromatic waves. The real free-surface velocity potential of monochromatic waves is defined as
(30)Φ=Re{ϕe-iωt}
in which i=-1 is the pure imaginary number. The symbol ω indicating the angular frequency is defined as ω=2π/T and T 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 Figure8, 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, ϕS:
(31)lim⁡r→∞(∂ϕS∂r-ikϕS)=0.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,ϕ=ϕI+ϕS, and (31) can be implemented at an artificial open boundary far away from the scatter as
(32)∂ϕ∂r-ikϕ=∂ϕI∂r-ikϕI,atr=r0≫0,
where r=x2+y2 and ϕI=aIeiS is the pseudo free surface velocity potential of incident wave in which S denotes its phase function, and aI denotes wave amplitude. For waves propagating toward the +x direction, the phase S of the incoming wave is thus ikx. The no-flux boundary condition given at the cylinder surface is
(33)∂ϕ∂r=0,atr=a,
where a 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
(34)ϕ=aI∑n=0∞εnin{Jn(kr)Hn(kr)×Jn′(ka)Hn′(ka)}cos⁡n(θ-θI),
where
(35)ε={1,ifn=02,ifn>0
and θI is the incident angle of waves, and Jn and Hn are the Bessel and Hankel functions of the first kind of order n. 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 r and θ, respectively:
(36)Ur=∂ϕ∂r=gaIω×∑n=0∞εnin-1{kJn′(kr)-kHn′(kr)kJn′(ka)kH′(ka)}×cos⁡n(θ-θI)Uθ=1r∂ϕ∂θ=gaIωr×∑n=0∞εnin-1{-Jn(kr)+kHn′(kr)kJn′(ka)Hn′(ka)}×nsinn(θ-θI).A case of water depth of 10 m and incident waves with periodT=20 s, which corresponds to wave number k=0.0308 m−1, is chosen for demonstration. The radius of the cylinder is chosen as 30 m while the open outer boundary is located at r0=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.Contours of equal-amplitude of waves scattered by a cylinder in constant water depth,T=20 s. (a) MFPM results; (b) analytical solutions.
(a)
(b)Contours of equal-phase of waves scattered by a cylinder on constant water depth,T=20 s. (a) MFPM results; (b) analytical solutions.
(a)
(b)Contours of magnitude of velocity of waves scattered by a cylinder on constant water depth,T=20 s. (a) MFPM results; (b) analytical solutions.
(a)
(b)Comparisons of MFPM results and analytical solutions of amplitude, phase, and magnitude of velocity around atr=a.
(a)
(b)
(c)


## 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 (x-y 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 (x-y coordinate) to the computational domain (ξ-η coordinate) can be reduced to a problem governed by two Laplace equations:(37a)∂2ξ∂x2+∂2ξ∂y2=0;(37b)∂2η∂x2+∂2η∂y2=0with the boundary conditions(38a)n⃑·∇η=0on⁡ξ=0andξ=ξ0(38b)n⃑·∇ξ=0on⁡η=0andη=η0,wheren⃑ is the outward normal vector of the boundaries. After (37a) and (37b) are solved separately, the obtained x-y 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:(39a)∂2x∂ξ2+∂2x∂η2=0;(39b)∂2y∂ξ2+∂2y∂η2=0.The relations between the partial derivatives are listed below:(40a)∂ξ∂x=1J∂y∂η;(40b)∂ξ∂y=-1J∂x∂η(40c)∂η∂x=-1J∂y∂ξ;(40d)∂η∂y=1J∂x∂ξ,whereJ is the Jacobian and is defined as
(41)J=∂x∂ξ∂y∂η-∂x∂η∂y∂ξ.A sketch of the transformation is illustrated in Figure14. A detailed derivation of these equations could be found in [58].Sketch of a transformation between two orthogonal coordinates from the physical domain to the transformed domain.
(a)
(b)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 ξ0 and η0 are all set to be 1 because of the symmetry. For forward transformation from x-y 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:(42a)AREO=∑j=1N(∇ξ·∇η)x⃑=x⃑j2∑j=1N(|∇ξ|2|∇η|2)x⃑=x⃑j(42b)MREO=max⁡{(∇ξ·∇η)x⃑=x⃑j2,j=1~N}(1N∑j=1N(|∇ξ|2|∇η|2)x⃑=x⃑j).The subscript “O” 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-5 and 7.53×10-4, 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 ∂ξ/∂x, ∂ξ/∂y, ∂η/∂x, and ∂η/∂y 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
(43)∇·(CCg∇η)+k2CCgη=0,
where C=ω/k and Cg=(1/2)(1+2kh/sinh⁡2kh)C represent wave celerity and group celerity, respectively; k 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
(44)∂η∂n=0.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 Figure19, the depth between the shoreline and outer boundary of the shoal is described as
(45)h=h0(rr0)2,a≤r≤r0,
where r is the distance measured from the center of the island, a is radius of the circular island, r0 is the radius of paraboloidal shoal, and h0 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 fora≤r≤r0, derived by Homma [71], are presented as the following for completeness:
(46)η=a0∑n=0∞kεni2n[Jn′(kr0)Hn(kr0)-Jn(kr0)Hn′(lr0)]×(αn+1αn-1a-2αnrαn-1+r-αn-1)cos⁡nθ+[(αn+1αn-1)(αn+1)(a-2αnr0αn-2-r0-αn-2)Hn(kr0)-kHn′(kr0)(αn+1αn-1)a-2αnr0αn-1+r0-αn-1].θI is the incident angle of propagating waves and Jn and Hn are the Bessel and Hankel functions of the first kind of order n, as in the fourth benchmark problem. The superscript of the symbol of prime indicates their derivatives, respectively.Differentiating (47) with respect to radial (r) and angular (θ) variables, respectively, the analytical solutions of radial and tangential velocities in the inner region of a≤r≤r0 can be obtained analytically with any difficulty.A case for a circular island with radius,a=10 km, as shown in Figure 19, is considered. The water depth is constant, h0=4 km, at distance greater than r0=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 t=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 x- and y-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.Contours of wave amplitude lines around a circular island on a paraboloidal shoal,T=720 s. (a) Numerical result; (b) analytic results.
(a)
(b)Contours of wave phase lines around a circular island on a paraboloidal shoal,T=720 s, at t=0. (a) Numerical results; (b) analytic results.
(a)
(b)Contours of magnitude of velocity around a circular island on a paraboloidal shoal,T=720 s, t=0. (a) Numerical results; (b) analytic results.
(a)
(b)Numerical and analytic results of amplitude, phase, and magnitude of velocity around atr=a.
(a)
(b)
(c)


#### 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 Figure24 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 Figure25. 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:
(47)ηH=a0ND∑n=0∞εn2πJ0(nδ)kaJn′(ka)Jn(kr)cos⁡n(θ-θI),
where ηH is water free-surface displacement within the harbor, a0 is incident wave amplitude, δ is half of the harbor entrance opening angle, θ0 is the angle from the positive x-coordinate to the centerline of entrance opening, θI is direction of incident waves with respect to positive x-axis, a is radius of the circular harbor, k is incident wave number, and N is weighted average of the forcing pressure:
(48)N=∑n=0∞2iπka2εnJ0(nδ)Hn′(ka)×(cos⁡nπ2cos⁡nθIcos⁡θ0+isinnπ2sinnθIsinnθ0),Figure 25
Definition of a circular harbor with protruding breakwaters.andD is term of the coefficient:
(49)D=∑n=0∞εn2πJn(ka)J02(nδ)kaJn′(ka)-∑n=0∞εn2πHn(ka)J02(nδ)kaHn′(ka)(1+cos⁡nπcos⁡2nθ0).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θI of 225° is presented. Present numerical results for an idealized circular harbor of radius a=0.75 ft with harbor mouth 60° opening (2δ).Response curves atr=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 ka=0.5 and the amplification factor, R, 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 ka=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 atr=0.7 ft, θ=0° of a circular harbor of radius, a=0.75 ft, 60° opening for θI=225°.Figure 27
Contours of wave amplitude of a circular harbor of radius,a=0.75 ft, 60° opening, with protruding breakwaters at ka=0.5 for θI=225°.Velocity vectors of a circular harbor of radius,a=0.75 ft, 60° opening, with protruding breakwaters at ka=0.5. (a) t=0, (b) t=T/2 for θI=225°.
(a)
(b)


#### 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]:
(50)∇·(PG∇η)+λ2PG(1-f2ω2)η=0,
where f 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, ω2=gλtanhλh, which may be real or imaginary values; P and G represent, respectively, the wave and group celerity as
(51)P=ωλG=12(1+2λhsinh⁡2λh)P,
and k is the wave number which could be expressed as
(52)k2=λ2(1-f2ω2).For real value ofλ, k is a real value when f<ω, 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)∂u∂t-fv=-g∂ς∂x∂v∂t+fu=-g∂ς∂y∂w∂t=-g∂ς∂z.In this section, Coriolis effects on long waves propagating in an infinite channel are investigated. The channel has two parallel side vertical walls of widthL, which is shown in Figure 29. Water depth is H0 and waves propagating from the left to the right.Figure 29
An infinite channel of widthL, with Coriolis parameter f.According to Pedlosky [75], the analytic solution in Kelvin wave could be expressed as
(54)η=Re[ηoe-fy/Coeiλx}u=η0H0C0e-fy/C0cos⁡(kx±ωt+φ),
where C0=(gH0)1/2; H0 is constant depth of the channel; g is the gravity; η0 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 withF(z)=cosh⁡λ(z+h)/cosh⁡λh, the fluid displacement, ς, can be expressed as
(55)ς(x,y,z;t)=η(x,y)F(z)e-iωtη(x,y)=A(x,y)ei∫kn·dx⃑=A(x,y)ei∫(kxdx+kydy),
where A is amplitude of free surface displacement, η; k is wave number; kx and ky are the components of k in x- and y-direction, respectively.Assuming that the flow velocities associated with wave action,(u,v), also have the same vertical distribution and with (U,V) as the components at the free surface:
(56)(u,v)=(U,V)F(z).From (53) and (55)-(56), the velocity components formula with Coriolis effects can be derived and expressed as
(57)U=gfω2-f2∂η∂y-igωω2-f2∂η∂xV=-gfω2-f2∂η∂x-igωω2-f2∂η∂y.Then from the no-flux condition,U→·n→=0, the solid boundary condition could be expressed as
(58)∂η∂n=±ifω∂η∂s,
where n and s 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
(59)∂ηS∂n-iληS=0.In order to illustrate the Coriolis effects of waves propagating in a channel with widthL=60 km and water depth H0=3 m have been used. A wave period T=46800 sec, amplitude = 1.0 m, and Coriolis parameter f=0.0000935638 sec−1 (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 c/f, the Rossby radius of deformation. The distribution of the wave velocity is parallel to the walls of the channel and also proportional to c/f.The results of analytical solutions: (a) distribution of wave amplitude; (b) distribution of wave velocity.
(a)
(b)For numerical computations, the incident wave amplitude is chosen to beηoe-fy/Co 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 x-direction between the numerical and analytical results at cross section x = 10,000 m are further compared as shown in Figure 32. Very good agreement between numerical and analytical results is observed. The velocity of y-direction is equal to 0 with a negligible error. Very good agreement between numerical results and analytical solutions is observed.The results of numerical model: (a) distribution of wave amplitude; (b) distribution of wave velocity.
(a)
(b)Comparisons of amplitude and velocity atx-diretion between numerical and analytical results at the cross section x = 10,000 m.
(a)
(b)
(c)


### 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
(60)∇2ϕ=0,
where ϕ(x⃑,t) is the velocity potential with x⃑ being the position vector and t being the time. The velocity vector is defined as v⃑=∇ϕ ((∂ϕ/∂x)i⃑+(∂ϕ/∂y)j⃑ for 2-D problems and (∂ϕ/∂x)i⃑+(∂ϕ/∂y)j⃑+(∂ϕ/∂z)k⃑ for 3-D problems). The kinematic free-surface boundary condition (KFSBC) and dynamic free-surface boundary condition (DFSBC) can be expressed as
(61)∇ϕ=dx⃑dt,x⃑∈theentiredomain,iiiiiiiiiiiiiincludingthefreesurface,dϕdt=-gz+12∇ϕ·∇ϕ,x⃑∈thefreesurface,
where g is the gravitational acceleration, z is the vertical coordinate, and C(t) 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)n⃑·∇ϕ=n⃑·v⃑b,
where n⃑ is the unit normal vector outward from the body surface, and v⃑b 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
(63)(ϕ|x⃑=x⃑j)(n)=(ϕ|x⃑=x⃑j)(n-2)+2Δt[(-gz+12∇ϕ·∇ϕ)x⃑=x⃑j](n-1),
where x⃑j denotes the position of the jth 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 nth 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,” x⃑j(n). 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 x⃑j(n):
(64)x⃑j(n)=x⃑j(n-2)+2Δt(∇ϕ|x⃑=x⃑j)(n-1).In the case of strong nonlinearity, the free-surface particle trajectoryx⃑j(n) should be modified by applying the Crank-Nicolson formula for a better numerical stability:
(65)x⃑j(n)=x⃑j(n-1)+Δt2[(∇ϕ|x⃑=x⃑j)(n)+(∇ϕ|x⃑=x⃑j)(n-1)].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 atx⃑j(n) 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)dξdt=u-(ξ,t)=Cη|x=ξh+η|x=ξ,
where ξ(t) is the position of the wave paddle at time t, while t is elapsed time since the start of the motion, C is the wave celerity, η is the free surface displacement, and h is the still water depth. In [82] it was proposed to apply Boussinesq’s solitary wave solution [83] into (66). That is,
(67)η=HS2,
where H is the wave height, and
(68)S=cosh⁡-1[KX]
in which K is the boundary outskirt decay coefficient, and
(69)X=(ξ-Ct-x0),
where x0 is the initial position of the wave crest. Because the solitary wave is generated in a quiescent water flume, x0 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
(70)K=3H4h3,
where the wave celerity is determined as
(71)C=g(h+H).Figure33 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].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.
(a)
(b)
(c)
(d)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)η=h((Hh)η1+(Hh)2η2+⋯+(Hh)9η9).The formulae to determine the values ofη1 to η9 could be found in [84]. In an analogous way, the boundary outskirt decay coefficient and the wave celerity are determined as
(73)K=3H4h3(1+(Hh)K1+(Hh)2K2+⋯(Hh)9K9),C2=gh(1+(Hh)C1+(Hh)2C2+⋯+(Hh)9C9),
where the formulae to determine K1,…,K9 and C1,…,C9 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 ofH0=0.14 m and h=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 Figure36.Figure 36
Initial nodal distribution of the numerical model of the liquid sloshing in a swaying tank.Figure37 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 t=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 t=6.0T~6.5T. 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:
(74)t=6.0T~6.5T.Comparison of the numerical results with the experimental data.
(a)
(b)
(c)Snapshots of traced fluid particles in the time interval oft=6.0T~6.5T.
(a)
(b)
(c)
(d)
(e)
(f)


## 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 (x-y 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 (x-y coordinate) to the computational domain (ξ-η coordinate) can be reduced to a problem governed by two Laplace equations:(37a)∂2ξ∂x2+∂2ξ∂y2=0;(37b)∂2η∂x2+∂2η∂y2=0with the boundary conditions(38a)n⃑·∇η=0on⁡ξ=0andξ=ξ0(38b)n⃑·∇ξ=0on⁡η=0andη=η0,wheren⃑ is the outward normal vector of the boundaries. After (37a) and (37b) are solved separately, the obtained x-y 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:(39a)∂2x∂ξ2+∂2x∂η2=0;(39b)∂2y∂ξ2+∂2y∂η2=0.The relations between the partial derivatives are listed below:(40a)∂ξ∂x=1J∂y∂η;(40b)∂ξ∂y=-1J∂x∂η(40c)∂η∂x=-1J∂y∂ξ;(40d)∂η∂y=1J∂x∂ξ,whereJ is the Jacobian and is defined as
(41)J=∂x∂ξ∂y∂η-∂x∂η∂y∂ξ.A sketch of the transformation is illustrated in Figure14. A detailed derivation of these equations could be found in [58].Sketch of a transformation between two orthogonal coordinates from the physical domain to the transformed domain.
(a)
(b)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 ξ0 and η0 are all set to be 1 because of the symmetry. For forward transformation from x-y 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:(42a)AREO=∑j=1N(∇ξ·∇η)x⃑=x⃑j2∑j=1N(|∇ξ|2|∇η|2)x⃑=x⃑j(42b)MREO=max⁡{(∇ξ·∇η)x⃑=x⃑j2,j=1~N}(1N∑j=1N(|∇ξ|2|∇η|2)x⃑=x⃑j).The subscript “O” 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-5 and 7.53×10-4, 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 ∂ξ/∂x, ∂ξ/∂y, ∂η/∂x, and ∂η/∂y 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
(43)∇·(CCg∇η)+k2CCgη=0,
where C=ω/k and Cg=(1/2)(1+2kh/sinh⁡2kh)C represent wave celerity and group celerity, respectively; k 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
(44)∂η∂n=0.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 Figure19, the depth between the shoreline and outer boundary of the shoal is described as
(45)h=h0(rr0)2,a≤r≤r0,
where r is the distance measured from the center of the island, a is radius of the circular island, r0 is the radius of paraboloidal shoal, and h0 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 fora≤r≤r0, derived by Homma [71], are presented as the following for completeness:
(46)η=a0∑n=0∞kεni2n[Jn′(kr0)Hn(kr0)-Jn(kr0)Hn′(lr0)]×(αn+1αn-1a-2αnrαn-1+r-αn-1)cos⁡nθ+[(αn+1αn-1)(αn+1)(a-2αnr0αn-2-r0-αn-2)Hn(kr0)-kHn′(kr0)(αn+1αn-1)a-2αnr0αn-1+r0-αn-1].θI is the incident angle of propagating waves and Jn and Hn are the Bessel and Hankel functions of the first kind of order n, as in the fourth benchmark problem. The superscript of the symbol of prime indicates their derivatives, respectively.Differentiating (47) with respect to radial (r) and angular (θ) variables, respectively, the analytical solutions of radial and tangential velocities in the inner region of a≤r≤r0 can be obtained analytically with any difficulty.A case for a circular island with radius,a=10 km, as shown in Figure 19, is considered. The water depth is constant, h0=4 km, at distance greater than r0=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 t=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 x- and y-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.Contours of wave amplitude lines around a circular island on a paraboloidal shoal,T=720 s. (a) Numerical result; (b) analytic results.
(a)
(b)Contours of wave phase lines around a circular island on a paraboloidal shoal,T=720 s, at t=0. (a) Numerical results; (b) analytic results.
(a)
(b)Contours of magnitude of velocity around a circular island on a paraboloidal shoal,T=720 s, t=0. (a) Numerical results; (b) analytic results.
(a)
(b)Numerical and analytic results of amplitude, phase, and magnitude of velocity around atr=a.
(a)
(b)
(c)


### 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 Figure24 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 Figure25. 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:
(47)ηH=a0ND∑n=0∞εn2πJ0(nδ)kaJn′(ka)Jn(kr)cos⁡n(θ-θI),
where ηH is water free-surface displacement within the harbor, a0 is incident wave amplitude, δ is half of the harbor entrance opening angle, θ0 is the angle from the positive x-coordinate to the centerline of entrance opening, θI is direction of incident waves with respect to positive x-axis, a is radius of the circular harbor, k is incident wave number, and N is weighted average of the forcing pressure:
(48)N=∑n=0∞2iπka2εnJ0(nδ)Hn′(ka)×(cos⁡nπ2cos⁡nθIcos⁡θ0+isinnπ2sinnθIsinnθ0),Figure 25
Definition of a circular harbor with protruding breakwaters.andD is term of the coefficient:
(49)D=∑n=0∞εn2πJn(ka)J02(nδ)kaJn′(ka)-∑n=0∞εn2πHn(ka)J02(nδ)kaHn′(ka)(1+cos⁡nπcos⁡2nθ0).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θI of 225° is presented. Present numerical results for an idealized circular harbor of radius a=0.75 ft with harbor mouth 60° opening (2δ).Response curves atr=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 ka=0.5 and the amplification factor, R, 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 ka=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 atr=0.7 ft, θ=0° of a circular harbor of radius, a=0.75 ft, 60° opening for θI=225°.Figure 27
Contours of wave amplitude of a circular harbor of radius,a=0.75 ft, 60° opening, with protruding breakwaters at ka=0.5 for θI=225°.Velocity vectors of a circular harbor of radius,a=0.75 ft, 60° opening, with protruding breakwaters at ka=0.5. (a) t=0, (b) t=T/2 for θI=225°.
(a)
(b)


### 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]:
(50)∇·(PG∇η)+λ2PG(1-f2ω2)η=0,
where f 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, ω2=gλtanhλh, which may be real or imaginary values; P and G represent, respectively, the wave and group celerity as
(51)P=ωλG=12(1+2λhsinh⁡2λh)P,
and k is the wave number which could be expressed as
(52)k2=λ2(1-f2ω2).For real value ofλ, k is a real value when f<ω, 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)∂u∂t-fv=-g∂ς∂x∂v∂t+fu=-g∂ς∂y∂w∂t=-g∂ς∂z.In this section, Coriolis effects on long waves propagating in an infinite channel are investigated. The channel has two parallel side vertical walls of widthL, which is shown in Figure 29. Water depth is H0 and waves propagating from the left to the right.Figure 29
An infinite channel of widthL, with Coriolis parameter f.According to Pedlosky [75], the analytic solution in Kelvin wave could be expressed as
(54)η=Re[ηoe-fy/Coeiλx}u=η0H0C0e-fy/C0cos⁡(kx±ωt+φ),
where C0=(gH0)1/2; H0 is constant depth of the channel; g is the gravity; η0 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 withF(z)=cosh⁡λ(z+h)/cosh⁡λh, the fluid displacement, ς, can be expressed as
(55)ς(x,y,z;t)=η(x,y)F(z)e-iωtη(x,y)=A(x,y)ei∫kn·dx⃑=A(x,y)ei∫(kxdx+kydy),
where A is amplitude of free surface displacement, η; k is wave number; kx and ky are the components of k in x- and y-direction, respectively.Assuming that the flow velocities associated with wave action,(u,v), also have the same vertical distribution and with (U,V) as the components at the free surface:
(56)(u,v)=(U,V)F(z).From (53) and (55)-(56), the velocity components formula with Coriolis effects can be derived and expressed as
(57)U=gfω2-f2∂η∂y-igωω2-f2∂η∂xV=-gfω2-f2∂η∂x-igωω2-f2∂η∂y.Then from the no-flux condition,U→·n→=0, the solid boundary condition could be expressed as
(58)∂η∂n=±ifω∂η∂s,
where n and s 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
(59)∂ηS∂n-iληS=0.In order to illustrate the Coriolis effects of waves propagating in a channel with widthL=60 km and water depth H0=3 m have been used. A wave period T=46800 sec, amplitude = 1.0 m, and Coriolis parameter f=0.0000935638 sec−1 (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 c/f, the Rossby radius of deformation. The distribution of the wave velocity is parallel to the walls of the channel and also proportional to c/f.The results of analytical solutions: (a) distribution of wave amplitude; (b) distribution of wave velocity.
(a)
(b)For numerical computations, the incident wave amplitude is chosen to beηoe-fy/Co 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 x-direction between the numerical and analytical results at cross section x = 10,000 m are further compared as shown in Figure 32. Very good agreement between numerical and analytical results is observed. The velocity of y-direction is equal to 0 with a negligible error. Very good agreement between numerical results and analytical solutions is observed.The results of numerical model: (a) distribution of wave amplitude; (b) distribution of wave velocity.
(a)
(b)Comparisons of amplitude and velocity atx-diretion between numerical and analytical results at the cross section x = 10,000 m.
(a)
(b)
(c)


## 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
(43)∇·(CCg∇η)+k2CCgη=0,
where C=ω/k and Cg=(1/2)(1+2kh/sinh⁡2kh)C represent wave celerity and group celerity, respectively; k 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
(44)∂η∂n=0.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 Figure19, the depth between the shoreline and outer boundary of the shoal is described as
(45)h=h0(rr0)2,a≤r≤r0,
where r is the distance measured from the center of the island, a is radius of the circular island, r0 is the radius of paraboloidal shoal, and h0 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 fora≤r≤r0, derived by Homma [71], are presented as the following for completeness:
(46)η=a0∑n=0∞kεni2n[Jn′(kr0)Hn(kr0)-Jn(kr0)Hn′(lr0)]×(αn+1αn-1a-2αnrαn-1+r-αn-1)cos⁡nθ+[(αn+1αn-1)(αn+1)(a-2αnr0αn-2-r0-αn-2)Hn(kr0)-kHn′(kr0)(αn+1αn-1)a-2αnr0αn-1+r0-αn-1].θI is the incident angle of propagating waves and Jn and Hn are the Bessel and Hankel functions of the first kind of order n, as in the fourth benchmark problem. The superscript of the symbol of prime indicates their derivatives, respectively.Differentiating (47) with respect to radial (r) and angular (θ) variables, respectively, the analytical solutions of radial and tangential velocities in the inner region of a≤r≤r0 can be obtained analytically with any difficulty.A case for a circular island with radius,a=10 km, as shown in Figure 19, is considered. The water depth is constant, h0=4 km, at distance greater than r0=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 t=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 x- and y-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.Contours of wave amplitude lines around a circular island on a paraboloidal shoal,T=720 s. (a) Numerical result; (b) analytic results.
(a)
(b)Contours of wave phase lines around a circular island on a paraboloidal shoal,T=720 s, at t=0. (a) Numerical results; (b) analytic results.
(a)
(b)Contours of magnitude of velocity around a circular island on a paraboloidal shoal,T=720 s, t=0. (a) Numerical results; (b) analytic results.
(a)
(b)Numerical and analytic results of amplitude, phase, and magnitude of velocity around atr=a.
(a)
(b)
(c)


## 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 Figure24 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 Figure25. 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:
(47)ηH=a0ND∑n=0∞εn2πJ0(nδ)kaJn′(ka)Jn(kr)cos⁡n(θ-θI),
where ηH is water free-surface displacement within the harbor, a0 is incident wave amplitude, δ is half of the harbor entrance opening angle, θ0 is the angle from the positive x-coordinate to the centerline of entrance opening, θI is direction of incident waves with respect to positive x-axis, a is radius of the circular harbor, k is incident wave number, and N is weighted average of the forcing pressure:
(48)N=∑n=0∞2iπka2εnJ0(nδ)Hn′(ka)×(cos⁡nπ2cos⁡nθIcos⁡θ0+isinnπ2sinnθIsinnθ0),Figure 25
Definition of a circular harbor with protruding breakwaters.andD is term of the coefficient:
(49)D=∑n=0∞εn2πJn(ka)J02(nδ)kaJn′(ka)-∑n=0∞εn2πHn(ka)J02(nδ)kaHn′(ka)(1+cos⁡nπcos⁡2nθ0).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θI of 225° is presented. Present numerical results for an idealized circular harbor of radius a=0.75 ft with harbor mouth 60° opening (2δ).Response curves atr=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 ka=0.5 and the amplification factor, R, 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 ka=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 atr=0.7 ft, θ=0° of a circular harbor of radius, a=0.75 ft, 60° opening for θI=225°.Figure 27
Contours of wave amplitude of a circular harbor of radius,a=0.75 ft, 60° opening, with protruding breakwaters at ka=0.5 for θI=225°.Velocity vectors of a circular harbor of radius,a=0.75 ft, 60° opening, with protruding breakwaters at ka=0.5. (a) t=0, (b) t=T/2 for θI=225°.
(a)
(b)


## 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]:
(50)∇·(PG∇η)+λ2PG(1-f2ω2)η=0,
where f 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, ω2=gλtanhλh, which may be real or imaginary values; P and G represent, respectively, the wave and group celerity as
(51)P=ωλG=12(1+2λhsinh⁡2λh)P,
and k is the wave number which could be expressed as
(52)k2=λ2(1-f2ω2).For real value ofλ, k is a real value when f<ω, 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)∂u∂t-fv=-g∂ς∂x∂v∂t+fu=-g∂ς∂y∂w∂t=-g∂ς∂z.In this section, Coriolis effects on long waves propagating in an infinite channel are investigated. The channel has two parallel side vertical walls of widthL, which is shown in Figure 29. Water depth is H0 and waves propagating from the left to the right.Figure 29
An infinite channel of widthL, with Coriolis parameter f.According to Pedlosky [75], the analytic solution in Kelvin wave could be expressed as
(54)η=Re[ηoe-fy/Coeiλx}u=η0H0C0e-fy/C0cos⁡(kx±ωt+φ),
where C0=(gH0)1/2; H0 is constant depth of the channel; g is the gravity; η0 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 withF(z)=cosh⁡λ(z+h)/cosh⁡λh, the fluid displacement, ς, can be expressed as
(55)ς(x,y,z;t)=η(x,y)F(z)e-iωtη(x,y)=A(x,y)ei∫kn·dx⃑=A(x,y)ei∫(kxdx+kydy),
where A is amplitude of free surface displacement, η; k is wave number; kx and ky are the components of k in x- and y-direction, respectively.Assuming that the flow velocities associated with wave action,(u,v), also have the same vertical distribution and with (U,V) as the components at the free surface:
(56)(u,v)=(U,V)F(z).From (53) and (55)-(56), the velocity components formula with Coriolis effects can be derived and expressed as
(57)U=gfω2-f2∂η∂y-igωω2-f2∂η∂xV=-gfω2-f2∂η∂x-igωω2-f2∂η∂y.Then from the no-flux condition,U→·n→=0, the solid boundary condition could be expressed as
(58)∂η∂n=±ifω∂η∂s,
where n and s 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
(59)∂ηS∂n-iληS=0.In order to illustrate the Coriolis effects of waves propagating in a channel with widthL=60 km and water depth H0=3 m have been used. A wave period T=46800 sec, amplitude = 1.0 m, and Coriolis parameter f=0.0000935638 sec−1 (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 c/f, the Rossby radius of deformation. The distribution of the wave velocity is parallel to the walls of the channel and also proportional to c/f.The results of analytical solutions: (a) distribution of wave amplitude; (b) distribution of wave velocity.
(a)
(b)For numerical computations, the incident wave amplitude is chosen to beηoe-fy/Co 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 x-direction between the numerical and analytical results at cross section x = 10,000 m are further compared as shown in Figure 32. Very good agreement between numerical and analytical results is observed. The velocity of y-direction is equal to 0 with a negligible error. Very good agreement between numerical results and analytical solutions is observed.The results of numerical model: (a) distribution of wave amplitude; (b) distribution of wave velocity.
(a)
(b)Comparisons of amplitude and velocity atx-diretion between numerical and analytical results at the cross section x = 10,000 m.
(a)
(b)
(c)


## 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
(60)∇2ϕ=0,
where ϕ(x⃑,t) is the velocity potential with x⃑ being the position vector and t being the time. The velocity vector is defined as v⃑=∇ϕ ((∂ϕ/∂x)i⃑+(∂ϕ/∂y)j⃑ for 2-D problems and (∂ϕ/∂x)i⃑+(∂ϕ/∂y)j⃑+(∂ϕ/∂z)k⃑ for 3-D problems). The kinematic free-surface boundary condition (KFSBC) and dynamic free-surface boundary condition (DFSBC) can be expressed as
(61)∇ϕ=dx⃑dt,x⃑∈theentiredomain,iiiiiiiiiiiiiincludingthefreesurface,dϕdt=-gz+12∇ϕ·∇ϕ,x⃑∈thefreesurface,
where g is the gravitational acceleration, z is the vertical coordinate, and C(t) 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)n⃑·∇ϕ=n⃑·v⃑b,
where n⃑ is the unit normal vector outward from the body surface, and v⃑b 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
(63)(ϕ|x⃑=x⃑j)(n)=(ϕ|x⃑=x⃑j)(n-2)+2Δt[(-gz+12∇ϕ·∇ϕ)x⃑=x⃑j](n-1),
where x⃑j denotes the position of the jth 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 nth 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,” x⃑j(n). 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 x⃑j(n):
(64)x⃑j(n)=x⃑j(n-2)+2Δt(∇ϕ|x⃑=x⃑j)(n-1).In the case of strong nonlinearity, the free-surface particle trajectoryx⃑j(n) should be modified by applying the Crank-Nicolson formula for a better numerical stability:
(65)x⃑j(n)=x⃑j(n-1)+Δt2[(∇ϕ|x⃑=x⃑j)(n)+(∇ϕ|x⃑=x⃑j)(n-1)].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 atx⃑j(n) 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)dξdt=u-(ξ,t)=Cη|x=ξh+η|x=ξ,
where ξ(t) is the position of the wave paddle at time t, while t is elapsed time since the start of the motion, C is the wave celerity, η is the free surface displacement, and h is the still water depth. In [82] it was proposed to apply Boussinesq’s solitary wave solution [83] into (66). That is,
(67)η=HS2,
where H is the wave height, and
(68)S=cosh⁡-1[KX]
in which K is the boundary outskirt decay coefficient, and
(69)X=(ξ-Ct-x0),
where x0 is the initial position of the wave crest. Because the solitary wave is generated in a quiescent water flume, x0 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
(70)K=3H4h3,
where the wave celerity is determined as
(71)C=g(h+H).Figure33 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].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.
(a)
(b)
(c)
(d)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)η=h((Hh)η1+(Hh)2η2+⋯+(Hh)9η9).The formulae to determine the values ofη1 to η9 could be found in [84]. In an analogous way, the boundary outskirt decay coefficient and the wave celerity are determined as
(73)K=3H4h3(1+(Hh)K1+(Hh)2K2+⋯(Hh)9K9),C2=gh(1+(Hh)C1+(Hh)2C2+⋯+(Hh)9C9),
where the formulae to determine K1,…,K9 and C1,…,C9 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 ofH0=0.14 m and h=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 Figure36.Figure 36
Initial nodal distribution of the numerical model of the liquid sloshing in a swaying tank.Figure37 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 t=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 t=6.0T~6.5T. 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:
(74)t=6.0T~6.5T.Comparison of the numerical results with the experimental data.
(a)
(b)
(c)Snapshots of traced fluid particles in the time interval oft=6.0T~6.5T.
(a)
(b)
(c)
(d)
(e)
(f)


## 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)dξdt=u-(ξ,t)=Cη|x=ξh+η|x=ξ,
where ξ(t) is the position of the wave paddle at time t, while t is elapsed time since the start of the motion, C is the wave celerity, η is the free surface displacement, and h is the still water depth. In [82] it was proposed to apply Boussinesq’s solitary wave solution [83] into (66). That is,
(67)η=HS2,
where H is the wave height, and
(68)S=cosh⁡-1[KX]
in which K is the boundary outskirt decay coefficient, and
(69)X=(ξ-Ct-x0),
where x0 is the initial position of the wave crest. Because the solitary wave is generated in a quiescent water flume, x0 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
(70)K=3H4h3,
where the wave celerity is determined as
(71)C=g(h+H).Figure33 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].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.
(a)
(b)
(c)
(d)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)η=h((Hh)η1+(Hh)2η2+⋯+(Hh)9η9).The formulae to determine the values ofη1 to η9 could be found in [84]. In an analogous way, the boundary outskirt decay coefficient and the wave celerity are determined as
(73)K=3H4h3(1+(Hh)K1+(Hh)2K2+⋯(Hh)9K9),C2=gh(1+(Hh)C1+(Hh)2C2+⋯+(Hh)9C9),
where the formulae to determine K1,…,K9 and C1,…,C9 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 ofH0=0.14 m and h=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 Figure36.Figure 36
Initial nodal distribution of the numerical model of the liquid sloshing in a swaying tank.Figure37 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 t=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 t=6.0T~6.5T. 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:
(74)t=6.0T~6.5T.Comparison of the numerical results with the experimental data.
(a)
(b)
(c)Snapshots of traced fluid particles in the time interval oft=6.0T~6.5T.
(a)
(b)
(c)
(d)
(e)
(f)


## 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.



---

*Source: 350364-2014-07-24.xml*