The Analytical Expression for the Mass Flux in the Wet/Dry Areas Method

Accurate prediction of the distribution of shear stress is essential for the numerical investigation of the interaction between the wind and water waves on ocean surface since the shear stress plays a key role in this type of interfacial flows. The numerical velocity distribution provided by the computational fluid dynamics should have high accuracy as the shear stress is computed by the derivative of the numerically predicted velocity. The recently developed wet/dry areas method based on the conservative integral form of the Navier-Stokes equations mathematically reveals that the convection terms in the Navier-Stokes equations should be calculated by the areas exposed to the water and air. In this paper, an analytical expression for the mass flux in the wet/dry areas method is derived, the discussion is focused on why this analytical expression should be used for more accurate numerical simulation of interfacial free surface flows.


Introduction
In many practical problems the water is under the conditions of fast moving air.The couplings of air-water in a severe weather is more important than quiet weather for the designs of shipment, structures on the coastal zone and in the ocean.The coupling processes between surface gravity waves and adjacent winds and currents are of global significance for climate.This interaction also plays a significant role in the transportation of sediment, wave growth,. . .and so forth.Air and water interact through the shear stress and pressure in fluid, therefore, an accurate prediction of the shear stress and pressure is a necessity for the various applications.In a CFD simulation, the velocity distribution is predicted by using numerical algorithms, and the shear stress is computed by the derivative of the velocity.In order to produce a correct shear stress, a spurious velocity distribution should be avoided in the numerical simulation since the numerical results with spurious velocity do not provide correct distributions of the velocity, pressure, and shear stress in the fluid.

ISRN Applied Mathematics
The accuracy of numerical schemes used in the numerical simulations of interfacial flows has caused concerns in previous study.For example, it was revealed that the simple average density causes spurious velocities near the interface 1 .Long before that, Rudman 2 found that it is crucial to accurately calculate the convection terms mass flux in the Navier-Stokes equations.For the staggered grid, his method requires that the surfaces of the f-control volume are also the surfaces of the u-control and v-control volumes.To achieve this, the mesh for volume fraction f must be twice as fine.Rudman's method was extended to 3D and other recent two-phase flow problems 3-5 .Bussmann et al. 6 pointed out that when the density ratio is large a careful calculation of the convection terms in the Navier-Stokes equations is needed.They developed a numerical model based on a collocated mesh by extending Rudman's method.For a standard staggered mesh, how to accurately calculate the mass flux is still a question.Recently, the author in this paper developed a numerical method for interfacial free surface flow on a standard staggered grid 7 .In this paper an analytical expression for the mass fluxes passing through the surfaces of the uand v-control volumes is derived.The analysis in this paper reveals that the analytical formula produces unique and accurate value for the mass fluxes, whereas the widely used average-density method produces large error for the mass flux and spurious velocity.
In Section 2 of this paper, the governing equations in the conservative integral form are presented and discretized by following the standard discretization procedure of the control volume method.In Section 3, the focus is concentrated on the calculation of the mass fluxes passing through the surfaces of the control volume and a simple but accurate formula for the mass flux is derived.In Section 4, the numerical results and comparisons are presented.

Governing Equation and Discretization
The governing equations in conservative integral form are given in 8 .For incompressible fluid and a control volume with fixed shape, Ω, the continuity equation for volume conservation is written as

2.1
The continuity equation for mass conservation is The momentum equation for u is given by and the momentum equation for v is given by where m Ω ρdΩ is the total fluid mass within the control volume.The equation for volume fraction is given by where v ui vj is the velocity vector, n n x i n y j is the unit vector outward from the surface of control volume, p is the pressure, ρ is the density, μ is the viscosity, g is the gravitational acceleration, S is the surface of control volume, Ω is the volume of control volume, and f is the volume fraction of water.The staggered grid arrangement for u, v, p, and f are shown in Figure 1 a .The control volume for the velocity component u is shown in Figure 1

2.7
If the surface tension is not included, velocity u, pressure p, and shear stress μ ∂u/∂n are continuous functions in the solution domain and on the interface, and then they may be assumed to be constants on each face of the control volume.
Following the calculations described in Patankar 9 and Shyy 10 , we subtract 2.6 multiplied by u P from the left-hand side of 2.7 , and we finally get the algebraic equation for the velocity u as follows:

2.10
Term u b contains the terms from the high-order scheme for the convection terms, if u b is zero then we have the first-order upwind scheme.Similarly applying momentum equation 2.4 and continuity equation 2.2 to the v-control volume we have the algebraic equation for velocity v as follows:

