Numerical strategies to speed up CFD computations with free surface—Application to the dynamic equilibrium of hulls

Abstract This article presents two numerical procedures to speed up computations when dealing with a Reynolds Averaged Navier Stokes (RANS) solver based on the Volume of Fluid (VoF) or multifluid method to treat the free surface. The first one is a time-splitting procedure for the volume fraction equation, enabling the use of larger time steps for the resolution of the flow, without penalizing accuracy. However, these large time steps destabilize the coupling with the ship motion simulation when computing a dynamic equilibrium position in marine applications. The second procedure is therefore a quasi-static approach to solve the coupled problem of dynamic equilibrium. A comparison of these procedures with classical simulations shows that numerical solutions of realistic problems can be obtained up to four times faster.


Introduction
Free-surface capturing methods, such as the Volume of Fluid (VoF) or level-set formulation, have become more and more popular among the CFD developers involved in viscous marine hydrodynamics. The reason for this increasing interest is that these approaches are more robust than those based on a freesurface fitting methodology since no regridding is necessary. Moreover, the merging or breakup of the interface is also handled in a natural way. To achieve computations with VoF-like approach, specific compressive discretisation schemes are used to solve the volume fraction transport equation and keep the sharpness of the interface (Ubbink, 1997;Queutey and Visonneau, 2007). Even if the capacities and the flexibility of such an approach are unquestionable, two drawbacks can be highlighted: The formulation is intrinsically unsteady since the volume fraction is convected by the flow. Up to now, no steady formulation has ever been successful for non-academic test cases, to the knowledge of the authors.
The compressive property imposes a numerically severe Courant number limitation.
It seems inappropriate and wasteful to use such a complete unsteady approach when dealing with physically steady cases. It is all the more a pity since for an implicit solver, the volume fraction equation is the only equation to have such a Courant number limitation! This issue was first underlined by Ubbink (1997) in the conclusions of his PhD thesis: ''the Courant limitation is not insurmountable because it should be possible to apply a technique of sub-cycling where the time step of the main loop is divided into smaller steps in order to advect the volume fractionsy''.
In this paper, to reduce this Courant number limitation, such an original time-splitting (also called time subcycling) procedure for the volume fraction equation is developed and validated for steady state cases. For instance, it enables to increase the global time step while keeping the Courant number constant. The procedure creates no problem for fixed bodies. However, when coupling the flow solver with a resolution of Newton's laws for the ship body, in order to reach the dynamic equilibrium position of a hull, the large time steps lead to a divergent flow-motion coupling, due to the added mass effects (Söding, 2001). To avoid this problem, a quasi-static approach has been implemented to govern the motion of the body in time. This technique has been successfully combined with the time-splitting procedure for the volume fraction, removing the flow/motion instabilities.
Compared to a classical unsteady approach using the resolution of Newton's law, this new numerical procedure to deal with steady cases for hydrodynamics applications enables to reduce significantly the CPU time.
After a brief description of the RANS solver in which this work has been implemented (Section 2), this article outlines the timesplitting procedure (Section 3) as well as the quasi-static approach (Section 4). Then, in Section 5 a test-case is shown demonstrating the capability and the efficiency of these techniques.

