Unsteady 2D and 3D Navier-Stokes Solver with Application of Multigrid Scheme to Pressure Poisson Fractional Step on Arbitrary Unstructured Grids in Various Applications with Emphasis on Ship Motion

,


Introduction
Iterative methods have many applications in engineering problems especially in Computational Fluid Dynamics (CFD). It is found that, in CFD problems, the accuracy is directly related to the number of grids in computational domains. Unfortunately, increasing the mesh number is not always possible, because it directly increases the computational time in common iterative methods. In recent decays, many extensive researches have been conducted to accelerate the speed of common iterative methods. One of them is the multigrid method which is widely used in CFD codes.
Common single-grid iterative solvers, like SOR and its family, have poor convergence characteristics, especially in finer grids. ey are told to produce a smooth error. In other words, they remove high frequencies error easily, but low frequency reduces very slowly. e main idea of the multigrid method is to change the solution to a coarse grid, on which smooth error is rough and low frequency acts like higher ones.
Many researchers used the multigrid method to solve linear systems of equations in order to speed up the solvers. In 1995, Zang and Street [1] applied a 3D finite volume multigrid solver, based on a second-order fractional step scheme, to simulate Navier-Stokes equations. For data transferring, they used biquadratic interpolation between subdomains. In 1996, Meza and Tuminaro [2] used the multigrid method as a preconditioner combined with a conjugated method to solve drift-diffusion equations. Lin [3] in 1994 developed a multigrid method for computing the unsteady inviscid flows on unstructured grids. On the other hand, Liu et al. [4] used a dual time stepping multigrid scheme on structured grids to simulate incompressible flows. Ginzburg and Wittum [5] in 2001 used the V cycle multigrid method as a preconditioner for the biconjugate stabilized method to solve Navier-Stoke equations combined with the VOF method. ey used SIMPLE and Gauss-Seidel transforming operators. Subsequently, Tang et al. [6] in 2003 used the same way and applied a V cycle multigrid method to their 3D unsteady incompressible overset-grid solver. rough this method, the data transferring between subdomains takes place on the finest grids. In 2003, ai and Zhao [7] presented a parallel Navier-Stokes solver based on an unstructured finite volume method, using a multigrid domain decomposition approach. A year later, they [8] successfully applied a high-order finite volume scheme to simulate unsteady incompressible flows, by introducing an unstructured multigrid method. Lambropoulos et al. [9] presented a parallelization of the agglomeration multigrid technique to solve Favre-average Navier-Stokes equations on unstructured grids. e discretization method was finite volume method and FAS (full approximation scheme) was implemented on 2D/3D cases such as airfoil and wing to demonstrate the efficiency of the model. In 2014, Li et al. [10] tried to show that the multigrid method is more effective when the residual restriction satisfies the flux conservation principle. ey conducted a study on four convective-diffusion cases and found that the multigrid acceleration is most desirable when the residual restriction satisfies the flux conservation principle. Later on, Lv et al. [11] in 2006 implemented this approach to simulate 3D unsteady compressible flows with arbitrary moving objects. ey applied a unique combination of a parallel multigrid method with low Mach number preconditioning. In 2007, Ben Cheikh et al. [12] investigated the efficiency and accuracy of the multigrid method to solve Navier-Stokes/Boussinesq equations. ey applied the AFMG (Accelerated Full Multigrid Method) for the solution of an 8 : 1 differentially heated cavity. ey showed that AFMG could speed up the CPU time up to 251 times, related to RBSOR (Red and Black Overrelaxation). Griffith [13] presented a 2D unsplit, staggered discretization of incompressible Navier-Stokes equations with an efficient solver for a system of linear equations. He used the Krylov subspace method with the projection method preconditioner. e preconditioner employed CG (conjugated gradient) method with a geometric multigrid method. He validated his model with several problems such as the lid-driven cavity at various Reynolds numbers. In 2010, Santhosh et al. [14] presented a three-dimensional, transient, and incompressible Navier-Stokes solver. e discretization method is a fractional step method which is second-order time accuracy and the multigrid method assists its Poisson solver. To show the accuracy and efficiency of the code, 3D lid-driven cavity was simulated by the 4-level V cycle multigrid method. Zong-zhe et al. [15] in 2011 developed a new agglomeration multigrid method to generate coarse grids, automatically. ey implement their algorithm to solve compressible Navier-Stokes equations and test that on NACA0012 at Reynolds number Re � 5000. In 2013, San and Staples [16] presented a coarse grid projection method to accelerate incompressible flow. is method was introduced to apply to problems including the Poisson equation. ey applied the V cycle of their computational method to several 2D/3D test cases in order to show its efficiency. ey found that computational time decreases by increasing the distortion ratio on non-Cartesian grids. Langer [17] in 2014 developed a finite volume model with an agglomeration multigrid method to calculate the steady solution of compressible RANS equations on unstructured grids. He used the Gauss-Seidel method to solve the linear systems of equations. Implementing the presented model on two/threedimensional examples, he showed a wide range of applicability of the solution. In 2014, Cheng and Samtaney [18] developed a high-resolution code to simulate large eddy simulation of incompressible turbulent boundary layers over a flat plate. ey used a fourth-order multigrid method to solve the modified Helmholtz equation and the Poisson equation. In 2016, Stiller [19] developed a multigrid method for two-dimensional grids in order to apply the Galerkin formulation of the Poisson equation. He proposed his model to be used in more complicated applications such as CFD. Lygidakis et al. [20] in 2016 presented an agglomeration multigrid methodology which was applied on two academic CFD codes: GALATEA and GALATEA-I. ese two codes are able to simulate compressible and incompressible flows on unstructured grids. Comparing different approaches of multigrid methodology, they found that, in general, higher acceleration would archive in incompressible flow compared with the compressible flow. at is due to the elliptical nature of incompressible equations. In 2017, Matthioudajis et al. [21] developed a fourth-order finite-difference model on cellcentered grids to calculate Navier-Stokes equations. ey introduced a new high-order transfer operator to apply the multigrid method on pressure correction. In their work, a comparison was made between different multigrid cycles. ey showed that the accuracy of the method is retained in time and space by increasing the Reynolds number. Abide and Zeghmati [22] in 2019 investigated a multigrid defect correction to solve a fourth-order compact scheme of a singular Poisson equation. ey applied the presented method to the Hodge-Helmholtz decomposition test case and showed better accuracy related to second-order discretization. In 2019, Liu et al. [23] proposed a mode multigrid method and applied it to a steady-state solver for unstructured grids. ey used the solver for 2D and 3D foils and showed that the method is 3 to 6 times faster than the baseline method while ensuring computational accuracy. Ha et al. [24] in 2019 presented an interpolation operator on a multigrid method based on distance weight. e performance of the method was compared with that of an interpolation operator by area (volume) intersection for the finite volume discretization. By solving heat conduction problems on 2D and 3D domains, they concluded that the distance weight operator is faster and more reliable than the area intersection approach.

