# The Numerical Analysis of the Flow on the Smooth and Nonsmooth Boundaries by IBEM/DBIEM

Authors: Gang Xu; Guangwei Zhao; Jing Chen; Shuqi Wang; Weichao Shi
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The value of the tangential velocity on the Boundary Value Problem (BVP) is inaccurate when comparing the results with analytical solutions by Indirect Boundary Element Method (IBEM), especially at the intersection region where the normal vector is changing rapidly (named nonsmooth boundary). In this study, the singularity of the BVP, which is directly arranged in the center of the surface of the fluid computing domain, is moved outside the computational domain by using the Desingularized Boundary Integral Equation Method (DBIEM). In order to analyze the accuracy of the IBEM/DBIEM and validate the above-mentioned problem, three-dimensional uniform flow over a sphere has been presented. The convergent study of the presented model has been investigated, including desingularized distance in the DBIEM. Then, the numerical results were compared with the analytical solution. It was found that the accuracy of velocity distribution in the flow field has been greatly improved at the intersection region, which has suddenly changed the boundary surface shape of the fluid domain. The conclusions can guide the study on the flow over nonsmooth boundaries by using boundary value method.


---

## Body

## 1. Introduction

In the field of interaction between fluid and structures, Finite Element Method (FEM) and Boundary Element Method (BEM) [1, 2] have been applied to predict the hydrodynamic force acting on the ocean structures. The FEM is usually used to solve the Boundary Value Problem (BVP). Compared to the Boundary Element Method (BEM), the coefficient matrix of the FEM is banded and has a smaller memory requirement. Ma et al. [3], Hu et al. [4], and Wang and Wu [5, 6] obtain a variety of fully nonlinear and weakly nonlinear problems solution by FEM. However, the major challenge with the FEM is the mesh generation. For the flow over complex structures, a complicated mesh generator is commonly required to follow the motion of the structure and the fluid flow. Hence, the BEM to solve the BVP is employed in this work.BEM has been widely used to solve the BVP in potential theory since Jaswon [7] and Symm [8]. In the BEM, mixed distribution of sources and dipoles integral equations and sources, only distribution in the integral equations, are widely used, where the former is named Direct Boundary Element Method (DBEM) as it originates from Green’s third formula [9–11], while the latter is often called Indirect Boundary Element Method (IBEM) as it comes from the singularity analysis [12]. Hess and Smith [13] first introduced this method as the solution of the hydrodynamic problems. However, it possesses some fundamental defects which limit its application to the fluid’s and structure’s interaction problems. DBEM is superior to IBEM from the view of accuracy for solving the velocity on the boundary surface as IBEM cannot solve the induced velocity with satisfactory accuracy for nonsmooth boundaries. Newman and Lee [14] and Choi et al. [15] show that the nonlinear force acting on a box is sensitive to the panel resolution near sharp edges by IBEM. Compared to the time-domain prediction based on the DBEM/IBEM, the frequency domain approaches [16] have been widely adopted to solve hydrodynamic problems. Some commercial software, such as WAMIT [17] and HYDROSTAR [18], can evaluate in practical offshore engineering within the frequency domain. This however cannot make it suitable for nonlinear problems. In time-domain prediction, many researchers prefer to simulate problems using the Transient Green Function (TGF) [19, 20] and the Rankine Source Method (RSM) [21, 22]. Although the RSM is simple to evaluate, it requires the discretization of the surfaces of the fluid domain. An artificial damping zone [23, 24] should also be given to absorb the wave energy. TGF satisfies the linear free surface boundary condition and the far field radiation condition. However, the TGF is complex and difficult to obtain, especially for the nonlinear flow problem.In order to improve the computation accuracy of IBEM for the solution of nonsmooth boundary problems, the Desingularized Boundary Integral Equation Method (DBIEM) is used/ will be employed to analyze the fluid flow problems, which has been used previously in solving wave-structure interaction problems such as in the work by Zhang et al. [25, 26], Beck [27], Kim et al. [28], Celebi [29], Kara et al. [30], and Xu et al. [31]. Compared with the DBEM/IBEM, the integral kernels of the DBIEM are no longer singular as the singularities are placed slightly outside the fluid domain. This is particularly advantageous when the direct differentiation is applied to the integral equation to obtain the velocity by IBEM. In this work, the DBIEM is selected to obtain high accuracy velocity distribution for simulating nonsmooth boundary problems. The main advantage of DBIEM, compared with FEM, also lies in only having to discretize the surfaces of the fluid domain. When the boundary of the fluid domain is confined and the number of the discretized elements is limited, the DBIEM may offer a better computational efficiency and less memory requirement, even its matrix is fully populated.This paper shows the numerical analysis of the IBEM/DBIEM for solving smooth and nonsmooth boundary problems. Mathematical formulation of DBEM, IBEM, and DBIEM is reviewed in Section2. The numerical implementations for the methods are discussed in Section 3. The numerical results and discussions are presented in Section 4.


## 2. Mathematical Formulation

### 2.1. Direct Boundary Element Method (DBEM)

As a common assumption in the potential flow theory, the fluid is assumed to be incompressible, inviscid, and flow irrotational. Then, the fluid motion can be described by the velocity potentialϕ and the fluid velocityv=▽ϕ, which is governed by the Laplace equation in the fluid domainD as shown in Figure 1; the equation for the velocity potentialϕ becomes(1)∇2ϕ=0inD,Based on Gauss formula and three-dimensional Green’s third formula, we can find the valueϕ(p) by DBEM in the flow domainD andϕ(p) on the surfaceS=S1+S2+S3+S4+S5+S6 of the flow domainD by(2)Cpϕp=∬SGp,q∂ϕq∂nq-ϕq∂∂nqGp,qdsq(3)Gp,q=1rHere,p(x,y,z) is a field point andq(ξ,η,ζ) is a source point on the surfaceS of the fluid domainD. The geometrical coefficientC(p) =2π when the field point is located on the surfaceS; otherwise,C(p) =4π.G(p,q) is a Green’s function, wherer is the distance between field pointp and source pointq.Figure 1
The sketch of coordinate system and fluid domain.


### 2.2. Indirect Boundary Element Method (IBEM) to Solve a Uniform Flow over a Sphere

If we consider a uniform flow over a sphere, we only need to use half of the domain to represent a uniform flow over the whole sphere as long as we set reasonable boundary condition on the surfaces of the half domain. The reference system of Cartesian coordinates is defined by letting (x,y) plane coincide with the symmetry surface of the sphere and parallel with flow velocityU,z points vertically upwards from the symmetry surface of the sphere as shown in Figure 2. The sphere surface is denoted bySs and its unit normal vector directed outward from the fluid region is denoted by n. The boundary surfaceS of the Model (a) and Model (b) isS=Ss +Sa-b + Sc-d + Sd-e +Se-f + Sf-a.Figure 2
Computational domain and Green’s function for the potential flow model.
(a)
G without symmetric condition (b)
G with symmetric conditionBased on the potential flow theory, the velocity potentialϕ satisfies the Laplace equation and the velocity potential is subject to the following boundary conditions on the boundary surfaces as shown in Model (a), where the Green’s function does not consider the symmetric condition over symmetric plane surfaceSa-b-c-d:(4)∂ϕ∂n=Uonthesurfaced-e(5)∂ϕ∂n=-Uonthesurfacef-a(6)∂ϕ∂n=0onthesurfacea-b-c-dande-f(7)∂ϕ∂n=0onthesurfaceofthesphereSsUsing the Direct Boundary Element Method (DBEM) in Section 2.1, we have(8)ϕp=14π∬sGp,q∂ϕq∂nq-ϕq∂∂nqGp,qdsqintheflowdomainand the source strengthσ(9)∬sσq∂Gp,q∂npdsq=∂ϕp∂npThen we can use Indirect Boundary Element Method (IBEM) to find the velocity distribution of the fluid domain by(10)∇∬sσqGp,qdsq=∇ϕpHere, the Green’s functionG for the Model (a), which does not consider the symmetry surface, satisfies(11)Gp,q=1r1p,qFrom (4)-(11), we use half domain simulation to represent the uniform flow over the whole sphere as mentioned above. For Model (a) in Figure 2, the mesh discretization on the surfaces of the half domain is needed, so Model (a) has a mutation of the normal vector at the junction of the hemispherical and symmetry surfaceSa-b-c-d. This means the boundary at junction is nonsmooth. Model (b) however uses the symmetric condition of Green’s function G on the symmetry surfaceSa-b-c-d. Thus the mesh’s discretization on the symmetry surfacesSa-b and Sc-d of Model (b) is not needed and there is no normal vector mutation at the junction of hemispherical and symmetry surface. This means the boundary at the fictitious junction is smooth. Here, the Green’s function G satisfies(12)Gp,q,q′=1r1p,q+1r2p,q′where(13)r1p,q=x-ξ2+y-η2+z-ζ2(14)r2p,q′=x-ξ2+y-η2+z+ζ2Here,p (x,y,z) is the field point in the considered domain, and the source pointq’(ξ,η,-ζ) is the mirror point of the source pointq(ξ,η,ζ) on the discretized surfaceS using thexoy plane as the symmetry surface.