The ISIS-CFD flow solver
ISIS-CFD is developed by the EMN group (Equipe Modé lisation Numé rique) of the Fluid Mechanics Laboratory of the Ecole Centrale Nantes. The solver is based on the Finite-Volume method to build a second-order accurate discretisation of the RANS equations. The flow equations are constructed face by face, which means that cells having an arbitrary number of arbitrarily shaped faces can be accepted. This enables the simulation of flows around complex geometries. The velocity field is obtained from the momentum conservation equations and the pressure field is extracted from the mass conservation constraint, or continuity equation, transformed into a pressure equation. In the case of turbulent flows, additional transport equations for modelled turbulent parameters are solved in a form similar to the momentum equations: they are discretised and solved using the same principles.
Free-surface flows are computed through an interface capturing method: the flow phases are modelled with a transport equation for the volume fraction of water c w in each cell: c w ¼1 means that the cell is completely filled with water, c w ¼0 means that only air is present in the considered cell. The interface between air and water is represented by the numerical discontinuity in the solution for c w . The effective flow physical properties (dynamic viscosity m and density r) are obtained from the properties for each phase, (m w , m a ) and (r w , r a ), respectively, for water and air, with the following constitutive relations: Special attention has to be paid to preserve the sharpness of the interface when solving the transport equation of c w . Therefore, this equation is discretised with a scheme using specific antidiffusive properties, which are fully activated when the amount of fluid drained out of a cell during one time step does not exceed the cell volume. This implies a constraint in the local Courant number (the Courant number is an adimensional parameter roughly defined by: DtV=Dx, where V is the velocity through the considered cell, Dx is the size of the cell, and Dt is the global time step of the temporal discretisation). The latter should not exceed 1 and the usual target value retained is 0.3 for 3D flows. Here, the Blended Reconstructed Interface Capturing Scheme (BRICS) scheme is used (Wackers et al., 2011).
An ALE (Arbitrary Lagrangian Eulerian) approach is used to deal with moving bodies (Hirt et al., 1974;Donea et al., 2004;Leroyer, 2004). All configurations of motion (up to six solved or imposed DOF) can be applied. Analytical weighted deformation techniques have been developed to preserve a mesh fitted to the body during its motion (Leroyer, 2004;Leroyer and Visonneau, 2005).

Description
With an interface capturing approach, we have to solve the ALE convection equation for the volume fraction (c w is simply denoted by c in the following): where V is the domain of interest, or control volume, bounded by the closed surface S moving at the velocity U ! d with a unit normal vector n ! directed outward. The time derivative following the moving grid is written d=dt.
Finite-Volume discretisation leads to the following discretised form: where t p and t c mean, respectively, the previous and the current instant, and Dt ¼ t c Àt p is the current time step. c f denotes the reconstruction of the volume fraction at the centre of the face, whereas F U ! and F U d ! represent, respectively, the velocity flux and the grid displacement velocity flux through the considered face S f . Since steady configurations are investigated here, accurate resolution in time is not required, so the basic first order Euler implicit scheme is applied for the time derivative. The spirit of the time-splitting approach is to enhance the fulfilment of the Courant-Friedrichs-Lewy (CFL) condition related to Eq. (2), by using a specific time step for the volume fraction, which is a fraction of the time step associated with the global simulation. In other terms, the global time step Dt is split into a sequence of smaller ones leading naturally to smaller Courant numbers ( Fig. 1). As a consequence, the volume fraction equation is solved sequentially several times during a single global time step, for only one update of the other flow equations; each resolution is called a ''subcycle''. If we note N the number of ''subcycles'', the split time step Dt i is equal to Dt=N. The intermediate volumes V i , for which no mesh is reconstructed, are linearly interpolated between t p and t c . Thus, the convection equation for c becomes Instead of solving Eq.
(2), we solve N times Eq. (3), going forward in time progressively. From one equation with a typical Courant number Co for the classical approach, the time-splitting approach leads to N equations to solve but with typical Courant numbers around Co=N. Finally, by summing Eq. (3) from i equal 1 to N,w e obtain the following Eq. (4), similar to Eq. (2): The free-surface moves little by little during these small time steps. As the CPU time related to the resolution of the volume fraction equation is not large compared with other parts of the solver (especially the computation of the pressure), the global CPU time of the simulation is strongly reduced. When the steady state of the flow is reached, the solution is obviously the same as the one obtained with a classical approach, since the temporal derivative term vanishes. Then, solving either once or several times a convection equation with a null temporal derivative does not change anything. With the time splitting approach, the way to reach a steady state remains physically correct (even if we try to remove the transitory state quickly without solving it accurately). This is probably the reason why this procedure is very robust.

