BOUSSINESQ MODELLING OF NEARSHORE WAVES UNDER BODY FITTED COORDINATE*

FANG Ke-zhao

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China

Key Laboratory of Coastal Disaster and Defence, Ministry of Education, Hohai University, Nanjing 210098, China, E-mail: kfang@dlut.edu.cn

ZOU Zhi-li, LIU Zhong-bo, YIN Ji-wei

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China

BOUSSINESQ MODELLING OF NEARSHORE WAVES UNDER BODY FITTED COORDINATE*

FANG Ke-zhao

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China

Key Laboratory of Coastal Disaster and Defence, Ministry of Education, Hohai University, Nanjing 210098, China, E-mail: kfang@dlut.edu.cn

ZOU Zhi-li, LIU Zhong-bo, YIN Ji-wei

State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China

A set of nonlinear Boussinesq equations with fully nonlinearity property is solved numerically in generalized coordinates, to develop a Boussinesq-type wave model in dealing with irregular computation boundaries in complex nearshore regions and to facilitate the grid refinements in simulations. The governing equations expressed in contravariant components of velocity vectors under curvilinear coordinates are derived and a high order finite difference scheme on a staggered grid is employed for the numerical implementation. The developed model is used to simulate nearshore wave propagations under curvilinear coordinates, the numerical results are compared against analytical or experimental data with a good agreement.

Boussinesq equation, waves, body fitted coordinate

Introduction

The several decade’s development of the Boussinesq wave theory with the help of modern computer technology shows that, the Boussinesq-type wave models are effective and attractive tools to simulate wave propagations in nearshore regions. Since the pioneer work of Peregrine[1], a variety of numerical models have been developed with satisfactory results as compared against analytical solutions, or laboratory and field data, and recent literature reviews can be found in Ref.[2,3].

Most of the numerical implementations of Boussinesq-type wave models invoke the use of the Finite Difference Method (FDM) on rectangular grids[4,5]. Though it is easier for understanding and programming, some difficulties would arise with respect to the following two situations: (1) to deal with the irregular computation boundaries, (2) to enhance models’ efficiency by using different levels of griddiscretization density. One would come across these two situations in the applications of a nearshore wave model. Firstly, in the nearshore region, the shoreline geometries are usually irregular and the computational accuracy would be impaired when they are coarsely treated as stair-stepped boundaries by rectangular grids. On the other hand, the uniform grids would mean a same resolution level (as measured by PPW-points per wavelength in convention) in the entire computation domain. In dealing with a broad spectrum of wind waves, the time-domain Boussinesq model with uniform grids may become too expensive to be used in large near-shore regions, and waves would be over-resolved in deep water and under-resolved in shallow water[6].

The Finite Element Method (FEM) or the Finite Volume Method (FVM) based numerical solvers to Boussinesq-type equations adopt some approaches to overcome these problems such as to use unstructured grid discretization and hence they are more versatile than the finite-difference method in dealing with geometrically complex domains, e.g., FEM based models developed by Sorensen et al.[6]and Zhao et al.[7,8]and FVM based models developed by Shiach and Mingham[9-11]. However in terms of programming and computational cost, FEM and FVM are not superior to FDM. Especially higher order nonlinear terms or dispersive terms, appearing in most newly derived higher order Boussinesq equations and playing a key role in describing the nearshore wave evolution, are quite hard to be numerically dealt with and are usually neglected. Contrarily, FDM with structured curvilinear grids has the advantage of both the geometric flexibility of the FEM and the simplicity of the FDM and thus seems to be a good choice for these problems[9,12,13]. However, in most of the models listed above (except for Shi et al.’s model) only a weak nonlinearity is considered, which will prevent the models from accurately predicting the wave motions with relatively strong nonlinear characteristics.

In the present study, a curvilinear wave model is developed, based on a set of fully nonlinear Boussinesq equations. The contravariant velocity components are employed as the independent variables and the governing equations in curvilinear coordinates are derived. Then this set of equations are discretized on a staggered gird system and following Chen et al.[14], Shi et al.[9], a higher order finite difference scheme is developed to numerically solve the equations. The accuracy of the program is verified by using it to simulate a Gaussian hump of the water oscillation in a closed wave basin. Then, the present curvilinear model is used to deal with complex boundaries and for the gird refinement in the interested region.