The Wet/Dry Areas Method
Now we have to calculate the mass fluxes F e , F w , F n , and F s given by 2.10 .In order to illustrate numerical method, we now focus on the mass flux on the northern face shown in Figure 2. It is a common practice that the mass flux on face A n is calculated by where ρ n is the average density on face A n , Δx is the area of face A n .When the staggered grid is used, the interpolation techniques have been used to calculate the average density ρ n , but different interpolation methods may produce different values for the average density ρ n .Thus the question is how the mass flux on the surface of control volume can be accurately calculated.This question is answered by the following analysis.
In order to accurately calculate the mass flux we go further to write 3.1 as where ρ w is the density of the water, ρ a is the density of the air, L a is the dry area exposed to the air, L w is the wet area exposed to the water, and Δx is the area of face A n .Equation 3.2 is an analytical, therefore accurate, expression for the mass flux passing through northern face of the control volume.It reveals an important principle that the mass flux passing through the faces of control volume should be calculated by the area exposed to the water wet area and the area exposed to the air dry area .We also see that the density used in 3.2 should be the density of water ρ w and the density of air ρ a not the average density ρ n in 3.1 .Figure 2 and the integration in 3.2 also reveal that when the whole area of A n is exposed to the air, we have L a Δx, L w 0; when the whole area of A n is exposed ISRN Applied Mathematics to the water, we have L w Δx, L a 0; when a part of A n is exposed to the water and another part is exposed to the air then neither L w nor L a is zero but L w L a Δx.In a two-dimensional problems L w and L a are the lengths on the edges of the control volume, for a three-dimensional problem, they are the areas on the surfaces of the control volume.Equation 3.2 was obtained by a slightly different approach in 7 .
It must be emphasised that 3.2 is mathematically accurate since term ρ w L w ρ a L a in that equation is naturally produced by the conservative integral form of the governing equations, and most importantly 3.2 states that the wet and dry areas should be calculated in the numerical discretization process and the densities in these equations are known water and air densities which are very different from the average-density ρ n in 3.1 .
It should be also emphasised that 3.1 and 3.2 produce very different values for the mass flux.For example, in the case shown in Figure 3, if f 1 0.9, f 2 0, the air density ρ a 1 and the water density ρ w 1000, then the average volume fraction f 0.5 f 1 f 2 0.45 and the average density ρ n fρ w 1 − f ρ a 0.45 × 1000 1 − 0.45 × 1 450.55, the mass flux F n 450.5v n Δx; whereas using the wet/dry areas method we have L w 0 and L a Δx, thus from 3.2 the mass flux F n v n Δx.From this example we have seen that the difference between the mass fluxes produced by the average density method and the wet/dry areas method is very large and clearly that the average density method has produced too large error for the mass flux.From 2.7 , we can see that the time rate of change of the momentum on the control volume and the total momentum fluxes on the surfaces of control volume on the left hand of the equation are balanced by the total pressure forces and total shear stresses acting on the surfaces of the control volume on the right hand of the equation .Therefore, the large discretization error in mass flux is fed into momentum equation 2.7 and this causes large errors in the pressure and shear stress.The numerical example next section will show this large error introduced in the average density method inevitably causes spurious velocity distribution near the interface.
When applying the wet/dry areas method we have to calculate the areas L w and L a .For a two-dimensional problem, the formula for L w and L a were derived in 7 by reconstructing the interface based on the solution of the volume fraction f.The SIMPLEC algorithm is applied to obtain the velocity and pressure distribution.The solution of the volume fraction f in 2.5 is obtained by a high-resolution scheme CICSAM Compressive Interface Capturing Scheme for Arbitrary Meshes 11 .

Application to the Progressive Waves in Two-Layer Viscous Fluids
In order to validate the numerical method, the wet/dry areas method is applied to the numerical simulation of the progressive wave of the air-water flow when the wind speed is zero shown in Figure 4.

Initial and Boundary Conditions
Although the solution of the linear wave is not the solution of the Navier-Stokes equations, it produces a reasonable approximation to the problem.In the solution of a linear progressive wave in a two-layer inviscid fluid 12 , the surface elevation, η, of the wave travelling in the x-direction is η a sin kx − σt , 4.1 The velocity components, pressure, and volume fraction in the air are given by The dispersion equation takes the following form: The initial surface elevation of the water wave, velocity distributions in the air and water and volume fraction f are evaluated by substituting t 0 into 4.1 -4.9 .When t > 0, the location of the interface, the velocity for the water and air and volume fraction f at the inlet are calculated by substituting x 0 into 4.1 -4.3 , 4.5 , 4.6 -4.7 , and 4.9 .At the top of the domain, a symmetric condition for the velocity is applied.At the bottom and on the step beach, a no-slip wall condition is applied.At the outlet, the normal derivative of velocity is set to be zero.