Example
To illustrate the time splitting procedure, a 2D configuration of a submerged NACA-0012 hydrofoil is considered. This wellknown test-case experimentally investigated by Duncan (1983) is often used for the validation of numerical methods developed to compute free-surface flows. In the experimental setup, the hydrofoil, whose chord is c ¼20.3 cm, is towed in a tank with speed U¼0.8 m/s, with an angle of attack a ¼ 51 (see Fig. 2). The relevant non-dimensional parameters based on the chord length and the free-stream velocity are Fr¼ 0.5672 and Re¼1.423e5. The distance between the profile and the bottom of the basin is kept fixed (H¼17.5 cm from the mid-chord of the profile), whereas the depth of submergence s is varying. Here, we are only interested in the case s¼23.6 cm.
The simulations are performed using a mesh of 38 000 cells. The two-equation k2o SST closure (Menter, 1993) is used to take into account the turbulence phenomena. Boundary conditions are imposed as follows: on the top and on the outlet, a hydrostatic pressure is imposed. On the inlet, the velocity is imposed at its far field value. A slip condition is applied on the bottom of the tank, whereas wall-function boundaries are applied on the whole surface of the foil.
Different simulations are performed with various numbers of subcycles N. For all the simulations, the time step is such that the Courant number Co for the split volume fraction equation remains around 0.3. Table 1 summarizes the numerical setup and results in terms of CPU time. As expected, the CPU time is greatly reduced using the time-splitting approach, and the final solutions (forces on the body and free-surface deformation) are nearly identical: Fig. 3 shows that the converged value of the lift is the same for all computations, but the path to convergence is slightly different. The classical computation reaches the mean value quickly but with oscillations, which take a long time to damp out. When the number of subcycles is increased, the physical time to obtain the converged value increases too, but the gain in CPU time remains significant (the CPU time for 10 subcycles is four times shorter than for classical approach). Above N¼ 10, the efficiency of the approach saturates but robustness is kept even if a large break is visible at t % 9 s for N¼20. Fig. 4 confirms that the converged free surface is not influenced by the time-splitting approach.