1. Mathematical model

1.1 Governing equations in curvilinear coordinate

The second order Boussinesq equations with fully nonlinear characteristics derived by Zou[15]are expressed in terms of the free surface elevation η and the depth-averaged velocity u as

where the subscript t denotes the partial derivative with respect to t,?=(?x,?y) is the 2-D gradient operator, h is the still water depth andd=h+η is the local water depth, g the gravity acceleration.B1and B2are two free parameters B1=29/885,B2= 2/59. The above equations involve a Padé[2,2]approximation of the exact linear dispersion and the shoaling property is suitable for wave propagations in a medium water depth. Specifically, all the nonlinear terms of second order are retained, therefore, Eqs.(1)-(3) have a fully nonlinear characteristics.

The coordinate transformation in the general form is expressed as

where ξ1, ξ2are the new independent coordinates in the transformed image domain. A sketch of the coordinate transformation is presented in Fig.1, where the irregular computation grids in the physical domain are transformed into the regular ones in the image domain with corresponding transformations of boundaries.

Fig.1 The sketch of the transformation

In the curvilinear coordinates the contravariant components of the velocity (u1,v1) vector are employed as the dependent variables to make the implementation of lateral boundary conditions more convenient[9]. The contravariant components of the velocity vector are defined as

where ui,u (i=1,2) are contravariant and covai riant components of the velocity, respectively. The following relations hold in the tensor space for a scalar variable f or a vector variable u

andijg is the contravariant metric tensor and could be calculated viais the Kronecker delta).in Eq.(8) is defined as

is called Christoffel symbol of the second kind.

By using the above expressions, Eqs.(1) and (2) are written as

in which k,l,m=(1,2) denote the first or second components of a variable, for example, x1=ξ1, x2= ξ2and u·u is given by the last expression in Eq.(8).

1.2 Numerical schemes

The governing Eqs.(12)-(13) will be rearranged for a clear numerical implementation. And the contravariant velocity components u1and v1will be denoted simply as u and v by dropping the superscripts hereafter for simplification.

The mass conservation Eq.(12) may be expressed as

The momentum Eq.(13) in1ξ can be written as

Similarly, the2ξ direction momentum Eq.(13) can be written as

The numerical schemes are similar to those in Shi e t al.’s model[9]. Eqs.(15)-(26) are discretized on a staggered gird with scale variables (e.g., η) defined in the gird center, u locates at Δξ1/2 ahead of ηwhile v locates at Δξ2/2 ahead of η. A five point difference formula is used to approximate the first-order spatial derivative terms, for a fourth-order accuracy. While dispersive terms are finite-differenced only to the second-order accuracy. The fourthorder Adams-Bashforth-Moulton predictor-corrector scheme is employed to perform time marching. The predictor step is the third-order explicit Adams-Bashforth scheme, given by

whileth e fourth-order Bas hforth-Mount scheme is given by

where f could be interpreted asη,u,vandΔf is set tobe 0.0001 in the present study.

1.3 Boundary conditions

The vertical solid wall is assumed to enclose the entirecomputation domain, and along these boundaries, the following conditions hold

where n is the unit normal vector pointing outward of the wall. Besides, the sponge layers are also placed in front of the solid wall when necessary to achieve the transmissive boundary condition and the internal wave generation method are also implemented in the present study, with the numerical implementation following the work of Shi et al.[9].

2. Numerical results

2.1 Wave sloshing in arectangle tank

The 2-D model shown in Eqs.(1)-(3) involves a number of mixed derivative terms (especially in higher order dispersive terms and nonlinear terms), and they inevitably become very lengthy after being expressed in curvilinear coordinates. The complexity of the 2-D equations hence requires more careful programming and verifications to avoid possible coding errors. One simple yet efficient testing case used for verification is the evolution of waves in a rectangular basin bounded by four vertical side walls, as shown in the first panel in Fig.2, where an initial wave hump with Gaussian-type shape is released to evolve in a closed basin under the force of gravity. The surface elevation should be symmetric about x and y axes due to the symmetric boundary and initial conditions. Since no water is allowed to escape from the basin, the total water volume inside the domain should be conserved. The symmetry and conservative properties will be used to examine the corresponding results obtained from the model. The computation domain is 7.5 m× 7.5 m and is discretized under the curvilinear coordinates, as shown in the bottom panel in Fig.2. The water depth is 0.45 m and the maximum value in the initial Gaussian shape surface is 0.045 m, the simulation duration is 30 s with time step of 0.01 s.