### 2.3. Desingularized Boundary Integral Element Method (DBIEM) to Solve a Uniform Flow over a Sphere

In this study, the unknown velocity potentialϕ also can be obtained by the Desingularized Boundary Integral Element Method (DBIEM) for the Boundary Value Problem (BVP). This method obtains the solution by distributing Rankine sources over the surfaceS outside the fluid domainD. This surface is at a small distance away from the corresponding real boundary of the fluid domain, as presented in Figure 3. The velocity potential in the fluid domainD can be written as follows:(15)ϕp=∬SσqGp,qdsqwhereq(ξ,η,ζ) is the integration point on the integration surfaceS outside the fluid domain,p(x,y,z) is the field point where the potential is evaluated,σ is the unknown source strength distribution over the surfaceS.Figure 3
Schematic of isolated source and node locations.For the problem considered in this work, we construct the solution using a constant-strength source point within each mesh element over the integration boundary. Applying the relevant boundary conditions, the desingularized boundary integral equations that must be solved to determine the unknown source strengths are(16)∬SσxsGp,qdΩ=ϕ0xcxc∈Γd(17)∬Sσxs∂Gp,q∂ndΩ=χxcxc∈Γnwherexs=(ξ,η,ζ): a source point on the integration surface,xc=(x,y,z): a field point on the real boundary,ϕ0=the given velocity potential value atxc,Гd=surface on whichϕ0 is given,χ=the given normal velocity of velocity potential atxc,Гn =surface on whichχ is given,G(p,q)=1/r(p,q), wherer is the distance between field pointp and the source pointq.In the desingularized method, the source distribution is outside the fluid domain so that the source points never coincide with the field points and therefore the integrals are nonsingular. In addition, because of the desingularization, we can use simple isolated Rankine sources and obtain the equivalent accuracy. This greatly reduces the complexity of the form of the influence coefficients that make up the elements of the kernel matrix. Then the integral equations in (16) and (17) can be replaced by a discrete summation ofN-isolated singularities located at a small distance away from the corresponding control point on the boundaries,(18)∑i=1NσqiGp,qi=ϕ0p(19)∑i=1Nσqi∂∂npGp,qi=∂ϕ0p∂npThe desingularized distance between the/an isolated source point and corresponding control point is given by(20)Ld=ldDmβwhereld and β are constants andDm is a measure of the local mesh size (typically the square root of the local mesh area). The accuracy and convergence of the solutions are sensitive to the choices ofld and β. Therefore, appropriateld and β values need to be determined after numerical investigation. The recommended values areld =0.5-1.0 and β=0.5. A detailed study for three-dimensional flow with regard to the performance of DBIEM with the desingularization parameters (Cao et al. [28]) has been done in this work.Once the above integral equations using isolated Rankine source are solved, the fluid velocity in (18) can be calculated from direct derivatives,(21)V=∇ϕp=∑i=1Nσqi∇Gp,qi


## 2.1. Direct Boundary Element Method (DBEM)

As a common assumption in the potential flow theory, the fluid is assumed to be incompressible, inviscid, and flow irrotational. Then, the fluid motion can be described by the velocity potentialϕ and the fluid velocityv=▽ϕ, which is governed by the Laplace equation in the fluid domainD as shown in Figure 1; the equation for the velocity potentialϕ becomes(1)∇2ϕ=0inD,Based on Gauss formula and three-dimensional Green’s third formula, we can find the valueϕ(p) by DBEM in the flow domainD andϕ(p) on the surfaceS=S1+S2+S3+S4+S5+S6 of the flow domainD by(2)Cpϕp=∬SGp,q∂ϕq∂nq-ϕq∂∂nqGp,qdsq(3)Gp,q=1rHere,p(x,y,z) is a field point andq(ξ,η,ζ) is a source point on the surfaceS of the fluid domainD. The geometrical coefficientC(p) =2π when the field point is located on the surfaceS; otherwise,C(p) =4π.G(p,q) is a Green’s function, wherer is the distance between field pointp and source pointq.Figure 1
The sketch of coordinate system and fluid domain.


## 2.2. Indirect Boundary Element Method (IBEM) to Solve a Uniform Flow over a Sphere

If we consider a uniform flow over a sphere, we only need to use half of the domain to represent a uniform flow over the whole sphere as long as we set reasonable boundary condition on the surfaces of the half domain. The reference system of Cartesian coordinates is defined by letting (x,y) plane coincide with the symmetry surface of the sphere and parallel with flow velocityU,z points vertically upwards from the symmetry surface of the sphere as shown in Figure 2. The sphere surface is denoted bySs and its unit normal vector directed outward from the fluid region is denoted by n. The boundary surfaceS of the Model (a) and Model (b) isS=Ss +Sa-b + Sc-d + Sd-e +Se-f + Sf-a.Figure 2
Computational domain and Green’s function for the potential flow model.
(a)
G without symmetric condition (b)
G with symmetric conditionBased on the potential flow theory, the velocity potentialϕ satisfies the Laplace equation and the velocity potential is subject to the following boundary conditions on the boundary surfaces as shown in Model (a), where the Green’s function does not consider the symmetric condition over symmetric plane surfaceSa-b-c-d:(4)∂ϕ∂n=Uonthesurfaced-e(5)∂ϕ∂n=-Uonthesurfacef-a(6)∂ϕ∂n=0onthesurfacea-b-c-dande-f(7)∂ϕ∂n=0onthesurfaceofthesphereSsUsing the Direct Boundary Element Method (DBEM) in Section 2.1, we have(8)ϕp=14π∬sGp,q∂ϕq∂nq-ϕq∂∂nqGp,qdsqintheflowdomainand the source strengthσ(9)∬sσq∂Gp,q∂npdsq=∂ϕp∂npThen we can use Indirect Boundary Element Method (IBEM) to find the velocity distribution of the fluid domain by(10)∇∬sσqGp,qdsq=∇ϕpHere, the Green’s functionG for the Model (a), which does not consider the symmetry surface, satisfies(11)Gp,q=1r1p,qFrom (4)-(11), we use half domain simulation to represent the uniform flow over the whole sphere as mentioned above. For Model (a) in Figure 2, the mesh discretization on the surfaces of the half domain is needed, so Model (a) has a mutation of the normal vector at the junction of the hemispherical and symmetry surfaceSa-b-c-d. This means the boundary at junction is nonsmooth. Model (b) however uses the symmetric condition of Green’s function G on the symmetry surfaceSa-b-c-d. Thus the mesh’s discretization on the symmetry surfacesSa-b and Sc-d of Model (b) is not needed and there is no normal vector mutation at the junction of hemispherical and symmetry surface. This means the boundary at the fictitious junction is smooth. Here, the Green’s function G satisfies(12)Gp,q,q′=1r1p,q+1r2p,q′where(13)r1p,q=x-ξ2+y-η2+z-ζ2(14)r2p,q′=x-ξ2+y-η2+z+ζ2Here,p (x,y,z) is the field point in the considered domain, and the source pointq’(ξ,η,-ζ) is the mirror point of the source pointq(ξ,η,ζ) on the discretized surfaceS using thexoy plane as the symmetry surface.