Simulation Parameters
The numerical simulations are performed for an air-water system on a laboratory scale.The water and air depth was set at h 0.12 m and h 0.78 m, respectively.The wavelength for the water wave is set at L 0.2 m, wave amplitude a 0.006 m which leads to a deep water wave with h/L 0.6 and a large wave steepness 2a/L 0.06 ak 0.1885 .The wave phase speed, c, is calculated by 4.11 .The length of the computational domain is 4 m and the length of the main test domain is 3 m.The length of step-shape beach is 0.65 m with an average slope of 1/20.The values of the parameters used in the simulation are g 9.8 m/s 2 , water viscosity μ w 1.0 × 10 −3 Ns/m 2 , air viscosity μ a 1.0 × 10 −5 Ns/m 2 , ρ w 1000 kg/m 3 , and ρ a 1.0 kg/m 3 .
A rectangular nonuniform Cartesian mesh is used in the whole solution domain.In order to resolve the boundary layer at the bed Δy 0.0001, 0.0003, and 0.0005 m were tested for the first grid size.It is found that Δy 0.0001 and 0.0003 m produce accurate solution, then Δy 0.0003 m is used.Above the first grid, the grid size increases to the middle of the water wave where the grid size starts to decrease.Near the wave, uniform grids with 8, 12, 16 and 20 grid points were tested to cover the wave height.It was found that both 12 and 16 grid points are sufficient to produce accurate results and, therefore, 16 grid points are distributed vertically around the wave height.Above the peak of the wave, the grid size gradually increases to the top of the domain.A total of 60, 90, and 120 grid points in y direction were tested, and 90 grid points was found to be sufficient to produce accurate results.In the horizontal direction, a uniform mesh is used with the grid size of Δx 3Δy wave , where Δy wave is the grid size in the wave height.The time of numerical simulations lasts until the wave generated at the inlet at t 0 propagates to the middle of the beach.The time step varies from 0.00002 s to 0.00005 s when the wind speed is increased from U 0 to U c. It was found that further increase of the wind speed leads to difficulties in convergence, and the time step has to be reduced, leading to expensive calculations.

Results and Discussion
When the wind speed U 0, Figures 5, 6, and 7 show the profiles of the surfaces of water waves and velocity vectors produced by the potential flow, the wet/dry areas method and the average density method in one typical wave length.In the solution of the potential flow shown in Figure 5 the velocity distribution in the water and air is the classic orbiting velocity of potential flow.In the solution of the viscous fluid shown in Figure 6, the velocity distribution produced by the wet/dry areas method is very similar to that of the potential flow in most of the solution domain, especially in the water.The major difference between the solutions of the potential flow and the viscous flow takes place in the air near the water surface where the wet/dry method produces two re-circulations, one is just above the wave peak, another is just above the wave trough.However, the velocity distribution produced by the average density method shown in Figure 7 is very different from those in Figures 5 and  6.Most difference is that the average density method produces a strong jet flow just above the wave trough and a recirculation just under the wave trough.The recirculation under the wave trough is a spurious velocity distribution since it is impossible for the water to rotate when the wind speed is zero.
In order to quantitatively compare the results the u-velocity varying with the ycoordinate shown is plotted at different locations x/L 0.25 and 0.5.At the wave peak, x/L 0.25, Figure 8 shows that the positive u-velocities of the potential and the viscous flows increase to maximum from the bed to the water surface.Near the water surface of the wave peak the water velocity of the viscous flow is slightly higher than the potential flow.Across the interface, the u-velocity of the potential flow switches to a negative maximum within a distance of zero because of the absence of viscosity.In contract, the u-velocity of the viscous flow is continuous but rapidly decreases to zero from the interface to the centre of the recirculation in the air and reaches a negative maximum at the position above the centre of the recirculation.Above this, the u-velocity produced by the potential flow and viscous flow decrease to small values.A remarkable feature is that over most of the profiles the difference between the velocity magnitudes of the potential and viscous flows is less than 1%.However, just above the wave peak, the average density method produces very strong jet flow with large positive velocity which is several times higher than those produced by the potential flow and the wet/dry areas method.At the wave trough, x/L 0.75, Figure 9 shows that both potential flow and wet/dry method produce negative u-velocities in the water.The negative u-velocities increase to a maximum on the water surface.In contrast, the average density method produces a positive u-velocity under the water surface because of the recirculation under the water surface.On the air side of the interface the u-velocity of the potential flow switches to a positive maximum from a negative maximum, whereas the u-velocity produced by th wet/dry areas method changes continuously from the negative maximum to a positive maximum at the position above the centre of the recirculation in the air.The difference between the potential and viscous flows is also small over most of the profiles.Again the average density method produces very large velocity just above the water.
In order to compare the analytical solution of the potential flow with the solution of the viscous flow produced by the wet/dry areas method the streamlines defined by ψ/ ack are plotted in Figure 10.It is observed that these two solutions have almost identical flow patterns in the water, but the streamlines of the viscous flow are centred inside the air near the peak and trough, whereas the streamlines of the potential flow are centred on the interface at the peak and trough.A remarkable difference between the potential flow and viscous flow is that there is no recirculation in the potential flows but there are two recirculations in the air above the peak and trough of water wave.Figure 11 shows the streamlines and velocity vectors produced by the average density method.It shows that the strong jet above the water surface is very different from the velocities produced by the potential flow and the wet/dry areas method, especially the recirculation under the water surface is clearly visible, indicating the spurious flow pattern.The nondimensional vorticity is defined by ω ∂v/∂x − ∂u/∂y 2akσ .

