SIMULATING BREAKING WAVES WITH THE REYNOLDS STRESS TURBULENCE MODEL

BACKGROUND AND INTRODUCTION In the past two decades, previous computational fluid dynamics (CFD) studies of breaking waves have shown a marked tendency to severely over-estimate turbulence levels, both preand post-breaking (e.g. Brown et al. 2016). Larsen and Fuhrman (2018) have proved that this problem is due to the unconditional instability of twoequation turbulence closure models (both k-ω and k-ε types) in the potential flow core region beneath surface waves. A method for formally stabilizing two-equation models was derived in their work, and it was shown that stabilized two-equation models lead to pronounced improvement in the predicted turbulence and undertow velocity profiles prior to breaking and in the outer surf zone. However, even the stabilized two-equation models in Larsen and Fuhrman (2018) were still rather inaccurate in the inner surf zone (i.e. closer to the shoreline), thus seemingly requiring yet more advanced methods of achieving turbulence closure.


BACKGROUND AND INTRODUCTION
In the past two decades, previous computational fluid dynamics (CFD) studies of breaking waves have shown a marked tendency to severely over-estimate turbulence levels, both pre-and post-breaking (e.g. Brown et al. 2016). Larsen and Fuhrman (2018) have proved that this problem is due to the unconditional instability of twoequation turbulence closure models (both k-ω and k-ε types) in the potential flow core region beneath surface waves. A method for formally stabilizing two-equation models was derived in their work, and it was shown that stabilized two-equation models lead to pronounced improvement in the predicted turbulence and undertow velocity profiles prior to breaking and in the outer surf zone. However, even the stabilized two-equation models in Larsen and Fuhrman (2018) were still rather inaccurate in the inner surf zone (i.e. closer to the shoreline), thus seemingly requiring yet more advanced methods of achieving turbulence closure.
The present study takes the Reynolds stress turbulence model (RSM) in terms of the Wilcox stress-ω model (Wilcox, 2006), as a prime candidate for simulating breaking waves. RSMs resolve all components of the Reynolds stress, eliminating e.g. the assumed isotropy of turbulence inherent within two-equation models. Therefore, RSMs are expected to produce more accurate results for predicting breaking waves as they capture the physics more realistically, comparing to two-equation models. Meanwhile, comparing to the large eddy simulation (LES) and the direct numerical simulation (DNS), RSMs can largely reduce grid numbers and computational costs.
However, it was found in Brown et al. (2016) that using RSM model in terms of the LRR stress-ε model (Launder et al., 1975) significantly over-estimated the turbulent kinetic energy for breaking waves at both pre-and postbreaking zones to an even larger degree than using standard two-equation closures. The present work proves that the reason is the lack of buoyancy production terms in the LRR stress-ε model in Brown et al. (2016). Different from two-equation models, the RSMs (both LRR stress-ε model and Wilcox stress-ω model) are unconditionally stable according to the present stability analysis and progressive wave train tests. The present study proves that the Boussinesq eddy viscosity assumption in twoequation models is the fundamental reason for the turbulence kinetic energy to exponentially grow in the potential flow beneath surface waves.
In the present study, first, novel stability analysis of RSMs in terms of the LRR stress-ε model and the Wilcox stressω model (Wilcox, 2006) is performed. The RSMs are then implemented in OpenFOAM (version v1912) and applied for simulating progressive wave trains to confirm their stability in the potential flow regions. Finally, the Wilcox stress-ω model, which demonstrates a higher stability and better wall treatment, is applied for simulating spilling breaker hydrodynamics. The results are compared with the experiments of Ting and Kirby (1994) and the stabilized k-ω model of Larsen and Fuhrman (2018).