## 2.3. Desingularized Boundary Integral Element Method (DBIEM) to Solve a Uniform Flow over a Sphere

In this study, the unknown velocity potentialϕ also can be obtained by the Desingularized Boundary Integral Element Method (DBIEM) for the Boundary Value Problem (BVP). This method obtains the solution by distributing Rankine sources over the surfaceS outside the fluid domainD. This surface is at a small distance away from the corresponding real boundary of the fluid domain, as presented in Figure 3. The velocity potential in the fluid domainD can be written as follows:(15)ϕp=∬SσqGp,qdsqwhereq(ξ,η,ζ) is the integration point on the integration surfaceS outside the fluid domain,p(x,y,z) is the field point where the potential is evaluated,σ is the unknown source strength distribution over the surfaceS.Figure 3
Schematic of isolated source and node locations.For the problem considered in this work, we construct the solution using a constant-strength source point within each mesh element over the integration boundary. Applying the relevant boundary conditions, the desingularized boundary integral equations that must be solved to determine the unknown source strengths are(16)∬SσxsGp,qdΩ=ϕ0xcxc∈Γd(17)∬Sσxs∂Gp,q∂ndΩ=χxcxc∈Γnwherexs=(ξ,η,ζ): a source point on the integration surface,xc=(x,y,z): a field point on the real boundary,ϕ0=the given velocity potential value atxc,Гd=surface on whichϕ0 is given,χ=the given normal velocity of velocity potential atxc,Гn =surface on whichχ is given,G(p,q)=1/r(p,q), wherer is the distance between field pointp and the source pointq.In the desingularized method, the source distribution is outside the fluid domain so that the source points never coincide with the field points and therefore the integrals are nonsingular. In addition, because of the desingularization, we can use simple isolated Rankine sources and obtain the equivalent accuracy. This greatly reduces the complexity of the form of the influence coefficients that make up the elements of the kernel matrix. Then the integral equations in (16) and (17) can be replaced by a discrete summation ofN-isolated singularities located at a small distance away from the corresponding control point on the boundaries,(18)∑i=1NσqiGp,qi=ϕ0p(19)∑i=1Nσqi∂∂npGp,qi=∂ϕ0p∂npThe desingularized distance between the/an isolated source point and corresponding control point is given by(20)Ld=ldDmβwhereld and β are constants andDm is a measure of the local mesh size (typically the square root of the local mesh area). The accuracy and convergence of the solutions are sensitive to the choices ofld and β. Therefore, appropriateld and β values need to be determined after numerical investigation. The recommended values areld =0.5-1.0 and β=0.5. A detailed study for three-dimensional flow with regard to the performance of DBIEM with the desingularization parameters (Cao et al. [28]) has been done in this work.Once the above integral equations using isolated Rankine source are solved, the fluid velocity in (18) can be calculated from direct derivatives,(21)V=∇ϕp=∑i=1Nσqi∇Gp,qi


## 3. Numerical Implementation

### 3.1. Numerical Implementation for DBEM

For the DBEM, (3) is solved using/applying the following numerical procedure, in which the surfaceSis discretized into finite numbers of facets. The corresponding values of ∂ϕ/∂n andϕ are taken as constant over each facet and applied at the facet center. Equation (3) is thereby represented by the following set of simultaneous equations:(22)∑j=1nS1+nS2+nS3+nS4+nS5+nS6Dijϕj+∑j=1nS1+nS2+nS3+nS4+nS5+nS6Sij∂ϕ∂nj=0i=1,2,…,nS1+nS2+nS3+nS4+nS5+nS6In whichnS are the numbers of facets on domain surface, respectively.The matrix coefficientsSij and Dij correspond to integrals of the Green’s function and its normal derivative over the area ∆Q of thejth facet, respectively. Sij and Dij are written as(23)Sij=∫ΔQjGijdsj,(24)Dij=∫ΔQj∂Gij∂njdsj,i≠j;2π,i=j.


### 3.2. Numerical Implementation for IBEM

In the IBEM, the source strengthσ of (9) is solved by the following numerical implementation for Model (a) and Model (b), respectively:Model (a)  (25) ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D i j σ j = ∂ ϕ ∂ n i  i = 1,2 , … ,   n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - aModel (b)  (26) ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D i j σ j = ∂ ϕ ∂ n i  i = 1,2 , … ,   n S s + n S d - e + n S e - f + n S f - aThen the velocity distribution for thex,y, andz direction in the flow field can be obtained byModel (a)  (27) ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D D i j x σ j = ∂ ϕ ∂ x i ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D D i j y σ j = ∂ ϕ ∂ y i ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D D i j z σ j = ∂ ϕ ∂ z i i = 1,2 , … ,   n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - aModel (b)  (28) ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D D i j x σ j = ∂ ϕ ∂ x i ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D D i j y σ j = ∂ ϕ ∂ y i ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D D i j z σ j = ∂ ϕ ∂ z i i = 1,2 , … ,   n S s + n S d - e + n S e - f + n S f - aIn (25)-(28),nS are the numbers of facets on domain surface. The matrix coefficientsDij correspond to the normal derivative of Green’s function over the area ∆Q of thejth face, as presented in (24) and the matrix coefficientsDDij are presented in (29).(29)DDijx,DDijy,DDijz=∫ΔQj∂Gij∂xjdsj,∫ΔQj∂Gij∂yjdsj,∫ΔQj∂Gij∂zjdsj


### 3.3. Numerical Implementation for DBIEM

For the flow over sphere, all vertical velocities on the boundary surfaces are known. Then only (19) needs to be solved in DBIEM by the following numerical implementation:(30)∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijσj=∂ϕ∂nii=1,2,…,nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aand the velocity can be obtained by(31)∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijxσj=∂ϕ∂xi∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijyσj=∂ϕ∂yi∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijzσj=∂ϕ∂zii=1,2,…,nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aThe matrix coefficientsNDij are presented in (32).(32)NDijx,NDijy,NDijz=∂Gij∂xj,∂Gij∂yj,∂Gij∂zj


## 3.1. Numerical Implementation for DBEM

For the DBEM, (3) is solved using/applying the following numerical procedure, in which the surfaceSis discretized into finite numbers of facets. The corresponding values of ∂ϕ/∂n andϕ are taken as constant over each facet and applied at the facet center. Equation (3) is thereby represented by the following set of simultaneous equations:(22)∑j=1nS1+nS2+nS3+nS4+nS5+nS6Dijϕj+∑j=1nS1+nS2+nS3+nS4+nS5+nS6Sij∂ϕ∂nj=0i=1,2,…,nS1+nS2+nS3+nS4+nS5+nS6In whichnS are the numbers of facets on domain surface, respectively.The matrix coefficientsSij and Dij correspond to integrals of the Green’s function and its normal derivative over the area ∆Q of thejth facet, respectively. Sij and Dij are written as(23)Sij=∫ΔQjGijdsj,(24)Dij=∫ΔQj∂Gij∂njdsj,i≠j;2π,i=j.