Fig.2 The sketch of the initial wave surface and computation domain

The contour maps for the surface elevation at two instantaneous moments from the simulations are shown in Fig.3. A short time after the release of the wave hump, i.e., at t=2s, the collapsed water is symmetric about x and y axes. After a long time evolution, t=20s, this symmetric property stillapparently holds, indicating that the finite difference solver under the curvilinear coordinates to discretize the corresponding mixed derivative terms is consistent both inx and y directions. The mass conservation is checked by evaluating the relative error between total mass during time elapses to that at the initial stage, and the results demonstrate that this error is negligible as its value is within 0.023%.

Fig.3 Surface elevation contour at t=2s and t=20s

Fig.4 The circular channel and the grid mesh in the physical domain

2.2 Wave propagation in a circular channel

To demonstrate the present model’s capability in handling horizontally varying geometries, weconsider the wave propagation in a circular channel, where the combined effects of reflection and diffraction areincident from the right bottom into a circular channel with constant water depth h=4m and the inner radius r=75m and the outer radius R=200m.

Fig.5The instantaneous surface elevations at four instants

The physical domain has been discretized int o a 500×125-grid mesh, as shown in Fig.4. The grid sp acingin the radial direction is 1 m, while a constanta ngular grid is used along the channel length, resulting in a maximum tangential grid size of 1.26 m near the outer wall and a minimum of 0.47 m near the inner wall. Plane waves enter the channel from the bottom right and propagate around the bend in a counter-clockwise direction. An internal wave generation source is located at the right end of the channel and two sponge layers are placed in the straight channels extending from the right end to the left end of the circular channel (not shown in Fig.4). The incident wave period T, and the wave height H are 4 s and 0.025 m, respectively. This very small ratio H/h is chosen for the comparison between the numerical results and the linear analytic solution. The time step is chosen as 0.05 s in the present case.

Fig.6 The comparisons between the linear analytical solutions and the numerical results

To illustrate the process of wave propagation in the ch annel, Fig.5 shows the transient surface elevation of a wave train into the channel. It is shown that the waves initially propagate in a straight line (t= 20s), but as the channel bends, the waves startto diffract around the bend and simultaneously run into the curving channel sidewall and are reflected around the bend (t=80sand 140 s), the combined effects of diffraction and reflection finally induces a complex wave pattern in the channel (t=200s). To quantitatively show the accuracy of the model, the comparisons of the water surface variation along the outertion are shown in Fig.6 after a periodic steady state is achieved (where the angleα is defined counterclockwise from the right bottom). Good agreements are found.

Fig.7 The sketch of computation domain

2.3 Wave propagation over a shoal

A physical experiment for the wave transforma-experimental data provide a good test case for model validation. Figure 7 shows the plan view of the wave basin and the transects of the wave gage locations. The computation domain is 25 m×18.2 m with 321×257 grid points as shown in Fig.7. A finer grid was used in the shoal region to better resolve the wave focusing phenomenon. The number of gird points would have to be increased four times if a uniform grid is used with a grid size equal to that at the shallowest part of the shoal. The incident wave was sinusoidal, with wave height of 0.0118 m and wave period of 1.0 s. The time step is 0.01 s. In the numerical simulations, the wave field reaches a steady condition after running the model for 40 s.

Fig.8 The comparison of wave height bettwen computed results and experimental data

The wave height distribution along seven different transects are shown in Fig.8 and it is seen that the numerical results are in good agreements with the experimental data. Along the longitudinal section A-A, the model predicts well the wave shoaling and focusing and the decrease of the wave height after the shoal. The strong wave focusing over the shoal is quite clear, and as a result, the maximum wave height reaches about 2.7 times the incident wave. The shoaling and focusing process could also be seen from sections G-G, F-F to E-E. Note that the wave field behind the shoal is the combined effects of refraction and diffraction, the Boussinesq model also correctly predicts the transverse variation of this complex wave field, as seen at transects D-D, C-C and B-B.