Background and
is Paper. A three-dimensional Navier-Stokes solver is developed based on a fractional step method which was introduced by Kim and Choi [25] and Panahi et al. [26]. e finite volume method was used to discretize the equations on unstructured grids. e solver is combined with CICSAM [27] volume of fraction method to capture the free surface elevation in various applications, especially in ship motions. By applying the body-attached mesh and calculating the pressure field in the computational domain, hydrodynamic forces and moments can be calculated on the ship hull [28]. As pointed out, the common single-grid iterative solvers have poor convergence characteristics. In this paper, the agglomeration multigrid method is applied to the Poisson solver of fractional steps in order to speed up the iterative solver. For this purpose, a fully unstructured database is established for the solver. e finest grid is introduced as an input to the algorithm (it could be generated by an external mesh generated software like Gambit). e coarse grids are automatically generated through the presented agglomeration algorithm. is method was initially introduced by Smith [29] and Lallemand et al. [30]. In this paper, in order to generate coarse grids with better performance, a new coarsening method is applied and for more accuracy, a rather complicated data structure is designed to save all faces of the coarsened cell. ese nonstandard cells are used to save the geometrical details of the ship like rudder and brackets. In order to verify the accuracy of the solver and investigate its performance, various 2D and 3D test cases are simulated. e acceleration of V and W multigrid cycle is also compared in all simulations.

