THE NUMERICAL INVESTIGATION OF WAVE GENERATION BY LOW SPEED WINDS

In this study we investigate the process of wind-wave generation by using direct numerical simulation (DNS). Air and water domains are respectively solved by Navier-Stokes equations in 3D Cartesian coordinates in which air-water coupled boundary conditions are specified at interface. A shear wind on the top of the air domain is specified, and air and water domains are subsequently driven to fully-developed turbulence, allowing wave growth at interface. In this paper we improve the work published by Lin et al. (2008). Instead of simplified linear boundary conditions (BCs), we derive and impose the non-linear BCs for normal stress at interface, as well as non-linear curvature terms used to balance the discontinuity. The results show that at the linear (initial) stage, faster wave growth is seen with non-linear BCs than with linearized BCs. This is reversed during the exponential (developed) stage.


INTRODUCTION
Wind-generated waves are the most commonly observed phenomena when air flow blows over water surface.The process of wave growth is strongly related to the interaction of air and water.However, the detailed mechanism of how waves can be generated, amplified and evolved by winds is still unclear.Understanding how air-water interaction mechanisms accounts for wind-wave generation has been a long time investigation.The earliest theoretical studies have proposed wave growth mechanisms under relatively simple conditions; the sheltering mechanism was the earliest hypothesis proposed by Jeffreys in 1925 (Jeffreys, 1925).Later theoretical works suggest a systematical framework: linear growth during the initial generation (Phillips, 1957(Phillips, , 1977) ) and then exponential growth during the developed stage (Miles 1957;Belcher and Hunt 1993).
As computer speed has increased, numerical-based studies of wind-wave generation have become correspondingly tenable.By solving the Navier-Stokes equations, numerical modeling provides a detailed depiction of growth processes, allowing a fundamental understanding of more complicated flow conditions.Most numerical studies, however, emphasize either wave-dominated effect, i.e., wave effect on air motion; or wind-dominated effect, i.e., wind stress effect on waves.That is, only single phase flow is simulated numerically coupled with pre-specified quantities from the other phase.Only a few 'true' two-phase flow simulations studies have been published.One such study is that of Lin et al. (2008), who first proposed a numerical investigation using direct numerical simulation (DNS).The authors used 3D Cartesian Eulerian grid fixed on space, implying the small amplitude (linear) assumption.Additionally Yang andShen (2011a, 2011b) developed an improved numerical model using generalized curvilinear coordinates, allowing computational grids to move with the water surface.
In this work, we propose an improvement on the model of Lin et al. (2008) on the boundary conditions (BCs).Two of four air-water coupled interfacial BCs represent the force balance due to water surface tension, discontinuity of normal stress BC and continuity of tangential stress BC.Instead of using the linearized approximation, we derive nonlinear formulations for normal stress BC without any simplifications.In this paper, we address the methodology, formulation and model developments in detail.Using the resulting DNS model, we apply a shear wind on the top of air domain to drive both air and water flows.The waves are subsequently generated from an initially calm water surface.To study the numerical results of the wind-wave generation process and the effect of the nonlinear BC, we compare the results from linearized normal stress BC (Lin et al., 2008) with the results from nonlinear normal stress BC developed in this work.The conclusions are discussed in the last section.

WIND-WAVE GENERATION BY DNS Problem Description
As shown in Figure 1, the simulation of wind-wave generation using DNS can be solved by specifying two separate domains which bound an interface.Navier-Stokes equations are respectively solved for air and water domain, with air-water coupled boundary conditions imposed at the interface at each time step.Both domains are solved on fixed three-dimensional Cartesian Eulerian grids.The simulation of wind wave Re (a) and Re (w) are Renolds numbers respectively for air and water domains.For example, Re (a) is defined as where U is the reference velocity, L is a reference scale, and ν (a) is air kinematic viscosity.In this study, we choose U = U 0 and L = h.The system is then forced by applying a constant shear wind U 0 on the top of air domain.The shear wind drive the flow in the air domain; motion is then transferred to the water domain due to interfacial friction, thus allowing wave growth at the water surface.Numerically the force can be described as Dirichlet boundary conditions.The boundary conditions on the top of air domain (z = h) are therefore given by where U 0 = 3(m/s) in this study.To eliminate the effect of the presence of solid bottoms, free-slip boundary conditions are imposed at the bottom of water domain (z = −h) to emulate the infinite depth, which can be written as