In this study, a Boussinesq-type wave model is developed withfully nonlinear properties accurate to the second order in generalized coordinates. Contravaiant velocity components are adopted as the independent variables in the curvilinear coordinates and the corresponding governing equations are derived. After being discretized on a staggered grid system, the equations are solved by using a higher order finite difference formula with a composite fourth order Adams-Bashforth-Moulton time integration. Three examples are simulated and the corresponding numerical results are in good agreements with the available analytical or experimental data, which demonstrates the model’s capability and efficiency in dealing with the complex boundary and in the grid refinement in the interested region.

3. Conclusion

[1]PEREGRINE D. H. Long waves on a beach[J]. Journal of Fluid Mechanics, 1967, 27: 815-827.

[2]LAKHAN V. C. Advances in coastal modeling[M]. New York: Elsevier, 2002, 1-41.

[3]MA Q. W. Advances in numerical simulation of nonlinear water waves[M]. Singpore: World Scientific Publishing, 2010, 245-285.

[4]CHEN Q., KIRBY J. T. and DALRYMPLE R. A. et al. Boussinesq modeling of longshore currents[J]. Journal of Geophysical Research, 2003, 108(C11): 2601-2618.

[5]FANG K. Z., ZOU Z. L. and LIU Z. B. Numerical simulation of rip current generated on a barred beach[J]. Chinese Journal of Hydrodynamics, 2011, 26(4): 314-320(in Chinese).

[6]SORENSEN O. R., SCHAFFER H. A. and SORENSEN L. S. Boussinesq type modeling using an unstructured finite element technique[J]. Coastal Engineering, 2004, 50(4): 181-198.

[7]ZHAO M., TENG B. and CHENG L. A new form of generalized Boussinesq equations for varying water depth[J]. Ocean Engineering, 2004, 31(16): 2047-2072.

[8]SHIACH J. B., MINGHAM C. G. A temporally second-order accurate Godunov-type scheme for solving the extended Boussinesq equations[J]. Coastal Engineering, 2009, 56(1): 32-45.

[9]SHI F. Y., DALRYMPLE R. A. and KIRBY J. T. et al. A fully nonlinear Boussinesq model in generalized curvilinear coordinates[J]. Coastal Engineering, 2001, 42(4): 337-358.

[10]CIENFUEGOS R., BARTHELEMY E. and BONNETON P. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesqtype equations. Part I: Model development and analysis[J]. International Journal for Numerical Methods in Fluids, 2006, 51(11): 1217-1253.

[11]CIENFUEGOS R., BARTHELEMY E. and BONNETON P. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesqtype equations. Part II: Boundary conditions and validation[J]. International Journal for Numerical Methods in Fluids, 2007, 53(9): 1423-1455.

[12]LI Y. S., ZHAN J. M. Boussinesq-type model with boundary-fitted coordinate system[J]. Journal of Waterway, Port, Coastal Ocean Engineering, 2001, 127(3): 151-160.

[13]ZHOU Hua-wei, ZHANG Hong-sheng. A numerical model of nonlinearwave propagation in curvilinear coordinates[J]. China Ocean Engineering, 2010, 24(4): 597-610.

[14]CHEN Q., KIRBY J. T. and DALRYMPLE R. A. et al. Boussinesq modeling of wave transformation, breaking and runup. II:2D[J]. Journal of Waterway, Port, Coa- stal Ocean Engineering, 2000, 126(1): 48-56.

[15]ZOU Zhi-li. Water wave theories and their applications[M]. Beijing: Science Press, 2005, 597-610(in Chinese).

[16]FANG Ke-zao, ZOU Zhi-li. A 2D nuemrical wave model based on fouth order fully nonlinear Boussiensq equations[J]. The Ocean Engineering, 2011, 29(1): 32-39(in Chinese).

September 26, 2011, Revised January 12, 2012)

* Project supported by the National Natural Science Foundation of China (Grant Nos. 51009018, 51079024), the Founds for Creative Research Groups of China (Grant No. 50921001), the Key Laboratory of Coastal Disaster and Defence, Ministry of Education, Hohai University (Grant No. 200803) and the State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology (Grant No. LP1105).

Biography: FANG Ke-zhao (1980-), Male, Ph. D., Lecturer