Governing Equations
Two main governing equations for unsteady, incompressible Newtonian fluid motion are momentum and conservation equations. ese equations in three-dimensional space are written as follows: where u is the velocity, ρ the density, P the pressure, and g i the gravitational acceleration. In the cases where there are two-phase fluid simulations, one should solve the volume fraction transport equations as in In this equation, α is a fractional value which demonstrates the amount of water in the cell: Cell inside water, 0, Cell inside air, 0 < α < 1, for transitional area.
Generally, equations (1)-(3) should be solved simultaneously. For this purpose, the fluid density and viscosity of each cell in the computational domain can be expressed as a function of the volume fraction α as follows: where the subscripts 1 and 2 denote water and air, respectively.

Discretization of Navier-Stokes Equations.
As previously pointed out, the finite volume method is applied to discretize the governing equations including Navier-Stokes equations.
Integrating Navier-Stokes equations over each control volume can lead to For the fluid velocity u i , the diffusion term (the first term on the right-hand side of equation (6)) is discretized using overrelaxed interpolation Jassak [31]: In order to discretize the convection term (the second term in the left-hand side of equation (6)), fluid velocity at the control volume faces is needed to compute fluxes. is can be evaluated as follows: To calculate u f , the gamma interpolation scheme is used which was introduced by Jassak [31].
It must be mentioned that the Crank-Nicholson scheme is used for time discretization of diffusion and convection terms in equation (6). is scheme is of second-order accuracy and increases stability. e pressure term (second term on the right-hand side of equation (6)) is discretized as where A f is the direction component of the face area vector. Using the common linear interpolation between two neighboring control volume centers results in severe oscillations in the velocity field. is is of great importance, especially when there are two fluids with high-density ratios, for example, water and air in this paper. Here, a Piecewise Linear Interpolation is used which was developed by Panahi et al. [26].

Discretization of VOF Equation.
e finite volume discretization of volume fraction transport equation (3) is based on the integration over the control volumes and the time step: e first term is a simple integral form. By applying the Gauss theorem on the second term and assuming a small variation of F f , it can be shown that where (α f ) in equation (11) must be approximated on the face of each control volume. Simple interpolation results in nonphysical values. is leads to using a high-order composite interpolation. In this paper, CICSAM (Compressive Interface Capturing Scheme for Arbitrary Meshes) scheme is applied [27]. CICSAM uses CBC (Convection Boundedness Criteria) (α f CBC ) and UQ (Ultimate Quickest) (α f UQ ) interpolation by introducing a weighting factor c f [32,33]. is factor is evaluated according to the slope of the free surface and the direction of fluid motion [27]. Based on NVD (Normalized Variable Diagram), normal face value is obtained as follows [32]: NVD definitions are shown in Figure 1 [26] where donor, acceptor, and upwind cells are defined according to flow direction for each control volume face: Substituting equation (13) into equation (12) results in the estimation of (α f ).
In the VOF approach, CFL (Courant-Friedrichs-Lewy) number plays an important role in the stability of the solution. e Courant number is calculated according to equation (14) for each cell and is set to be less than 0.2 for all two-phase flow test cases [28]: As mentioned before, a body-attached mesh is used to simulate the body motions. In this way, the relative velocity should be applied in all equations as follows: where u → fluid is fluid velocity vector and u → body is the velocity of the moving body (which is equal to mesh velocity).

Numerical Algorithm
e main algorithm which is used to couple the pressure and velocity field is the modified fractional steps which were firstly proposed by Kim and Choi [25]. Consider the Navier-Stokes equations for the velocity component u i in its discretized form on a given control volume. e first step is to solve the first intermediate velocity (u i ) with the pressure gradient in the previous time step (G i (P n )): where en, the second intermediate velocity is calculated as Using equation (18) results in the pressure Poisson equation which gives a new pressure field P n+1 : By writing equation (19) for all cells in the computational domain, a system of equations in the form (AX � B) is obtained. is matrix is solved by applying the multigrid method which will be explained in the "Multigrid Method." After that, the velocity at the new time step is computed as follows: Since the new pressure is used to calculate (in other words, correct) velocities, global mass conservation is automatically satisfied. After the calculation of all velocity components in an outer iteration, the control volume face velocity vector is calculated for the next time step from Figure 1: Flow direction determines the donor, acceptor, and upwind cells [26]. equation (22). It includes the effect of pressure gradient on the calculation of the face velocity vector to overcome the check-board pressure in the nonstaggered (collocated) arrangement: e overall algorithm is shown in Figure 2. Also, more details on the algorithm can be found in a paper by Panahi et al. [28].