## 3.2. Numerical Implementation for IBEM

In the IBEM, the source strengthσ of (9) is solved by the following numerical implementation for Model (a) and Model (b), respectively:Model (a)  (25) ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D i j σ j = ∂ ϕ ∂ n i  i = 1,2 , … ,   n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - aModel (b)  (26) ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D i j σ j = ∂ ϕ ∂ n i  i = 1,2 , … ,   n S s + n S d - e + n S e - f + n S f - aThen the velocity distribution for thex,y, andz direction in the flow field can be obtained byModel (a)  (27) ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D D i j x σ j = ∂ ϕ ∂ x i ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D D i j y σ j = ∂ ϕ ∂ y i ∑ j = 1 n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - a D D i j z σ j = ∂ ϕ ∂ z i i = 1,2 , … ,   n S s + n S a - b + n S c - d + n S d - e + n S e - f + n S f - aModel (b)  (28) ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D D i j x σ j = ∂ ϕ ∂ x i ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D D i j y σ j = ∂ ϕ ∂ y i ∑ j = 1 n S s + n S d - e + n S e - f + n S f - a D D i j z σ j = ∂ ϕ ∂ z i i = 1,2 , … ,   n S s + n S d - e + n S e - f + n S f - aIn (25)-(28),nS are the numbers of facets on domain surface. The matrix coefficientsDij correspond to the normal derivative of Green’s function over the area ∆Q of thejth face, as presented in (24) and the matrix coefficientsDDij are presented in (29).(29)DDijx,DDijy,DDijz=∫ΔQj∂Gij∂xjdsj,∫ΔQj∂Gij∂yjdsj,∫ΔQj∂Gij∂zjdsj


## 3.3. Numerical Implementation for DBIEM

For the flow over sphere, all vertical velocities on the boundary surfaces are known. Then only (19) needs to be solved in DBIEM by the following numerical implementation:(30)∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijσj=∂ϕ∂nii=1,2,…,nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aand the velocity can be obtained by(31)∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijxσj=∂ϕ∂xi∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijyσj=∂ϕ∂yi∑j=1nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aNDijzσj=∂ϕ∂zii=1,2,…,nSs+nSa-b+nSc-d+nSd-e+nSe-f+nSf-aThe matrix coefficientsNDij are presented in (32).(32)NDijx,NDijy,NDijz=∂Gij∂xj,∂Gij∂yj,∂Gij∂zj


## 4. Numerical Results and Discussions

### 4.1. IBEM to Simulate the Uniform Flow over Sphere

Compared to the main dimension of the sphere, the computation domain is large enough, as shown in Figure4. The flow velocity of the left boundary surface (as presented in Figure 2) of the domain is set to beU (U=0.025m/s), and the velocity of the right boundary surface is alsoU. Based on the problem of uniform flow over the/a sphere (r=1.0m), the Green function of Model (a) satisfiesG=1/r1 and Model (b) satisfiesG=1/r1+1/r2.Figure 4
Computation domain and mesh of Model (a) (left) and Model (b) (right).We need to find the accuracy of the velocity distribution on the smooth/nonsmooth boundaries. The IBEM here has been used to solve the flow around the sphere whether the normal vector of the mesh has a mutation at the symmetry surface by using the symmetric condition in Green’s function. It means that the model has a mutation of the normal vector on the hemisphere without the use of the symmetric condition in Green’s function on the symmetry surface.Figures5–8 show the relative error with the analytical solution at different positions between them. As can be seen from the figure, the result of symmetric condition used is obviously closer to the analytical solution. Velocity distribution and the relative error of velocity vary little with the water depthz, and the velocity inx andz direction obtained without symmetric condition is very pronounced along with the water depth, especially for the velocity near the (x,y) plane, where the normal vector has tremendous mutation. It shows that, using the distributed source method, the accuracy of the tangential velocity obtained atz≈0 is poor and has large error. Figure 5 shows the relationship between the velocity potential of the sphere on thexoz plane with the change of coordinatez (x>0); Figure 6 shows the relation between the velocity inx direction of the sphere on thexoz plane with the change of coordinatez(x>0); Figure 7 shows the relation between the velocity inz direction of the sphere on thexozplane and the change of the coordinatez(x>0). Figure 8 shows the relative error of the velocity inxandzon thexoz plane regarding the mesh density. The curve trend of the figures indicates the increase in solution accuracy of the velocity on the sphere through increasing the mesh number, but this will reduce the computational efficiency of the IBEM.Figure 5
The velocity potential and its relative error.Figure 6
The velocity inx direction and its relative error.Figure 7
The velocity inz direction and its relative error.Figure 8
The relative error of velocity inx andz directions regarding the mesh density.


### 4.2. IBEM with Galerkin Technology to Simulate the Uniform Flow over Sphere

If the symmetric condition in Green’s function has been used and the IBEM method cannot obtain high accuracy velocity distribution on the surface of the sphere where the boundary surface is not smooth, the distribution source method based on Galerkin technology [29] has been presented to study the distribution source method. The computational domain and mesh are the same as in Section 4.1. The distributed source method based on Galerkin technology shows(33)ϕi=1ai∑j=1N∬Sidsi∬SjGσjdsj(34)2πσi+1ai∑j=1N∬Sidsi∬Sj∂G∂niσjdsj=VniNote that the Green’s functionG is similar to the above section. Since we will use quadrilateral mesh element same as before, the solution of (33) and (34) will be solved by the four point Gauss integral method. First, the weights of four Gauss points and their Jacobi determinants for each surface element are obtained, respectively, and then the influence coefficients of four Gauss points and the source points on each unit are multiplied, respectively. The weight coefficient and final summation can get the influence coefficient between the panels.Figures9–11 show whether it is based on the Galerkin technology. Whether the sphere uses the symmetric condition, the velocity and the relative error curve between the analytical solutions and the different positions at thexoz plane (x>0) have also been presented. It can be seen from the figures that the distribution source method based on Galerkin technology can obviously improve the precision of the velocity potential in most positions. The relative error of the numerical results using Galerkin technology is smaller than that of the result only using symmetric condition, but the velocity results after the use of symmetric condition are obviously closer to the analytical solution. No matter whether the Galerkin technology is used or not, the velocity inx andz direction is very intense withz when the symmetric condition has not been considered, especially for the positions where there has been normal vector mutation of the meshes, but the Galerkin technology can also improve the accuracy of velocity on the surface of the sphere, where there has not been normal vector mutation. Figure 9 shows the relationship between the velocity potential of the sphere on thexoz plane with the change of coordinatez  (x>0); Figure 10 shows the relation between the velocity inx direction of the sphere on thexoz plane with the change of coordinatez  (x>0); Figure 11 shows the relation between the velocity inz direction of the sphere on thexoz plane with the change of the coordinatez  (x>0).Figure 9
The velocity potential and its relative error.Figure 10
The velocity inx direction and its relative error.Figure 11
The velocity inz direction and its relative error.


### 4.3. DBIEM to Simulate the Uniform Flow over Sphere

In this section, DBIEM is further used to study the flow over the sphere and to improve the accuracy of the velocity distribution of the flow. Here, the symmetric condition in Green’s function has not been considered.


#### 4.3.1. The Influence of the Parameterβ on the Numerical Results