4.13
The vorticity distribution in the water and air is presented in Figure 12.It is observed that the magnitude of the positive and negative vorticity is much larger in the air than in the water.At the peak of the water wave the vorticity has its maximum at the centre of the recirculation in the air due to the anti-clockwise rotation revealed in Figure 6, whilst at the trough of the water wave the vorticity has its minimum at the centre of the recirculation in the air due to clockwise rotation revealed in Figure 6.Along the vertical direction up into air and down into water the vorticity reduces to zero in a oscillatory manner.This oscillatory variation can be clearly observed in the air where there is a band of yellow colour positive vorticity area and where μ is computed by the harmonic mean of the viscosity of water μ w and the viscosity of air μ a .
Figure 13 shows the distribution of the shear stress in the water and air.The shear stress under water is positive under the peak of water wave and negative under the tough of water wave.The maximum positive and negative shear stresses are just under the water ISRN Applied Mathematics surface.The future of the shear stress distribution under the water is very similar to that revealed in 13 , but the shear stress in the air is much smaller than in the water and the distribution of the shear stress in the air is more complicated than in the water.Well above the water wave, there are two regions: the region above the wave peak has positive shear stress, whilst the region above the wave trough has positive shear stress.

Conclusions
An analytical, therefore accurate, expression for the mass flux passing through the faces of control volume is derived in this paper based on the conservative integral form of the Navier-Stokes equations.This expression produces unique and correct mass flux passing through the surface of control volume.It can be used to check any numerical algorithm.If an algorithm produces different mass flux, then the algorithm contains error.Since the convection terms in the Navier-Stokes equations are accurately computed, the spurious velocity is avoided.The theoretical analysis and numerical example in this paper reveal that the average density method produces an average density with very large error.This error in the density leads to large error in the mass fluxes on the surfaces of the uand v-control volumes, consequently, the spurious velocity around the interface is produced.The analysis in this paper reveals why the average density method is not accurate and the numerical results confirm the spurious velocity phenomenon revealed by Wemmenhove et al. 1 .The wet/dry areas method has been validated by various problems 7 and the accuracy, mesh dependence, robustness and efficiency have been tested in 7 .The wet/dry areas method is a two-phase fluid model, its all equations are derived with generality, therefore, the method can be applied to general interfacial free surface flows.The results of the water-air flow with wind are under preparation for publication.

Figure 1 :
Figure 1: a Staggered location for U, V , P and f, b U-control volume.

Figure 2 :
Figure 2: Northern face is exposed to the water and air.

Figure 3 :
Figure 3: The wet and dry areas when face A n is exposed the the air.

Figure 4 :
Figure 4: Flow pattern of progressive wave in the ocean surface.

Figure 5 :Figure 6 :
Figure 5: Velocity distribution produced by the potential flow, the wind speed U 0.

Figure 7 :Figure 8 :
Figure 7: Velocity distribution produced by the average method, the wind speed U 0.

Figure 9 :
Figure9: u-velocity profile along the y-axis direction produced by the potential flow dotted lines , the wet/dry areas method solid lines , and the average density method dashed lines at location x/L 0.75 wave trough .

Figure 10 :Figure 11 :
Figure 10: Streamlines produced by the potential flow dotted lines and the Wet/Dry Areas Method solid lines , the wind speed U 0.