Calculating Forces, Moments, and Rigid Body Motions
As pointed out earlier, the body-attached mesh (in which mesh and body move together) is used in order to simulate the motion of rigid bodies, like ship hull ( Figure 3). Solving the Navier-Stokes equations, one can calculate the flow forces and moments by integrating the pressure and viscous stress over the body boundaries (like ship hull). ese forces ( F → ) and moments ( M �→ ) cause the body to move in 6 degrees of freedom. ese movements could be calculated by solving the linear and angular momentum equations (equations (23) and (24)): Body cells where (I, α and ω) are body inertia matrix, body angular acceleration, and body angular velocity, respectively.

Multigrid Method
As previously pointed out, conventional iterative single-grid methods like SOR for solving Navier-Stokes equations on fine grids suffer from poor convergence characteristics. ese solvers remove high frequencies and short wavelengths, so they are known to have smooth property. Multigrid methods effectively remove low-frequency error components on finer grids, by converting the results to coarse grids. In this paper, the finest grid for the computational domain can be generated by any commercial software such as Gambit or ANSYS. en, coarse grids would be generated based on the coarsening algorithm.

Coarsening Algorithm.
e coarsening algorithm is based on the method which Okomoto et al. [34] suggested in 1998. e main steps of the algorithm are described as follows: (1) e starting vertex in the fine grid should be selected (e.g., vertex no. 1 in Figure 4) (2) In the second step, all the cells which contain the sharing vertex are merged (3) e next vertex is selected (e.g., vertex no. 2 in Figure 4) (4) If one of the sharing cells (sharing cells of a vertex: the cells which the vertex is one of its nodes) of this vertex (selected in step 3) is used in previous coarsening, then go to step 3 and select another vertex; else go to step 2 Figure 4 shows the above algorithm for a simple 8-cells grid. In this example, after merging all the sharing cells of vertex number 1, vertex number 2 is passed and vertex number 3 is selected for the next coarsening. Figure 4 shows a very simple cell arrangement while in reality, it is not always the case.
A slightly more complicated cell arrangement is displayed in Figure 5. In this case, after coarsening, two larger cells are produced. Figure 5(a) shows the larger cells with simplified faces. Although standard unstructured cells (Hexahedron, Wedge, Tetrahedron, and Pyramid are referred to as standard meshes) are produced in this method, eliminating the details may compromise the accuracy, especially near the important boundaries like no-slip wall (ship hull). In Figure 5(b), the coarsening is done by creating Calculate 2nd intermediate velocities u * using equation (20) Calculate pressure using poisson equation (21) Calculate the velocity in new time step using equation (23) Calculate flux at faces using equation (24) Advance in time t n+1 = t n + dt  Mathematical Problems in Engineering 5 polyhedron cells and without eliminating the details. In this method, a more complicated data structure is needed in order to deal with arbitrary polyhedron cells with an arbitrary number of faces. (In coarse grids, there is more than one face between two specific cells.) Although it is already known that this type of data structured requires more memory allocations which is discussed in the following sections. is method is more applicable in complex three-dimensional geometries like simulating ships with all appendages. For better coarsening and efficient simulation, it is recommended that the vertexes not be chosen, randomly. erefore, before the coarsening, a list of vertexes that have more sharing cells must be made, and the coarsening should follow according to the list. With this recommendation, the highest efficiency in grid coarsening is achieved. Based on the explained algorithm, three or four levels of coarse grids can be generated due to the size of the computational domain.

Multigrid Algorithm.
As mentioned before, the multigrid technique is used to damp the low-frequency error components. Based on the number of coarse grids, the multigrid method can be applied. Here, for simplicity, a twogrid scheme is introduced. e implementation of the algorithm for more grids is very similar. Assume the linear system of equation (Ax � B). e two-grid (multigrid) scheme to solve this system of equation is as follows [35]: (1) Relax with (A 1 x 1 � B 1 ) on fine grid (with common iterative methods like SOR) and find the approximate solution v 1 (relax until the convergence rate decreases) (2) Compute the fine grid residual (r 1 � f 1 − A 1 v 1 ) and restrict it to the coarse grid (3) Solve (A 2 e 2 � r 2 ) on the coarse grid (4) Prolong the coarse grid error to the fine grid and correct the fine grid approximation guess v 1 e subscripts 1 and 2 denote the fine and coarse grid, respectively. is algorithm can be extended on more grids and in different cycles. Figure 6 illustrates the structure of one iteration step (cycle) of a multigrid method with a few pictures. e cases c � 1 and c � 2 are usually used. For obvious reason, the cases c � 1 and c � 2 refer to V and W cycles, respectively ( Figure 4). e number c is also called a cycle index [36]. In this paper, V and W cycles are used for comparison in all simulations.t