Figures12 and 13 and Tables 1 and 2 give the relation between the velocity inx,z direction in thexoz plane when β=0.25, 0.5, 1.0 and 1.5, and the relative error. It can be seen from the figures that the relative error of the velocity inxandz direction at β=0.25 and 0.5 changes little with the water depthz, and the relative error increases gradually when β=1.0 and 1.5, especially at the position near the (x,y) plane (i.e., the normal vector mutation). From the above study, it is found that β=0.25 and 0.5 can get better results for the problem of flow around the sphere, but in order to ensure that the desingularized point does not coincide with the corresponding control point of the mesh, this work proposes to take β=0.5.Table 1
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition β=0.25 β=0.5 β=1.0 β=1.5 -0.08682 226.37 5.35 12.6 14.19 124 225 -0.25783 25.02 5.37 2.62 4.33 11.3 57.3 -0.42101 13.89 5.41 3.82 5.05 12.5 38.5 -0.57139 11.43 5.47 2.9 4.52 10.6 32.6 -0.70442 10.63 5.56 3.48 4.65 10.1 30.3 -0.81603 10.33 5.68 2.59 4.28 9.71 29.7 -0.90286 10.2 5.8 3.55 4.63 9.64 30 -0.96225 9.97 5.98 1.5 3.17 9.83 31 -0.98987 8.58 -- 2.91 4.21 9.98 31Table 2
The relative error of velocity inz direction (%).
Z/m without symmetric condition with symmetric condition β=0.25 β=0.5 β=1.0 β=1.5 -0.08682 226.28 5.42 12.6 14.26 124 225 -0.25783 25.09 5.46 2.72 4.43 11.4 57.3 -0.42101 14.03 5.55 3.97 5.19 12.7 38.6 -0.57139 11.66 5.71 3.15 4.77 10.8 32.8 -0.70442 11.03 5.97 3.92 5.08 10.5 30.6 -0.81603 11.08 6.44 3.41 5.08 10.5 30.3 -0.90286 11.8 7.4 5.26 6.32 11.2 31.3 -0.96225 14.45 10.28 6.4 7.98 14.3 34.4 -0.98987 44.06 41.38 40.6 41.39 44.9 57.8Figure 12
The velocity inx direction and its relative error.Figure 13
The velocity inz direction and its relative error.


#### 4.3.2. The Influence of the Parameterld on the Numerical Results

Figures14 and 15 and Tables 3 and 4 give the variation of the velocity inx,z direction atld=0.1, 0.5, 1.0, 1.5, 2.0, 2.5 and 3, respectively, and the relative error. It can be seen from the figures that the relative error of the velocity inx andz direction is huge when the integration point and control point coincide withld =0.1 (ld close to zero). Asld =0.5-2.0 gradually increases, the relative error decreases rapidly, but the relative error at the symmetry surface of sphere has a sudden change, with the increase to 2.0-3.0, due to the singular distance. Its relative error increases rapidly. It is found that the relative error is smaller whenld =1.5, and the change near the symmetry surface of sphere is not obvious. Therefore, this work proposes to selectld =1.5 for correlation calculation and analysis.Table 3
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition l d=0.1 l d=0.5 l d=1.0 l d=1.5 l d=2.0 l d=2.5 l d=3.0 -0.08682 226.37 5.35 104 86.7 14.19 5.57 42.62 131.08 110 -0.25783 25.02 5.37 98.9 32.1 4.33 0.58 0.71 19.99 124 -0.42101 13.89 5.41 98 34.3 5.05 1.24 1.19 10.12 121 -0.57139 11.43 5.47 97.7 33.1 4.52 0.66 0.68 13.63 117 -0.70442 10.63 5.56 97.4 33 4.65 0.98 0.05 10.95 113 -0.81603 10.33 5.68 97.3 32.7 4.28 0.35 1.1 11.91 109 -0.90286 10.2 5.8 97.2 32.7 4.63 1.03 0.04 9.20 106 -0.96225 9.97 5.98 97.1 32.2 3.17 0.83 2.24 11.09 103 -0.98987 8.58 -- 97.1 32.5 4.21 0.48 0.76 8.7 102Table 4
The relative error of velocity in z direction (%).
Z/m without symmetric condition with symmetric condition l d=0.1 l d=0.5 l d=1.0 l d=1.5 l d=2.0 l d=2.5 l d=3.0 -0.08682 226.28 5.42 104.0 86.7 14.26 5.64 42.67 144.78 110 -0.25783 25.09 5.46 98.9 32.2 4.43 0.68 0.61 18.30 124 -0.42101 14.03 5.55 98.0 34.4 5.19 1.39 1.34 9.42 121 -0.57139 11.66 5.71 97.7 33.3 4.77 0.92 0.42 13.61 117 -0.70442 11.03 5.97 97.4 33.3 5.08 1.42 0.49 10.37 113 -0.81603 11.08 6.44 97.3 33.3 5.08 1.19 0.25 10.77 109 -0.90286 11.80 7.40 97.3 33.9 6.32 2.79 1.74 7.01 106 -0.96225 14.45 10.28 97.3 35.5 7.98 4.19 2.84 5.57 103 -0.98987 44.06 41.38 98.2 58.7 41.39 39.1 38.34 34.18 101Figure 14
The velocity inx direction and its relative error.Figure 15
The velocity inz direction and its relative error.


#### 4.3.3. The Influence of Grid Density on Numerical Results

Figures16 and 17 and Tables 5 and 6 show that we divide into 9 parts in the vertical direction (Nv=9), when the circumference is divided into 20, 40, 60, 80, and 100 parts of the sphere (Nc= 20, 40, 60, 80, and 100). The velocity of the sphere in thexozplane and its relative error can be seen from the figures with the grid encryption for the velocity inxandz direction. The error is relatively small. Although the overall error of the numerical results reduces relatively after the mesh density is gradually encrypted, the accuracy has a small increase, but the relative error at the symmetry surface of the sphere (the normal vector mutation) is gradually increased with the mesh encryption, thus the infinite encrypted mesh density cannot improve the accuracy where the normal vector of meshes has severe mutation. The numerical accuracy of the present method shows that the precision of the numerical results cannot be improved obviously by mesh encryption and it is not conducive to the application of practical engineering. In this section, it is found that when the circumferential mesh is divided into 40 parts, the accuracy of numerical results can basically meet the demand.Table 5
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition N c=20 N c=40 N c=60 N c=80 N c=100 -0.08682 226.37 5.35 12.27 14.19 17.15 19.77 22.03 -0.25783 25.02 5.37 6.02 4.42 3.61 3.61 3.21 -0.42101 13.89 5.41 6.32 5.11 4.96 4.96 4.96 -0.57139 11.43 5.47 5.47 4.49 4.33 4.25 4.17 -0.70442 10.63 5.56 5.43 4.67 4.51 4.51 4.46 -0.81603 10.33 5.68 4.65 4.29 4.21 4.21 4.2 -0.90286 10.2 5.8 4.87 4.65 4.58 4.55 4.51 -0.96225 9.97 5.98 3.02 3.17 3.23 3.25 3.31 -0.98987 8.58 -- 4.08 4.19 4.22 4.22 4.22Table 6
The relative error of velocity inz direction (%).
Z/m without symmetric condition with symmetric condition N c=20 N c=40 N c=60 N c=80 N c=100 -0.08682 226.28 5.42 11.42 14.2 17.28 20.06 22.22 -0.25783 25.09 5.46 5.25 4.39 4.07 3.85 3.64 -0.42101 14.03 5.55 5.52 5.17 5.24 5.31 5.38 -0.57139 11.66 5.71 4.89 4.72 4.72 4.72 4.72 -0.70442 11.03 5.97 4.96 5.07 5.12 5.12 5.12 -0.81603 11.08 6.44 4.579 5.09 5.2 5.26 5.31 -0.90286 11.8 7.4 5.7 6.32 6.46 6.46 6.46 -0.96225 14.45 10.28 7.03 7.94 8.15 8.25 8.35 -0.98987 44.06 41.38 40.8 41.37 41.56 41.56 41.56Figure 16
The velocity inx direction and its relative error.Figure 17
The velocity inz direction and its relative error.