Interfacial Boundary Conditions
Air-water coupled boundary conditions are most significant parts in this study.We assume air and water are immiscible.Accordingly, there are 4 types of boundary conditions for immiscible two-fluid interface.
1. Continuity of velocity (no-slip BC): No-slip BC ensures the continuity of velocity at the interface.
2. Kinematic free surface BC: Kinematic free surface BC is used to update water elevation η at the interface in terms of convection velocities (u, v, w) given by Navier-Stokes equations.The updated water elevation will then be substituted into the third interfacial BC (described below).Eventually the pressure BCs at air and water interface can be obtained and fed back to Navier-Stokes equations.Using the time series of water elevation, we can determine all properties regarding wave growth, no matter transient or statistical, for additional study.The following two boundary conditions are derived from the balance of surface tension forces.As shown in Figure 2, given two immiscible fluids, i.e., air and water in this study, the total stress tensors at interface are given by T (a,w) = T (a,w) i j = −p (a,w) δ i j + τ (a,w) i j (7) where p is the total pressure and τ i j is shear stress tensor.Total pressure is composed of hydrodynamic pressure p h and static pressure ρgη, which can be denoted as Considering the interface as a smooth surface, at any point on the surface we have a unit normal vector n, and two unit tangent vectors tx and ty , as well as water surface tension σ.According to our force balance, the following two BCs can be obtained.

Continuity of tangential stress
The two equations describe the continuity of tangential stress.The tangential stress from the air domain projected onto the surface is exactly identical to the tangential stress from the water domain projected onto the same point.This is done since the water surface tension σ is presently assumed constant everywhere in this system.If σ is not a constant value, the term ∇ • σ can be introduced in these two equations.
After fully expanding these two equations, we find numerical difficulties not only involved in treating the coupled conditions (implying inevitable internal iterations between equation ( 9) and ( 10) at each time step), but also involved in boundary value problems with convolution terms in pseudospectral method.To avoid the difficulties, in this study equation ( 9) and ( 10) are solved in linearized forms.This is the only approximation used in this work.In Cartesian coordinates, the linearized equations ( 9) and ( 10) are respectively given by where µ is dynamic viscosity.
4. Discontinuity of normal stress balanced by tension force: This BC shows that when the normal stress from air and from water project onto the same point at the surface, the difference of magnitude is balanced by the curvature and water surface tension force.So far all the BCs are the same as the work proposed by Lin, et al. (2008).However, in this BC, we are solving fully expanded non-linear forms.Therefore, all the numerical results will compare the wave growth with linearized normal stress BC and non-linear normal stress BC.
At the first, equation ( 13) can be expanded by substituting into ( 7) and ( 8), and written in the nondimensionalized form as where ρ (a) and ρ (w) are respectively air and water densities; and Re (w) = ρ (w) UL µ (w) are respectively air and water Reynolds numbers; We = ρ (w) U 2 L σ is the Weber number.
G (w) , G (a) and curvature ∇ • n differ between non-linear and linearized forms (Lin et al. 2008).Numerical results in this study will examine the differences.Note that in the following equations (l) can alternatively be (a) for air domain or (w) for water domain.

Model Development
Numerical MethodsA pseudospectral method is employed as the numerical scheme to solve the model equations.Fourier basis is used to solve x (streamwise) and y (spanwise) derivatives, which implies that periodic boundary conditions are applied along two horizontal directions.The z direction is non-periodic and we use 2-nd order finite difference schemes to solve all vertical derivatives.Additionally, a non-uniform staggered grid is also used along z direction.On-site grids are used for solving u, v, p while off-site grids are used for solving vertical velocity w.Because air and water domain are solved separately, inevitable 'ghost grids' are placed above water domain and below air domain.Since Navier-Stokes equations are solved entirely in each domain, BCs are required on the ghost points, which are the boundary points for each domain.Thus, the BCs derived in the previous section are used to provide the boundary values for the ghost points.
For the configuration of non-uniform computational grids along vertical direction, finer grids are collocated near air-water interface for both domains, while coarser grids are collocated far from interface.The number of grids is (N x, Ny, Nz) = (64, 64, 65).The dimension of computational domain (Lx, Ly, Lz) equals to (6h, 6h, h) in which the reference length scale L = h is equal to 4cm.Reference velocity is set to be U = U 0 = 300(cm/s), so the reference time is 0.01333s.The time step is 0.005 in non-dimensional units, which equals to 6.6667e − 5s in dimensional units.A fractional step method, one of the project methods, is employed to solve Navier-Stokes equations, and low storage 2nd Runge-Kutta method is used for the time-marching scheme.The proposed model is also parallelized using OpenMP and the numerical case is performed using 8 processors.
Initializing the completed systemAccording to our numerical experience, a successful simulation case should proceed in prescribed steps rather than naively triggered ab initio by only the shear wind.A complete simulation is sequentially performed in three stages, which are explained below.
1. Stage I.The first step is to assign the mean velocity profile analytically, and set a constraint for free surface to retain flat interface.Therefore, we spin up the turbulence by adding a buoyancy force in the z-momentum equation for 80 turnover time units (physically 1.0664s).
2. Stage II.At the second stage the buoyancy force in the z-momentum equation is turned off, but we still continue the spin-up simulation for another 2400 large-eddy turnover time units to reach a pure shear-driven state (physically 32s).
3. Stage III.In turn, we release the constraint to allow the flat interface to become freely deformable.
Waves are then generated according to the prescribed air and water flow conditions, as well as the coupled BCs.Therefore, based on the fully developed shear-driven turbulent provided by the previous Stage II, we officially start our simulation whenl the waves reach fully developed state.