Restriction and Prolongation Operators.
As mentioned in the previous section, the multigrid method requires some operators to transfer data between the grids. e prolongation operator transfers a set of variables from a coarse grid to a fine grid and the restriction operators transfer a set of variables from a fine grid to a coarse grid. Many different operators are defined for prolonging and restriction. In this paper, the operators which are introduced by Ferziger and Peric [37] are used:   Mathematical Problems in Engineering Restriction operator: Prolongation operator: In equations (25) and (26), r is the cell center of the control volume, (N f ) is the number of the fine grid control volume, and f and c denote the fine and coarse grid, respectively. Figure 7 clearly shows the operators.

Numerical Results
Based on the algorithm described above, a 3D numerical code is developed and its accuracy is verified and different simulations are performed. In each case, different multigrid cycles are applied and the computational times are compared to see which cycle provides better efficiency in each case. For better comparison, two different parameters are considered: residual and computational time. e residual of each iteration for the Poisson equation is defined as follows: 7.1. 2D Triangular Lid-Driven Cavity. Firstly, a two-dimensional test case is investigated. e lid-driven cavity is a very common test case to check the accuracy of Navier-Stokes equations (one-phase flow). In this case, the domain is a triangular shape, as shown in Figure 8. e upper side is moving with unit velocity while the other two sides are assumed as no-slip walls. Meanwhile, 28,731 triangular unstructured cells are used to discretize the computational domain. A section of the meshed domain (section z-z in Figure 8) is shown in Figure 9. e Reynolds numbers are set to be (R e � (Vl/v)). is is a steady-state test case, so the simulation should be performed until the velocity and pressure field become steady. Velocity fields are plotted along two horizontal and oblique lines. e horizontal line joins the center of AB to the center of AC and the oblique line joins point A to the center of BC. Results are compared with the computational results of Jagannathan et al. [38] in Figures 10 and 11. Coarse grid Fine grid   Table 2. In this case, VMG4 and WMG4 do not offer impressive performance in reducing residuals; it seems that their difference with VMG3 and WMG3 is very small.
Based on these two parameters, the W cycle with 4 grids offers better efficiency. is cycle decreases the simulation time twofold.

Backward-Facing
Step. In this case, 2D flow over a step in a pipeline is investigated known as BFS (Backward-Facing Step) problem. e schematic geometry of the problem is shown in Figure 13. At the inlet boundary condition, a fully developed flow is applied. e Reynolds number, in this case, is computed (Re � (VD/]) � 200). D is the pipe diameter. A fully structured grid is used as the finest grid level.
Due to the Reynolds number, two back warding flows are developed at the top and bottom of the pipe, as shown in Figure 13. ere are 30,400 structured cells used for simulations. When the fluid flow becomes steady, the velocity profile at two sections (x � 6, 14) is plotted and compared by similar works [39,40] (Figures 14 and 15).
For better comparison, it is also necessary to calculate the length of the backward flow in the pipeline (X 1 , X 2 , X 3 ).
ese lengths are compared with those in Ertrukʼs work [40] displayed in Table 3. Results indicate good agreement in simulation.
In this test case, a three-grid level V and W cycle is used for the simulations. Table 4 shows two coarsened grid size level. Since coarse grid level 2 consists of 771 cells, more coarsening is not effective in the solution. erefore, 3 grid levels are used in multigrid cycles.
Using these grid sizes, the residual reduction is shown in Figure 16.
e computational times are also compared in Table 5 for 50 time steps.
Usually, the WMG3 cycle reduces the error in fewer iterations, but in this particular case, the VMG3 operates with a higher speed. is is because the number of cells is very low and with prolonging once, error oscillation in    coarser cells gets sufficiently smooth (V cycle). More usage of coarser cells and transferring more information between the cells (in this particular case) increases the CPU time (W cycle). erefore, in this case, VMG3 is a better iterative cycle to be selected.