## 4.1. IBEM to Simulate the Uniform Flow over Sphere

Compared to the main dimension of the sphere, the computation domain is large enough, as shown in Figure4. The flow velocity of the left boundary surface (as presented in Figure 2) of the domain is set to beU (U=0.025m/s), and the velocity of the right boundary surface is alsoU. Based on the problem of uniform flow over the/a sphere (r=1.0m), the Green function of Model (a) satisfiesG=1/r1 and Model (b) satisfiesG=1/r1+1/r2.Figure 4
Computation domain and mesh of Model (a) (left) and Model (b) (right).We need to find the accuracy of the velocity distribution on the smooth/nonsmooth boundaries. The IBEM here has been used to solve the flow around the sphere whether the normal vector of the mesh has a mutation at the symmetry surface by using the symmetric condition in Green’s function. It means that the model has a mutation of the normal vector on the hemisphere without the use of the symmetric condition in Green’s function on the symmetry surface.Figures5–8 show the relative error with the analytical solution at different positions between them. As can be seen from the figure, the result of symmetric condition used is obviously closer to the analytical solution. Velocity distribution and the relative error of velocity vary little with the water depthz, and the velocity inx andz direction obtained without symmetric condition is very pronounced along with the water depth, especially for the velocity near the (x,y) plane, where the normal vector has tremendous mutation. It shows that, using the distributed source method, the accuracy of the tangential velocity obtained atz≈0 is poor and has large error. Figure 5 shows the relationship between the velocity potential of the sphere on thexoz plane with the change of coordinatez (x>0); Figure 6 shows the relation between the velocity inx direction of the sphere on thexoz plane with the change of coordinatez(x>0); Figure 7 shows the relation between the velocity inz direction of the sphere on thexozplane and the change of the coordinatez(x>0). Figure 8 shows the relative error of the velocity inxandzon thexoz plane regarding the mesh density. The curve trend of the figures indicates the increase in solution accuracy of the velocity on the sphere through increasing the mesh number, but this will reduce the computational efficiency of the IBEM.Figure 5
The velocity potential and its relative error.Figure 6
The velocity inx direction and its relative error.Figure 7
The velocity inz direction and its relative error.Figure 8
The relative error of velocity inx andz directions regarding the mesh density.


## 4.2. IBEM with Galerkin Technology to Simulate the Uniform Flow over Sphere

If the symmetric condition in Green’s function has been used and the IBEM method cannot obtain high accuracy velocity distribution on the surface of the sphere where the boundary surface is not smooth, the distribution source method based on Galerkin technology [29] has been presented to study the distribution source method. The computational domain and mesh are the same as in Section 4.1. The distributed source method based on Galerkin technology shows(33)ϕi=1ai∑j=1N∬Sidsi∬SjGσjdsj(34)2πσi+1ai∑j=1N∬Sidsi∬Sj∂G∂niσjdsj=VniNote that the Green’s functionG is similar to the above section. Since we will use quadrilateral mesh element same as before, the solution of (33) and (34) will be solved by the four point Gauss integral method. First, the weights of four Gauss points and their Jacobi determinants for each surface element are obtained, respectively, and then the influence coefficients of four Gauss points and the source points on each unit are multiplied, respectively. The weight coefficient and final summation can get the influence coefficient between the panels.Figures9–11 show whether it is based on the Galerkin technology. Whether the sphere uses the symmetric condition, the velocity and the relative error curve between the analytical solutions and the different positions at thexoz plane (x>0) have also been presented. It can be seen from the figures that the distribution source method based on Galerkin technology can obviously improve the precision of the velocity potential in most positions. The relative error of the numerical results using Galerkin technology is smaller than that of the result only using symmetric condition, but the velocity results after the use of symmetric condition are obviously closer to the analytical solution. No matter whether the Galerkin technology is used or not, the velocity inx andz direction is very intense withz when the symmetric condition has not been considered, especially for the positions where there has been normal vector mutation of the meshes, but the Galerkin technology can also improve the accuracy of velocity on the surface of the sphere, where there has not been normal vector mutation. Figure 9 shows the relationship between the velocity potential of the sphere on thexoz plane with the change of coordinatez  (x>0); Figure 10 shows the relation between the velocity inx direction of the sphere on thexoz plane with the change of coordinatez  (x>0); Figure 11 shows the relation between the velocity inz direction of the sphere on thexoz plane with the change of the coordinatez  (x>0).Figure 9
The velocity potential and its relative error.Figure 10
The velocity inx direction and its relative error.Figure 11
The velocity inz direction and its relative error.


## 4.3. DBIEM to Simulate the Uniform Flow over Sphere

In this section, DBIEM is further used to study the flow over the sphere and to improve the accuracy of the velocity distribution of the flow. Here, the symmetric condition in Green’s function has not been considered.


### 4.3.1. The Influence of the Parameterβ on the Numerical Results

Figures12 and 13 and Tables 1 and 2 give the relation between the velocity inx,z direction in thexoz plane when β=0.25, 0.5, 1.0 and 1.5, and the relative error. It can be seen from the figures that the relative error of the velocity inxandz direction at β=0.25 and 0.5 changes little with the water depthz, and the relative error increases gradually when β=1.0 and 1.5, especially at the position near the (x,y) plane (i.e., the normal vector mutation). From the above study, it is found that β=0.25 and 0.5 can get better results for the problem of flow around the sphere, but in order to ensure that the desingularized point does not coincide with the corresponding control point of the mesh, this work proposes to take β=0.5.Table 1
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition β=0.25 β=0.5 β=1.0 β=1.5 -0.08682 226.37 5.35 12.6 14.19 124 225 -0.25783 25.02 5.37 2.62 4.33 11.3 57.3 -0.42101 13.89 5.41 3.82 5.05 12.5 38.5 -0.57139 11.43 5.47 2.9 4.52 10.6 32.6 -0.70442 10.63 5.56 3.48 4.65 10.1 30.3 -0.81603 10.33 5.68 2.59 4.28 9.71 29.7 -0.90286 10.2 5.8 3.55 4.63 9.64 30 -0.96225 9.97 5.98 1.5 3.17 9.83 31 -0.98987 8.58 -- 2.91 4.21 9.98 31Table 2
The relative error of velocity inz direction (%).
Z/m without symmetric condition with symmetric condition β=0.25 β=0.5 β=1.0 β=1.5 -0.08682 226.28 5.42 12.6 14.26 124 225 -0.25783 25.09 5.46 2.72 4.43 11.4 57.3 -0.42101 14.03 5.55 3.97 5.19 12.7 38.6 -0.57139 11.66 5.71 3.15 4.77 10.8 32.8 -0.70442 11.03 5.97 3.92 5.08 10.5 30.6 -0.81603 11.08 6.44 3.41 5.08 10.5 30.3 -0.90286 11.8 7.4 5.26 6.32 11.2 31.3 -0.96225 14.45 10.28 6.4 7.98 14.3 34.4 -0.98987 44.06 41.38 40.6 41.39 44.9 57.8Figure 12
The velocity inx direction and its relative error.Figure 13
The velocity inz direction and its relative error.


### 4.3.2. The Influence of the Parameterld on the Numerical Results