STABILITY ANALYSIS
The present stability analysis of RSMs (both LRR stress-ε model and Wilcox stress-ω model) found that the asymptotic growth rate of turbulence kinetic energy Γ ∞ equals to zero in the potential region beneath surface waves. For both LRR stress-ε and Wilcox stress-ω models, (1) where Γ ∞ is the asymptotic growth rate of turbulence kinetic energy, is a constant in the form of the model closure coefficients, and 12 = 1 2 � + � is the strain rate tensor. Time-and depth-averaged 12 gives ≪ 12 ≫= 0 beneath surface waves. Therefore, the growth rate Γ ∞ is 0 for RSMs regardless of the value of the constant . Conversely, it was found in Larsen and Fuhrman (2018) that for two-equation models, the growth rate of turbulence kinetic energy is always positive in the potential flow region below surface waves. The present study proves that it is the Boussinesq eddy viscosity assumption used in two-equation models: = 2 − 2 3 (2) leading to the shear production of turbulent kinetic energy to be in the form of = 0 , 0 = 2 (5) so that the growth rate of turbulence kinetic energy is in the form of 0 Γ ∞ ′ = 1 ≈ ′� 0 (6) where ′ is 0.125 for standard k-ω model according to Larsen and Fuhrman (2018), and depth-and timeaveraged ≪ 0 ≫ is proved to be always positive for progressive waves. Therefore, the turbulence kinetic energy is growing exponentially beneath surface waves if using two-equation turbulence models, and the wave elevation will decay after propagating for a number of wave periods (Larsen and Fuhrman, 2018). In the equations above, is the Reynolds stress tensor, is mean strain rate tensor, is the turbulent kinetic energy, is shear production term for , is the kinematic viscosity, ε is dissipation rate of and = ε/( ) is the specific rate of dissipation.
In the present study, the RSMs are implemented in OpenFOAM v1912 and applied for simulating progressive waves to confirm their stability in the potential flow regions. Figure 1 shows progressive wave elevations simulated by the standard LRR model and the standard k-ε model in OpenFOAM, both with buoyancy production terms added in the model. The same phenomenon is found when applying the Wilcox stress-ω model and the k-ω model. The standard two-equation models fail to propagate stable wave trains due to exponential growth of the turbulence kinetic energy beneath waves, while the RSMs are showing stable performances for simulating progressive waves. Figure 1 -Wave elevations of progressive waves simulated by (a) standard LRR model with buoyancy production terms and (b) the standard k-ε model with buoyancy production terms. The buoyancy production terms added in the model can be referred to Larsen and Fuhrman (2018).

SPILLING BREAKER SIMULATION
The Wilcox stress-ω model, which demonstrates better stability and wall treatment, is applied for simulating spilling breaking wave hydrodynamics of Ting and Kirby (1994). Figure 2 shows that the Wilcox stress-ω model avoids the over-production of turbulence prior to breaking, different from the standard two-equation model behaviors that showed in Larsen and Fuhrman (2018).
Figure 2 -Snapshot of turbulence levels (here depicted as νt/ν) for the spilling breaking case simulated in OpenFOAM using the Wilcox stress-ω model.
As we mentioned before, even though Larsen and Fuhrman (2018) proposed stabilized two-equation models that predicted accurate turbulence and undertow velocity profiles prior to breaking and in the outer surf zone, their results were still rather inaccurate in the inner surf zone. Figures 3 present the predictions of undertow velocity profiles using the Wilcox stress-ω model, compared to those predicted by stabilized k-ω model in Larsen and Fuhrman (2018) and measured in Ting and Kirby (1994). It can be seen that the present RSM simulations demonstrate high accuracy of undertow predictions from pre-breaking all the way into the inner surf zone. Meanwhile, the stabilized k-ω model fails to accurately predict the undertow velocity in the inner surf zone (Fig. 3 c,d).

CONCLUSIONS
The present study formally proved that the RSMs are unconditionally stable in the potential flow regions and produce more accurate results from pre-breaking all the way into the inner surf zone for simulating breaking waves, especially for the undertow velocity profiles in the inner surf zone which even stabilized two-equation models fails to accurately predict.