3D Cubical Lid-Driven Cavity.
Consider the three-dimensional cubical lid-driven cavity with sideʼs length of (l � 1 m), as depicted in Figure 17. e top wall is moving to the right with constant velocity (v � 1(m/s)), while for other faces, the stationary no-slip wall condition is applied. e cavity is filled with a Newtonian fluid (constant viscosity and density). In this case, the simulations are performed at three Reynolds numbers (Re � (vl/]) � 100, 400, and 1000). e computational grid of 50 × 50 × 50 control volumes is considered as a fine grid with (Δt � 0.01 s). e iteration loop is continued until (Log(Residual) < − 9) is achieved. After the solution becomes steady, the velocity component in X direction (u) is plotted along y-axes at the midplane      [42] at different Re numbers (Figure 18). e results are in agreement and this shows that a multigrid scheme can predict velocity distribution with good accuracy. Figures 19-21 compare the pressure contours at different sections (x � 0.5, y � 0.5, and z � 0.5) with Jiang et al. [41] at Re � 100, 400, and 1000, respectively.
As pointed out earlier, for solving the pressure Poisson equation, the multigrid method is used. e number of cells in four different grid sizes is shown in Table 6. As observed, the coarsening method can efficiently merge the cells and produce grid level 1 with less than (1/4) cells of the finest grid. e residual reduction versus iteration for the second time step is illustrated in Figure 22. Table 7 also shows the CPU time and memory allocations of different iterative methods for 50 time steps.
Based on the numerical results of Table 8, W cycle with 4 grids is used for the entire simulation of a 3D cubical liddriven cavity. It is clear that, for this case, the multigrid approach can easily increase the speed of computations by 1.85 times.
In order to better evaluate the current solver, a more complicated 3D cavity, two-sided nonfacing lid-driven (3D-TSNFL) is investigated (Figure 23). In this test case, the top wall is moving to the right, while the left vertical wall is moving down with the same constant velocity. All the other parameters are exactly the same as the 3D lid-driven cavity which was mentioned earlier.
e computational domain consists of 39 × 39 × 39 control volumes as a fine grid with (dt � 0.01 s). After the solution becomes steady, the velocity components (u and w) are plotted along zand x-axes at the midplane (y � 0.5) and compared with similar work of Ben Beya and Lili [43] ( Figure 24). e velocity contours are also illustrated in Figures 25  and 26 and compared with that provided by Ben Beya and Lili [43]. e number of cells in four different grid sizes is shown in Table 9 which were used in the multigrid method. As observed, the coarsening method can efficiently merge the cells and produce grid level 1 with less than (1/4) cells of the finest grid. e residual reduction versus iteration for the second time step is shown in Figure 27. Table 8 also shows the CPU time and memory allocations of different iterative methods for 50 time steps.
Based on the numerical result in Table 8, W cycle with 4 grids is used for the entire simulation of 3D-TSNFL. It is clear that, for this case, the multigrid approach can easily increase the speed of computations by 1.8 times.