RESULTS AND DISCUSSIONS Flow Snapshots
Figure 3 and 4 respectively show the water elevation η and streamwise velocity u at interface at 2.6s, 15.37s, 26.44s, and 66.8s.It can be observed that at 26.44s waves are in the transition state in which random waves become uniform.Before the transition region, wind waves occupy the initial linear growth stage, which is exemplified by highly random waves with shorter crests.After the transition region waves are developed and enter the exponential growth stage in which gathered waves propagate toward uniform direction with long crest.

Wave Growth
The rate of wave growth can be defined as η 2 -mean-squre-of-η, which is the most significant indicator for the state of wind waves.The wave growth results from linearized normal stress BC and nonlinear normal stress BC are compared and shown in Figure 5.It is seen that at the linear growth stage (t < 35s) waves generated using the nonlinear normal stress BC amplify faster than those from linearized normal stress BC.However, the opposite occurs after the transition region: waves generated from the linearized normal stress BC grow at a faster rate at exponential growth stage.Therefore, there is a crossover point found at the transition stage at around t = 37s.
At the linear growth stage, we can compare the numerical results with the theoretical prediction using the formulation proposed by Phillips (1957):   stress τ s , calculated by The results show in Figure 6.(Phillips, 1957) and numerical solutions for linear wave growth at the initial stage(t < 20s).Blue: analytical solution.Red: numerical solution with nonlinear normal stress BC.Green: numerical solution with linearized normal stress BC.

Evolution of Interfacial Properties
Figure 7 shows the evolution of some interfacial air properties.For both linear and nonlinear normal stress BC cases, wind shear stress τ s remains almost constant in the linear growth stage and linearly increases in the exponential growth stage.The friction velocity u * a , which is related to τ and can be calculated by Eq (20), also has the same trend and gives the averaged value u * a ≈ 8.616 (cm/sec It is found that at exponential growth stage, D p for linear BC case is larger than nonlinear BC case, which is consistent with the results of wave growth in Figure 5. Accordingly, it can be concluded that at exponential growth stage using linearized normal stress BC causes over-estimated growth rate and form stress.

CONCLUSIONS
We offer several concluding remarks drawn for the work proposed by this paper: 1.A high-resolution numerical tool has been developed for high Reynolds number problem using pseudospectral method.It is successfully applied for the simulation of air-water coupled two-phase flow.
2. Linear and nonlinear normal stress BCs for windwave generation process have been compared and studied using DNS.
3. From the the results of wind-wave generation, we can conclude that (a) Linear growth (t < 40s) stage: for the growth rate, faster in the case of nonlinear normal stress BC but slower in the case of linear normal stress BC.(b) Exponential growth (t > 40s): for the growth rate, slower in the case of nonlinear normal stress BC but faster in the case of linear normal stress BC.Using linearized stress BC formulation, form stress D p is also over-estimated at this stage.

Figure 2 :
Figure 2: Surface tension at immiscible two-fluid interface

Figure 3 : 2
Figure 3: Water surface elevation η at time t = 2.6s, 15.37s, 26.44s, and 66.8s (top to bottom).Left column: results from linearized normal stress BC.Right column: results from nonlinear normal stress BC

Figure 4 :
Figure 4: Streamwise velocity u at interface z = 0 at time t = 2.6s, 15.37s, 26.44s, and 66.8s (top to bottom).Left column: results from linearized normal stress BC.Right column: results from nonlinear normal stress BC

Figure 6 :
Figure 6: Comparison of analytical solution(Phillips, 1957) and numerical solutions for linear wave growth at the initial stage(t < 20s).Blue: analytical solution.Red: numerical solution with nonlinear normal stress BC.Green: numerical solution with linearized normal stress BC.
Mean surface current U s

Figure 7 :
Figure 7: Interfacial air properties ).The root-mean-square of shear stress fluctuation τ s 2 shows similar trends for both linear and nonlinear cases.One reason why both cases can not be distinguished clearly for shear stress-related properties is because the shear stress BC used is the same as that of Lin et al. (2008).Due to different normal stress BCs, on the other hand, root-mean-square of pressure fluctuation p a 2 and form stress D p show different results between linear BC and nonlinear BC.The stress D p is related to p a and can be calculated by