Figures14 and 15 and Tables 3 and 4 give the variation of the velocity inx,z direction atld=0.1, 0.5, 1.0, 1.5, 2.0, 2.5 and 3, respectively, and the relative error. It can be seen from the figures that the relative error of the velocity inx andz direction is huge when the integration point and control point coincide withld =0.1 (ld close to zero). Asld =0.5-2.0 gradually increases, the relative error decreases rapidly, but the relative error at the symmetry surface of sphere has a sudden change, with the increase to 2.0-3.0, due to the singular distance. Its relative error increases rapidly. It is found that the relative error is smaller whenld =1.5, and the change near the symmetry surface of sphere is not obvious. Therefore, this work proposes to selectld =1.5 for correlation calculation and analysis.Table 3
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition l d=0.1 l d=0.5 l d=1.0 l d=1.5 l d=2.0 l d=2.5 l d=3.0 -0.08682 226.37 5.35 104 86.7 14.19 5.57 42.62 131.08 110 -0.25783 25.02 5.37 98.9 32.1 4.33 0.58 0.71 19.99 124 -0.42101 13.89 5.41 98 34.3 5.05 1.24 1.19 10.12 121 -0.57139 11.43 5.47 97.7 33.1 4.52 0.66 0.68 13.63 117 -0.70442 10.63 5.56 97.4 33 4.65 0.98 0.05 10.95 113 -0.81603 10.33 5.68 97.3 32.7 4.28 0.35 1.1 11.91 109 -0.90286 10.2 5.8 97.2 32.7 4.63 1.03 0.04 9.20 106 -0.96225 9.97 5.98 97.1 32.2 3.17 0.83 2.24 11.09 103 -0.98987 8.58 -- 97.1 32.5 4.21 0.48 0.76 8.7 102Table 4
The relative error of velocity in z direction (%).
Z/m without symmetric condition with symmetric condition l d=0.1 l d=0.5 l d=1.0 l d=1.5 l d=2.0 l d=2.5 l d=3.0 -0.08682 226.28 5.42 104.0 86.7 14.26 5.64 42.67 144.78 110 -0.25783 25.09 5.46 98.9 32.2 4.43 0.68 0.61 18.30 124 -0.42101 14.03 5.55 98.0 34.4 5.19 1.39 1.34 9.42 121 -0.57139 11.66 5.71 97.7 33.3 4.77 0.92 0.42 13.61 117 -0.70442 11.03 5.97 97.4 33.3 5.08 1.42 0.49 10.37 113 -0.81603 11.08 6.44 97.3 33.3 5.08 1.19 0.25 10.77 109 -0.90286 11.80 7.40 97.3 33.9 6.32 2.79 1.74 7.01 106 -0.96225 14.45 10.28 97.3 35.5 7.98 4.19 2.84 5.57 103 -0.98987 44.06 41.38 98.2 58.7 41.39 39.1 38.34 34.18 101Figure 14
The velocity inx direction and its relative error.Figure 15
The velocity inz direction and its relative error.


### 4.3.3. The Influence of Grid Density on Numerical Results

Figures16 and 17 and Tables 5 and 6 show that we divide into 9 parts in the vertical direction (Nv=9), when the circumference is divided into 20, 40, 60, 80, and 100 parts of the sphere (Nc= 20, 40, 60, 80, and 100). The velocity of the sphere in thexozplane and its relative error can be seen from the figures with the grid encryption for the velocity inxandz direction. The error is relatively small. Although the overall error of the numerical results reduces relatively after the mesh density is gradually encrypted, the accuracy has a small increase, but the relative error at the symmetry surface of the sphere (the normal vector mutation) is gradually increased with the mesh encryption, thus the infinite encrypted mesh density cannot improve the accuracy where the normal vector of meshes has severe mutation. The numerical accuracy of the present method shows that the precision of the numerical results cannot be improved obviously by mesh encryption and it is not conducive to the application of practical engineering. In this section, it is found that when the circumferential mesh is divided into 40 parts, the accuracy of numerical results can basically meet the demand.Table 5
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition N c=20 N c=40 N c=60 N c=80 N c=100 -0.08682 226.37 5.35 12.27 14.19 17.15 19.77 22.03 -0.25783 25.02 5.37 6.02 4.42 3.61 3.61 3.21 -0.42101 13.89 5.41 6.32 5.11 4.96 4.96 4.96 -0.57139 11.43 5.47 5.47 4.49 4.33 4.25 4.17 -0.70442 10.63 5.56 5.43 4.67 4.51 4.51 4.46 -0.81603 10.33 5.68 4.65 4.29 4.21 4.21 4.2 -0.90286 10.2 5.8 4.87 4.65 4.58 4.55 4.51 -0.96225 9.97 5.98 3.02 3.17 3.23 3.25 3.31 -0.98987 8.58 -- 4.08 4.19 4.22 4.22 4.22Table 6
The relative error of velocity inz direction (%).
Z/m without symmetric condition with symmetric condition N c=20 N c=40 N c=60 N c=80 N c=100 -0.08682 226.28 5.42 11.42 14.2 17.28 20.06 22.22 -0.25783 25.09 5.46 5.25 4.39 4.07 3.85 3.64 -0.42101 14.03 5.55 5.52 5.17 5.24 5.31 5.38 -0.57139 11.66 5.71 4.89 4.72 4.72 4.72 4.72 -0.70442 11.03 5.97 4.96 5.07 5.12 5.12 5.12 -0.81603 11.08 6.44 4.579 5.09 5.2 5.26 5.31 -0.90286 11.8 7.4 5.7 6.32 6.46 6.46 6.46 -0.96225 14.45 10.28 7.03 7.94 8.15 8.25 8.35 -0.98987 44.06 41.38 40.8 41.37 41.56 41.56 41.56Figure 16
The velocity inx direction and its relative error.Figure 17
The velocity inz direction and its relative error.


## 4.3.1. The Influence of the Parameterβ on the Numerical Results

Figures12 and 13 and Tables 1 and 2 give the relation between the velocity inx,z direction in thexoz plane when β=0.25, 0.5, 1.0 and 1.5, and the relative error. It can be seen from the figures that the relative error of the velocity inxandz direction at β=0.25 and 0.5 changes little with the water depthz, and the relative error increases gradually when β=1.0 and 1.5, especially at the position near the (x,y) plane (i.e., the normal vector mutation). From the above study, it is found that β=0.25 and 0.5 can get better results for the problem of flow around the sphere, but in order to ensure that the desingularized point does not coincide with the corresponding control point of the mesh, this work proposes to take β=0.5.Table 1
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition β=0.25 β=0.5 β=1.0 β=1.5 -0.08682 226.37 5.35 12.6 14.19 124 225 -0.25783 25.02 5.37 2.62 4.33 11.3 57.3 -0.42101 13.89 5.41 3.82 5.05 12.5 38.5 -0.57139 11.43 5.47 2.9 4.52 10.6 32.6 -0.70442 10.63 5.56 3.48 4.65 10.1 30.3 -0.81603 10.33 5.68 2.59 4.28 9.71 29.7 -0.90286 10.2 5.8 3.55 4.63 9.64 30 -0.96225 9.97 5.98 1.5 3.17 9.83 31 -0.98987 8.58 -- 2.91 4.21 9.98 31Table 2
The relative error of velocity inz direction (%).
Z/m without symmetric condition with symmetric condition β=0.25 β=0.5 β=1.0 β=1.5 -0.08682 226.28 5.42 12.6 14.26 124 225 -0.25783 25.09 5.46 2.72 4.43 11.4 57.3 -0.42101 14.03 5.55 3.97 5.19 12.7 38.6 -0.57139 11.66 5.71 3.15 4.77 10.8 32.8 -0.70442 11.03 5.97 3.92 5.08 10.5 30.6 -0.81603 11.08 6.44 3.41 5.08 10.5 30.3 -0.90286 11.8 7.4 5.26 6.32 11.2 31.3 -0.96225 14.45 10.28 6.4 7.98 14.3 34.4 -0.98987 44.06 41.38 40.6 41.39 44.9 57.8Figure 12
The velocity inx direction and its relative error.Figure 13
The velocity inz direction and its relative error.