Collapse of a Water Column with Obstacle.
is test case is a three-dimensional two-phase flow problem. e geometry and initial conditions are illustrated in Figure 28. In the first time step, the column of water breaks and flows into the domain. is problem is a very complicated flow and because of the lower density of air (related to water), it gets trapped in the water. Trapped air is subjected to a large buoyancy force and tends to rise up. e pressure which water exerts on the obstacle along with water height at two different stations of the reservoir is the main parameter which is investigated in this test case.
Based on the computational domain, a fully structured grid is used as the finest grid level, similar to BFS. According to previews cases, a multigrid method is used with a three-grid level for simulation. e grids level size is shown in Table 10.
As observed, the proposed coarsening algorithm has a better performance in fully unstructured grids. Comparison of the coarsening in an unstructured triangular cavity test case with the structured BFS and 3D dam break grids reveals that, in unstructured grids, the coarsening ratio is about 4.5, while in the structured grids, the ratio is almost 3.9 to 4.0.
is is due to the fact that, in structured grids, the number of sharing cells for each vertex is less than the unstructured grids. For example, in a triangular cavity, there are vertexes with 6 sharing cells, while this number is 4 in BFS or dam break test case. e best multigrid cycle is chosen due to the results of the residual reduction in the second time step Figure 29 and Mathematical Problems in Engineering also the simulation time of the problem for 100 time steps which is shown in Table 11. In this simulation, it is concluded that W cycle is slightly better than V cycle because of a relatively large number of cells. Hence, W cycle is selected for the entire simulation. e simulation is performed for 6.0 seconds. As shown in Figure 28, the water height is computed at points H1 and H2 and pressure is numerically measured on the obstacle at points P1 to P3. is test case is experimentally conducted at the Marine Research Institute of the Netherland (MARIN) and results were published by Kleefsman et al. in 2005 [44]. Figure 30 shows the time history of water height at H1 and H2 and is compared with experimental measurements.
In both pictures of Figure 30, the agreements are good until the water has returned from the back of the reservoir (after about 1.5 s). After that, although some differences are observed, global behavior is almost still the same. e pressure histogram is also plotted at points P1 to P3 in Figure 31. It is demonstrated that the presented solver can estimate the pressure with acceptable accuracy related to experimental data, especially at points P1 and P2. Ultimately, the free surface profile is compared to experimental data in Figure 32. It is shown that CICSAM scheme can capture the free surface shape almost the same as the actual case which makes it a reliable approach in VOF method.

Barge Resistance.
In this section, the resistance of a moving barge is calculated as a simple marine test case. Figure 33 and Table 12 show the main parameters of the barge geometry. Barge resistance was experimentally measured in the towing tank of the marine laboratory at Sharif  University of Technology [26] at a constant velocity (V � 0.807(m/s)).
Due to water forces, the draft and trim of the barge may change. erefore, the simulation is performed in two degrees of freedom motion which basically means that the barge is allowed to have pitch and heave motions. A total of 158,891 structured grids are used as the finest grid size to simulate barge motion ( Figure 34). No-slip boundary condition is used for barge and outflow boundary condition is applied for other 6 boundaries.  Using the multigrid method, the residual reduction for the second time step is illustrated in Figure 35.
e simulation is conducted with all types of cycles for 50 time steps and as a result, the W cycle with a four-grid level is adopted because of the better CPU speed-up ratio Table 14. Figure 36 shows the resistance of the barge versus time. After about five sinusoidal oscillations, the resistance  becomes steady. ese oscillations appear due to the heave and pitch motion of the barge. According to Figure 36, it takes about 5 seconds for barge motions to become stable. e computed resistance is compared with experimental data in Table 15. e result shows that the presented numerical solver can estimate barge resistance with good accuracy.
e free surface near the barge is also compared with the experimental photo in Figure 37.

Resistance and Motion of a High-Speed Catamaran.
As seen in the previous section, a low-speed barge has low heave and pitch motion. erefore, in order to examine the capability of the proposed method to solve ship motions, it is recommended to investigate a high-speed craft. For this reason, in this section, a high-speed catamaran is investigated at three different (constant) velocities. Table 16 shows the main characteristic of the catamaran.
A total of 180,992 structural cells are used to simulate the catamaran at three constant velocities (Vel � 5, 16.5, 24) ( Figure 38).
Like the previous test case, the catamaran moves with a constant velocity while it is allowed to exhibit heave and pitch motions (2-DOF). Each simulation is performed until three main hydrodynamic parameters of catamaran become steady: heave motion, trim angle, and resistance force. Figure 39 shows the catamaran in steady conditions at different speeds.
Heave, trim angle, and resistance force are compared with the work of Panahi et al. [28] in Figure 40.
As shown in Figure 40, the presented model could predict the heave and pitch motion with reasonable accuracy. e resistance which is the result of these two parameters is also well calculated. For this test case, the multigrid method is applied using the grids which are listed in Table 17.
e residual reduction for the second time step is plotted in Figure 41.
Like other simulations, all types of cycles for 50 time steps are applied to investigate their performance Table 18. Based on Table 18, W cycle with a four-grid level was chosen for all simulations.