The quasi-static approach
When dealing with potential flows, Delhommeau (1987) underlined that the loads on a hull can be divided into four parts: external forces (due to mechanical bindings or mooring lines), gravity force, hydrostatic loads, and hydrodynamic loads. A firstorder evaluation of the hydrostatic loads around an equilibrium position can be obtained as a function of the sink dT z , the roll dR x and the trim dR y , as well as geometic and inertia characteristics     (Delhommeau, 1987 The coefficients of the matrix S depend on the density of water (r), the magnitude of gravity (g), the geometry of the wetted surface of the hull (S w ), the immerged volume (V i ), the centre of gravity (coordinates (X G ,Y G ,Z G )) and the vertical position of the centre of buoyancy (Z C ).
By computing S À1 , Eq. (6) can be used to provide the hydrostatic position of any hull through an iterative process, assuming that the mass and the centre of gravity position are known: Coupled to a flow solver, it enables to predict a dynamic equilibrium position. In fact, let us consider L ¼ðF z ,M x ,M y Þ the vertical force and torque (evaluated at the centre of gravity G) along the X and Y-axis, respectively, acting on the hull (including gravity, fluid force and possible external forces), in a given imposed position. By computing the coefficients of the matrix S using this current position of the ship, a prediction of the new position can be obtained by solving the linear system equation (6), giving dP ¼ðdT z ,dR x ,dR y Þ as a result.
This amounts to extrapolating an equilibrium position with a first order method while keeping the hydrodynamic load and the external forces constant. This inferred new position is relaxed and then applied progressively during the kinematic transient period DT i . After some time steps, a new prediction of the equilibrium position using the current fluid forces applied to the hull can be obtained. Then the procedure can go on up to convergence. The procedure is summed up as follows (t e means the time when an evaluation is carried out, DT h is the temporal interval between two evaluations, with DT h ZDT i ): A t e ¼t c , B computation of the current forces L acting on the hull, C prediction of dP to reach the new equilibrium position using Eq. (6), D underrelaxation of the new position with a coefficient aA½0; 1 : P new ¼ P old þ adP, E gradual variation of the position to reach P new within DT i , F when t c Z t e þ DT h ,g ot oA.
5. An application test-case: DTMB 5415 with free trim and sinkage Both the techniques described previously are applied here to the bare hull DTMB5415 test case (see Fig. 5) in the following conditions: Fr ¼0.28 and Re¼1.26e7, trim and sinkage free, L pp ¼3.048 m. This test case was used for the CFD Workshop Tokyo 2005 (CFDWS2005, case 1.3) and is the subject of experimental and numerical studies (see Longo and Stern, 2005;Carrica et al., 2007 for example). A mesh of 1.6 M cells around the half body is used to compute the flow (see Figs. 6 and 7). Boundary conditions are imposed as follows: on the top and the bottom of the fluid domain, a hydrostatic pressure is imposed, with a Neumann condition for velocity. On the inflow, outflow and side boundary, the velocity is imposed at its far field value, as is the volume fraction (in this case, a Neumann condition is prescribed  for the pressure). Mirror symmetry is used for the plane Y ¼0. On the body, a wall function is used everywhere, except on the deck where a slip condition is applied. A classical computation Ref (no subcycling, time-accurate resolution of Newton's laws) and a computation TS ÀQS combining the time-splitting and the quasi-static approaches are performed. For the latter, 10 subcycles are applied, using a time step equal to 0.05 s, i.e. 10 times those of the classical computation. The ship is accelerated from rest to its nominal velocity ( ¼1.531 m/s) in 2 s and 4 s, respectively, for the classical and quasi-static computation. Acceleration ramp follows a quarter of sine law to reach the nominal velocity smoothly. The convergence time is not very sensitive to the duration of this ramp, within a certain range. However, to avoid convergence troubles, it is recommended not to go far beyond an acceleration of 1 m/s 2 for model scale hulls. For the quasi-static computation, we prefer to double this transient stage since the time step is far larger.
The parameters for the quasi-static approach are as follows: a new prediction of the equilibrium position is computed each DT h ¼ 0:3 s and the relaxation value a of the displacement is set to 0.3. While the sinkage is updated from the beginning, the trim is released only after reaching the nominal velocity, i.e. after 4 s. Linear laws for trim and sinkage are imposed to the ship to reach the new predicted position within a temporal interval of DT i . For this case, DT i is chosen equal to DT h , i.e. as soon as the ship reaches its new position, a new prediction occurs, without any latency period in which the boat remains fixed. This latency period (i.e. DT i oDT h ) is needed for stiffer cases like planing hulls, to let the hydrodynamic force stabilize before a new prediction of the equilibrium.
Figs. 8 and 9 show the evolution of trim and sinkage for both configurations. It takes more time to reach the equilibrium position for the simulation TS À QS, but since the time step is multiplied by 10, this is not really a problem. However, contrary to the classical approach in which the oscillations of the motion are difficult to damp out quickly (even using an artificial damping term), the quasi-static approach gives a very stable equilibrium position. The values are in quite good agreement with the experimental data: differences are comparable to those obtained in Carrica et al. (2007). Here again, the free-surface is the same for both results (see Fig. 10). The gain in CPU time is similar to that exposed in Section 3.2.

Conclusions
This article describes first an algorithm related to a timesplitting technique, in which the advance in time is not the same for the volume fraction and for the other variables: the equation of the volume fraction is solved using a smaller time step which reduces significantly the Courant number. A test-case is presented showing that the optimal number of splits is around 10. Then, a quasi-static approach is described to be used instead of the resolution of the Newton's laws to reach dynamic equilibrium. The convergence is carried out through a succession of predicted positions reached with an imposed motion. Compared to the dynamic approach, the transients are easier to remove, the dynamic equilibrium is far less oscillating. As a consequence, a criterion can be added to test automatically several velocities the ones after the others. Moreover, the large time steps associated with the time-splitting technique can then be used without any flow/motion instabilities. As a consequence, both techniques together enable to speed up computations for reaching a dynamic equilibrium position, without penalizing accuracy. A test-case is shown demonstrating that the results are similar to a classical approach, but with a gain of CPU time around 4.
These techniques have been applied on various kinds of hulls. For some of them (like planing boats), for which the resolution of the equilibrium position is more stiff, a latency period is needed between two predictions to let the fluid forces to be stabilized.
This approach can also be used together with an adaptive time step law to the Courant number (Hay et al., 2006). In this case, a target Courant number can be specified for both the split fraction volume equation and the global time step (the ratio of these both values give approximatively the number of subcycles). This avoids thinking about which time step should be used to have enough accuracy for the interface capturing.
The time-splitting method can also be used with a standard small time step to reduce the Courant number and improve the accuracy of the interface capturing. As a matter of fact, in real cases, the criterion of having the Courant number lower than 0.3 is often difficult to reach everywhere. The next step is to investigate this sub-cycling method for unsteady physical problems, even if the accuracy of such an approach needs to be addressed in this case.