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.
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 Section 2. 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 velocity v=▽ϕ, which is governed by the Laplace equation in the fluid domain D as shown in Figure 1; the equation for the velocity potential ϕ becomesBased on Gauss formula and three-dimensional Green’s third formula, we can find the value ϕ(p) by DBEM in the flow domain D and ϕ(p) on the surface S=S1+ S2+ S3+ S4+ S5+ S6 of the flow domain D byHere, p(x,y,z) is a field point and q(ξ,η,ζ) is a source point on the surface S of the fluid domain D. The geometrical coefficient C(p) =2π when the field point is located on the surface S; otherwise, C(p) =4π. G(p,q) is a Green’s function, where r is the distance between field point p and source point q.
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 velocity U, z points vertically upwards from the symmetry surface of the sphere as shown in Figure 2. The sphere surface is denoted by Ss and its unit normal vector directed outward from the fluid region is denoted by . The boundary surface S of the Model (a) and Model (b) is S= Ss + Sa-b + + + Se-f + .
(a) G without symmetric condition
(b) G with symmetric condition
Based 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 surface Sa-b-c-d:Using the Direct Boundary Element Method (DBEM) in Section 2.1, we haveand the source strength σThen we can use Indirect Boundary Element Method (IBEM) to find the velocity distribution of the fluid domain byHere, the Green’s function G for the Model (a), which does not consider the symmetry surface, satisfiesFrom (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 surface Sa-b-c-d. This means the boundary at junction is nonsmooth. Model (b) however uses the symmetric condition of Green’s function on the symmetry surface Sa-b-c-d. Thus the mesh’s discretization on the symmetry surfaces Sa-b and 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 satisfieswhereHere, p (x, y, z) is the field point in the considered domain, and the source point q’(ξ, η,-ζ) is the mirror point of the source point q(ξ, η, ζ) on the discretized surface S using the xoy 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 surface S outside the fluid domain D. 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 domain D can be written as follows:where q(ξ,η,ζ) is the integration point on the integration surface S 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 surface S.
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
where xs=(ξ,η,ζ): a source point on the integration surface, xc=(x,y,z): a field point on the real boundary, ϕ0=the given velocity potential value at xc, Гd=surface on which ϕ0 is given, χ=the given normal velocity of velocity potential at xc, =surface on which χ is given, G(p, q)=1/r(p, q), where r is the distance between field point p and the source point q.
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 of N-isolated singularities located at a small distance away from the corresponding control point on the boundaries,The desingularized distance between the/an isolated source point and corresponding control point is given bywhere ld and β are constants and Dm 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 of ld and β. Therefore, appropriate ld and β values need to be determined after numerical investigation. The recommended values are ld =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,
3. Numerical Implementation
3.1. Numerical Implementation for DBEM
For the DBEM, (3) is solved using/applying the following numerical procedure, in which the surface S is discretized into finite numbers of facets. The corresponding values of 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:In which nS are the numbers of facets on domain surface, respectively.
The matrix coefficients and correspond to integrals of the Green’s function and its normal derivative over the area ∆Q of the jth facet, respectively. and are written as
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)
Model (b)Then the velocity distribution for the x, y, and z direction in the flow field can be obtained by
Model (a)
Model (b)In (25)-(28), nS are the numbers of facets on domain surface. The matrix coefficients Dij correspond to the normal derivative of Green’s function over the area ∆Q of the jth face, as presented in (24) and the matrix coefficients DDij are presented in (29).
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:and the velocity can be obtained byThe matrix coefficients NDij are presented in (32).
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 Figure 4. The flow velocity of the left boundary surface (as presented in Figure 2) of the domain is set to be U (U=0.025m/s), and the velocity of the right boundary surface is also U. Based on the problem of uniform flow over the/a sphere (r=1.0m), the Green function of Model (a) satisfies G=1/r1 and Model (b) satisfies G=1/r1+1/r2.
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.
Figures 5–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 depth z, and the velocity in x and z 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 at z≈0 is poor and has large error. Figure 5 shows the relationship between the velocity potential of the sphere on the xoz plane with the change of coordinate z (x>0); Figure 6 shows the relation between the velocity in x direction of the sphere on the xoz plane with the change of coordinate z (x>0); Figure 7 shows the relation between the velocity in z direction of the sphere on the xoz plane and the change of the coordinate z (x>0). Figure 8 shows the relative error of the velocity in x and z on the xoz 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.
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 showsNote that the Green’s function G 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.
Figures 9–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 the xoz 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 in x and z direction is very intense with z 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 the xoz plane with the change of coordinate z (x>0); Figure 10 shows the relation between the velocity in x direction of the sphere on the xoz plane with the change of coordinate z (x>0); Figure 11 shows the relation between the velocity in z direction of the sphere on the xoz plane with the change of the coordinate z (x>0).
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
Figures 12 and 13 and Tables 1 and 2 give the relation between the velocity in x, z direction in the xoz 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 in x and z direction at β=0.25 and 0.5 changes little with the water depth z, 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.
4.3.2. The Influence of the Parameter ld on the Numerical Results
Figures 14 and 15 and Tables 3 and 4 give the variation of the velocity in x, z direction at ld=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 in x and z direction is huge when the integration point and control point coincide with ld =0.1 (ld close to zero). As ld =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 when ld =1.5, and the change near the symmetry surface of sphere is not obvious. Therefore, this work proposes to select ld =1.5 for correlation calculation and analysis.
4.3.3. The Influence of Grid Density on Numerical Results
Figures 16 and 17 and Tables 5 and 6 show that we divide into 9 parts in the vertical direction (=9), when the circumference is divided into 20, 40, 60, 80, and 100 parts of the sphere (= 20, 40, 60, 80, and 100). The velocity of the sphere in the xoz plane and its relative error can be seen from the figures with the grid encryption for the velocity in x and z 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.
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.
The result of IBEM has been obviously improved by using symmetric condition of Green’s function. The velocity in both x and z 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.
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.
The accuracy of velocity has been greatly improved by DBIEM when β=0.5 and ld =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.
The present work has not considered the flow on the complex structures, which forms part of our further work.
Data Availability
The data used to support the findings of this study are included within the article.
Disclosure
Permanent address: School of Naval Architecture and Ocean Engineering, Jiangsu University of Science and Technology, No.2, Mengxi Road, Zhenjiang, China.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (Grants No. 51709137 and 51579120), National Natural Science Foundation of Jiangsu (BK20151327), Key University Science Research Project of Jiangsu (18KJA130001), Marine Equipment and Technology Institute of Jiangsu University of Science and Technology Funding Project (1174871801-17), and Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX18_2343).