Resistance and Motion of Oceanography Research Vessel.
e final test case involves the calculation of the resistance and motion of an Oceanography Research Vessel (RV) in still water [45]. Table 19 shows the main characteristic of the vessel. e resistance of the vessel is estimated experimentally at the Vienna Model Basin [45]. e wooden model of 50 m RV is made in 1 : 12 scale. To achieve better accuracy, two shaft A-brackets and two rudders are applied to the model as shown in Figure 42.
Based on the towing tank test, the vessel resistance at different speeds is estimated and presented in Table 20.
Based on the geometrical full-scale model, a fully unstructured grid is used for numerical simulation. e discretized RV is shown in Figure 43.
It should also be noted that A-brackets and rudders are modeled in a numerical model like an experimental wooden model for better accuracy Figure 44. ese appendages play an important role in creating hydrodynamic resistance of the ship hull.
Based on the proposed prolongation algorithm, three different grid size levels are created in numerical simulations. Table 21 shows the grid size properties. It should be noted that the numerical solver could not solve the coarse grid level 3 because of details in the appendage part of RV. In this part, the coarse cells become numerically unstable. erefore, in this test case, three grid levels are used in the multigrid cycles.        As observed, in fully unstructured grids, the proposed coarsening algorithm has a better performance. Comparison of the coarsening in an unstructured triangular cavity test case with the structured BFS or dam break grids reveals that, in unstructured grids, the coarsening ratio is about 4.5, while in structured grids, the ratio is almost 4. is is due to the fact that, in structured grids, the number of sharing cells for each vertex is less than the unstructured grids. For example, in the triangular cavity, there are vertexes with 6 sharing cells, while this number is 4 in the BFS test case. On the other hand, it may have a coarsening limitation on unstructured grids. For example in this case, "coarse grid level 3" could not be created because, by applying the presented coarsening algorithm, a very complex cell would be generated which occupies the whole space between rudders and brackets.
In the forward motion of the ships, heave and pitch motions are dominantly related to four other motions. In        Table 13.                 Figure 46 compares the simulated resistance curve with experimental data [45]. e comparison shows good agreement and indicates that the proposed numerical simulation can predict ship resistance with reasonable accuracy.
In displacement ships, with rather a low operative speed, the heave and pitch motions are small. In this case, the pitch angle is almost zero and the heave motion is ignorable. e heave diagram is plotted in Figure 47 for different velocities.
In order to capture free surface elevation with good resolution, the grid sizes are chosen to be finer near the ship and free surface related to far field. Figure 48 shows a free surface near the ship hull at different velocities.
As all the simulations are performed by the proposed multigrid algorithm, the residual reduction curve can be investigated in order to choose the best multigrid cycle. Figure 49 shows the residual curves for different multigrid cycles in the second time step. Table 22 displays the CPU time and memory allocations of different iterative methods for 100 time steps.  As observed in Table 22, the W cycle is proved to be more effective in terms of convergence and CPU time in the simulation of the research vehicleʼs resistance.

Conclusions
An unsteady 3D Navier-Stokes solver combined with the VOF approach is developed to simulate different fluid phenomena with emphasis on the ship motion. Applying a fully unstructured grid to the solver makes it possible to study the ship motions with all hull appendages such as rudder and bracket. To speed up the convergence rate, an agglomeration multigrid method is used with a new geometrical agglomeration technique to create coarse grids. For this purpose, a fully nonstandard unstructured grid with multiple faces is applied to the solver. e new data structure makes it possible to increase the accuracy in coarsened grids, especially in shipʼs hull details like appendages. As the fractional step method is the technique for coupling velocity and pressure field for solving Navier-Stokes equations, the multigrid method is applied to pressure Poisson equations.
is combination with new nonstandard grids is one of the main novelties of the current article. A rather simple prolongation/restriction operator is used to transfer data between coarse/fine grids in both V and W cycles. In order to validate the accuracy and investigate the performance of the proposed method, two distinguishing parameters are studied: computational time and residual reduction. In each test case, these two parameters are compared to ascertain which cycle has better performance.
It is determined that W cycle reduces the residual better than V cycle. is is due to the fact that, in one cycle of W scheme, double data transferring to coarse grids appear. is action helps the solver to damp the lower frequencies, easily. However, this data transferring takes more CPU time compared to V scheme. erefore, one can see in the test cases that W cycle is more suitable for more complicated problems with a large number of grids. Consequently, V cycle is preferred to be used in simple cases like two-dimensional problems, but the most important issue is that the multigrid method can speed up the solver up to 3 times without reducing accuracy.

Data Availability
Data will be made available upon request.

Conflicts of Interest
e authors further declare that there are no conflicts of interest