## 4.3.2. The Influence of the Parameterld on the Numerical Results

Figures14 and 15 and Tables 3 and 4 give the variation of the velocity inx,z direction atld=0.1, 0.5, 1.0, 1.5, 2.0, 2.5 and 3, respectively, and the relative error. It can be seen from the figures that the relative error of the velocity inx andz direction is huge when the integration point and control point coincide withld =0.1 (ld close to zero). Asld =0.5-2.0 gradually increases, the relative error decreases rapidly, but the relative error at the symmetry surface of sphere has a sudden change, with the increase to 2.0-3.0, due to the singular distance. Its relative error increases rapidly. It is found that the relative error is smaller whenld =1.5, and the change near the symmetry surface of sphere is not obvious. Therefore, this work proposes to selectld =1.5 for correlation calculation and analysis.Table 3
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition l d=0.1 l d=0.5 l d=1.0 l d=1.5 l d=2.0 l d=2.5 l d=3.0 -0.08682 226.37 5.35 104 86.7 14.19 5.57 42.62 131.08 110 -0.25783 25.02 5.37 98.9 32.1 4.33 0.58 0.71 19.99 124 -0.42101 13.89 5.41 98 34.3 5.05 1.24 1.19 10.12 121 -0.57139 11.43 5.47 97.7 33.1 4.52 0.66 0.68 13.63 117 -0.70442 10.63 5.56 97.4 33 4.65 0.98 0.05 10.95 113 -0.81603 10.33 5.68 97.3 32.7 4.28 0.35 1.1 11.91 109 -0.90286 10.2 5.8 97.2 32.7 4.63 1.03 0.04 9.20 106 -0.96225 9.97 5.98 97.1 32.2 3.17 0.83 2.24 11.09 103 -0.98987 8.58 -- 97.1 32.5 4.21 0.48 0.76 8.7 102Table 4
The relative error of velocity in z direction (%).
Z/m without symmetric condition with symmetric condition l d=0.1 l d=0.5 l d=1.0 l d=1.5 l d=2.0 l d=2.5 l d=3.0 -0.08682 226.28 5.42 104.0 86.7 14.26 5.64 42.67 144.78 110 -0.25783 25.09 5.46 98.9 32.2 4.43 0.68 0.61 18.30 124 -0.42101 14.03 5.55 98.0 34.4 5.19 1.39 1.34 9.42 121 -0.57139 11.66 5.71 97.7 33.3 4.77 0.92 0.42 13.61 117 -0.70442 11.03 5.97 97.4 33.3 5.08 1.42 0.49 10.37 113 -0.81603 11.08 6.44 97.3 33.3 5.08 1.19 0.25 10.77 109 -0.90286 11.80 7.40 97.3 33.9 6.32 2.79 1.74 7.01 106 -0.96225 14.45 10.28 97.3 35.5 7.98 4.19 2.84 5.57 103 -0.98987 44.06 41.38 98.2 58.7 41.39 39.1 38.34 34.18 101Figure 14
The velocity inx direction and its relative error.Figure 15
The velocity inz direction and its relative error.


## 4.3.3. The Influence of Grid Density on Numerical Results

Figures16 and 17 and Tables 5 and 6 show that we divide into 9 parts in the vertical direction (Nv=9), when the circumference is divided into 20, 40, 60, 80, and 100 parts of the sphere (Nc= 20, 40, 60, 80, and 100). The velocity of the sphere in thexozplane and its relative error can be seen from the figures with the grid encryption for the velocity inxandz direction. The error is relatively small. Although the overall error of the numerical results reduces relatively after the mesh density is gradually encrypted, the accuracy has a small increase, but the relative error at the symmetry surface of the sphere (the normal vector mutation) is gradually increased with the mesh encryption, thus the infinite encrypted mesh density cannot improve the accuracy where the normal vector of meshes has severe mutation. The numerical accuracy of the present method shows that the precision of the numerical results cannot be improved obviously by mesh encryption and it is not conducive to the application of practical engineering. In this section, it is found that when the circumferential mesh is divided into 40 parts, the accuracy of numerical results can basically meet the demand.Table 5
The relative error of velocity inx direction (%).
Z/m without symmetric condition with symmetric condition N c=20 N c=40 N c=60 N c=80 N c=100 -0.08682 226.37 5.35 12.27 14.19 17.15 19.77 22.03 -0.25783 25.02 5.37 6.02 4.42 3.61 3.61 3.21 -0.42101 13.89 5.41 6.32 5.11 4.96 4.96 4.96 -0.57139 11.43 5.47 5.47 4.49 4.33 4.25 4.17 -0.70442 10.63 5.56 5.43 4.67 4.51 4.51 4.46 -0.81603 10.33 5.68 4.65 4.29 4.21 4.21 4.2 -0.90286 10.2 5.8 4.87 4.65 4.58 4.55 4.51 -0.96225 9.97 5.98 3.02 3.17 3.23 3.25 3.31 -0.98987 8.58 -- 4.08 4.19 4.22 4.22 4.22Table 6
The relative error of velocity inz direction (%).
Z/m without symmetric condition with symmetric condition N c=20 N c=40 N c=60 N c=80 N c=100 -0.08682 226.28 5.42 11.42 14.2 17.28 20.06 22.22 -0.25783 25.09 5.46 5.25 4.39 4.07 3.85 3.64 -0.42101 14.03 5.55 5.52 5.17 5.24 5.31 5.38 -0.57139 11.66 5.71 4.89 4.72 4.72 4.72 4.72 -0.70442 11.03 5.97 4.96 5.07 5.12 5.12 5.12 -0.81603 11.08 6.44 4.579 5.09 5.2 5.26 5.31 -0.90286 11.8 7.4 5.7 6.32 6.46 6.46 6.46 -0.96225 14.45 10.28 7.03 7.94 8.15 8.25 8.35 -0.98987 44.06 41.38 40.8 41.37 41.56 41.56 41.56Figure 16
The velocity inx direction and its relative error.Figure 17
The velocity inz direction and its relative error.


## 5. Conclusion

The problem of the flow on the smooth and nonsmooth boundaries has been modelled through Indirect Boundary Element Method (IBEM) and Desingularized Boundary Integral Equation Method (DBIEM). The symmetric condition of Green’s function on the symmetric plane surface of the sphere has been considered. Whether it has normal vector mutation at the junction of hemispherical and symmetry surface through this way, which means the boundary at the fictitious junction is smooth by using symmetric condition, has been considered, too. Then numerical simulation and convergent analysis have been carried out for the mentioned model, from which the following conclusions can be drawn.( 1 ) The result of IBEM has been obviously improved by using symmetric condition of Green’s function. The velocity in bothx andz direction obtained without using symmetric condition is very intense along with the water depth, especially in place where the normal vector has severe mutation and the boundary is nonsmooth.( 2 ) The distribution source method based on Galerkin technology can improve the precision of the velocity potential, and, in most positions, the relative error of the numerical results using Galerkin technology is smaller than that of the result only using symmetric condition, but the velocity precision after the use of symmetric condition is closer to the analytical solution no matter whether the Galerkin technology is considered or not.( 3 ) The accuracy of velocity has been greatly improved by DBIEM when β=0.5 andld =1.5 have been undertaken for the flow simulations, where the symmetric condition of Green’s function has not been involved in the simulation. Although the overall error of the numerical results is relatively reduced after the mesh density is gradually encrypted, the accuracy has a small improvement, but the relative error near the symmetric plane surface of the sphere (the normal vector has severe mutation) is gradually increased with the mesh encryption. This will affect the efficiency of the numerical calculation and is not beneficial for the application in practical engineering.( 4 ) The present work has not considered the flow on the complex structures, which forms part of our further work.



---

*Source: 4131683-2019-07-07.xml*