Lagrangian block hydrodynamics for environmental fluid mechanics simulations

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

627

\

2010, 22(5), supplement :650-655 DOI: 10.1016/S1001-6058(10)60009-1

Lagrangian block hydrodynamics for environmental fluid mechanics simulations Lai-wai Tan, Vincent H. Chu* Department of Civil Engineering and Applied Mechanics, McGill University Montreal, Quebec, Canada * E-mail: [email protected] ABSTRACT: The Lagrangian block hydrodynamics is formulated based on the block advection of fluid. By enforcing the mass and momentum conservations on the Lagrangian mesh, the numerical oscillation problem encountered in the classical Eulerian computational methods is circumvented. A large number of the previously computationally difficult problems in environmental fluid mechanics are successfully simulated using the method. Examples of these simulations are described in this paper. KEY WORDS: Lagrangian description; computational hydrodynamics; block advection; shallow-water waves; turbulence simulations.

1 INTRODUCTION Two methods have been employed to describe the fluid motion. The Lagrangian description, although is the method for small deformation in solid bodies, is not usually used to analyze the motion in fluid. The common method to track the large and continuous deformation of the fluid is the Eulerian description. Most of the classical computational fluid dynamic is formulated based on the Eulerian method. Central to the Eulerian formulation is the evaluation of the fluxes. In the classical finite volume formulation, the fluxes are evaluated on the faces of the control volume using the values at the nodal points on the Eulerian mesh. Unfortunately, the error of this evaluation produces spurious numerical oscillations. In shallow-water wave and turbulence computations, the sudden changes in the depth and velocity in region across the shock waves and at the advancing and recessing fronts between the wet-and-dry interfaces are sensitive to the spurious oscillations. The negative water depth and/or the negative energy initiated by the oscillations usually result in the collapse of the numerical

computation. One method to avoid the evaluation of the fluxes is to formulate the computational problem using the Lagrangian block hydrodynamics. Chu and Altai[1,2] used blocks of fluid and Lagrangian advection of the blocks to conduct turbulence simulations in a streamfunction and vorticity formulation. Tan and Chu[3−7] extended the Lagrangian block method for onedimensional (1D) simulation of the waves in shallow waters using the primitive variable formulation. The extension of the primitive variable formulation to twodimensional (2D) problem is developed in this paper for environmental fluid mechanics application. 2 LAGRANGIAN BLOCK HYDRODYNAMICS The Lagrangian block hydrodynamics (LBH) is based on block advection. The blocks are material elements attached to the fluid on a Lagrangian mesh. The requirements for the conservations of mass and momentum are enforced on the Lagrangian mesh. The mass and momentum in the fluid are transferred by the blocks as the Lagrangian mesh deforms with the fluid. The LBH does not evaluate the fluxes. Therefore, the computational instabilities associated with the fluxes are circumvented. To overcome the intrinsic difficulty associated with large deformation in fluid, the blocks on the Lagrangian mesh are reconstructed to coincide with the Eulerian mesh in every time increment. The method of the LBH is to be distinguished with other Lagrangian method such as the particle-in-cell (PIC) methods of Harlow[8] and the smoothed particle hydrodynamics (SPH) methods of Monaghan[9]. The PIC and SPH methods use artificial particles. The LBH method uses the real fluids in the form of blocks

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

628

as the computational elements. The kernel function used in PIC and SPH to calculate the force between the artificial particles is not required. Eulerian grid u(i,j) h(i,j)

ΔxL t+Δt

Eulerian grid

Δy

Δy

(Fx )i , j = − g ΔyL

t h(i,j) u(i,j)

Δx

Eulerian grid

t+Δt ΔyL Δy v(i,j) Δx

Eulerian grid

(c) x-momentum block

u(i,j) t Δy

ΔxL t+Δt ΔyL

v(i,j) Δx

(d) y-momentum block

Fig. 1 (a) Staggered Eulerian mesh, (b) Volume block, (c)-(d) x- and y-momentum blocks at the beginning and the end of a Lagrangian advection time step

The present LBH is formulated using a staggered system of blocks for shallow-water wave simulation. Figure 1(a) shows the staggered grid and the relative locations of the volume block (hLi,jΔxLi,jΔyLi,j) to the momentum blocks [(ρuLi,jhLi,jΔxLi,j), (ρvLi,jhLi,jΔyLi,j)] on the grid. The x- and the y-components of the velocity are defined at a distance of ½ Δx to the west and ½ Δy to the south of the depth node, respectively. Figure 1(b) shows the advection of the volume block and Figs. 1(c) and 1(d) show the momentum blocks. A block is defined by its water depth hLi,j and the block widths xLi+1,j − xLi,j = ΔxLi,j and yLi,j+1 − yLi,j = ΔyLi,j. At the beginning of the Lagrangian advection, at time t, the blocks fit the Eulerian mesh, that is xLi,j = xi,j and yLi,j = yi,j. At the end of the advection step, at time t + Δt, ΔxLi,jΔyLi,jhLi,j = Δxi,jΔyi,jhi,j for the volume conservation. The edge positions of the blocks xLi,j and yLi,j at time t + Δt are calculated by the momentum equations: Du L i , j = (Fx )i , j , Dt

Δx hi , j − hi , j −1

Δy

− g (S ox − S fx )i , j

(2)

− g (S oy − S fy )i , j

(3)

(b) volume block

ΔxL t u(i,j)

= −g

hi , j − hi −1, j

where So = bottom elevation, g = gravity, Sf = cf ui,j |ui,j|/(2gh) = friction slope, and cf = coefficient of friction.

Δx

(a) staggered grid

(F )

y i,j

v(i,j) v(i,j)

functions of time only. For the shallow water waves, the forces on the blocks are calculated by assuming hydrostatic pressure variation over the depth:

Dv L i , j = (Fy )i , j Dt

(1)

where uLi,j = x-velocity, vLi,j = y-velocity. In the Lagrangian reference frame, the positions xLi,j(t) and yLi,j(t), and the velocities uLi,j(t) and vLi,j(t) are

To prevent entanglement between the Lagrangian paths, the mass and the momentum in the blocks are re-constructed on the Lagrangian mesh to match the Eulerian mesh at each computational time step[1,2]. A solution is possible when the Courant number Co = min[(u2+v2)½Δt/(Δx2+Δy2)½, (gh)½Δt/(Δx2+Δy2)½] is less than unity. An extensive series of grid refinement studies have been carried out by Tan and Chu[3−7] to show the convergent of the LBH simulations toward exact solutions for the dam-break flood waves by Ritter[10], Stoker[11], Hogg[12], Ancey et al.[13], and for the runup and overtopping of collapsing bores by Shen and Meyer[14] and Peregrine and Williams[15]. The computational accuracy and stability of the LBH method is further demonstrated by the 2D simulation examples to be presented in this paper. 3 SIMPLE ADVECTION – EXACT SOLUTION OF LBH For an observer moving with a uniform velocity (u, v) in a Lagrangian reference frame, the concentration cL is constant according to the advection equation Dc L ∂c ∂c ∂c = +u +v = 0 Dt ∂t ∂x ∂y

(4)

Figures 2(a) and (b) show the result obtained using LBH method and Figs. 2(c) and (d) are the result obtained using the QUICK (quadratic upstream interpolation for convective kinetics) finite volume method. The numerical solution by the LBH method is exact. However, the solution obtained by QUICK is contaminated by the spurious numerical oscillations produced by the truncation error in the estimate of the fluxes. Various numerical strategies have been developed in the past to suppress the spurious numerical oscillations[16−19]. Most of these were developed for the simple advection equation. Application of these strategies to the full set of continuity and momentum equations does not generally lead to computational stability[20].

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

629

(a) t = 0.1 s

h ho

ho = 10 m

1 shock front

0.5

(d)

hs ho = 0.396 0 −2

hd = 1 m (a) LBH, t = 0

(c) QUICK, t = 0

A

−1

(e)

0 −2

−1

A

(uStoker − us ) uStoker

(b) LBH, t = 5.64 s

0 1 x (t gho )

The following sections show examples of LBH computations to capture oblique dam-break shock waves and to track the wetting-and-drying interface of the back and forth sloshing of water in a parabolic bowl. 4 SHOCK CAPTURE The simulation of the oblique dam-break wave is an example to show the shock-capture ability of LBH. The oblique shock waves are produced in a 100 m × 100 m square basin as shown in Fig. 3. The water in the basin is divided by a dam along the diagonal of the basin with an initial upstream depth ho = 10 m and a downstream depth hd = 1 m. The removal of the dam along the diagonal causes the water to advance and form oblique shock waves. Figures 3(a), 3(b), and 3(c) show the depth contours obtained from the LBH simulations at time t = 0.1 s, 2.5 s and 4.0 s, respectively. Figures 3(d) and 3(e) show the water depth and velocity profiles computed by the LBH method at time t = 2.5 s. Exact solution of Stoker[11] is denoted by the circles for comparison. The LBH captured accurately the oblique shock wave without using any additional shock capturing scheme. Figure 3(f) shows the convergence towards the exact solution hs/ho = 0.396 as the block size is refined following the first-order convergent law of Celik et al.[21].

2 (f)

0.01 0.001

(d) QUICK, t = 5.64 s

Fig. 2 (a) Initial concentration block defined in LBH, (b) concentration at t = 5.64 s by LBH, (c) initial concentration block defined in QUICK, (d) concentration at t = 5.64 s by QUICK

2

1 u gho 0.5

(b) t = 2.5 s

(c) t = 4 s

0 1 x (t gho )

0.0001 0.001

m= 1 0.01 Δx (t gho )

0.1

Fig. 3 (a-c) Simulated water depth contours of oblique dambreak waves on wet downstream depth using mesh size Δx = 0.04 m, (d-e) cross sectional profile along A-A at time t = 2.5 s and (f) the convergence of frontal velocity based on five block size of Δx = 1 m, 0.5 m, 0.2 m, 0.1 m and 0.04 m

5 TRACKING WET-AND-DRY INTERFACE The capability of the LBH to track the wetting and drying interface is demonstrated by the computation of the water sloshing in a parabolic bowl as shown in Fig. 4. This problem has an exact solution by Thacker[22]. Initiated by a parabolic mound of water, the water moves up and down inside the bowl under the influence of gravity. The wetting and drying on the surface of the bowl is a very challenging numerical problem. If Eulerian numerical methods were used, the spurious numerical oscillations would produce negative water depth at the wave front leading to computation breakdown. However, the LBH computation has accurately tracked the infinite cycles of advances and recesses of the water on the surface of the parabolic bowl with absolute computational stability. Figure 5 shows the potential and kinetic energies of the water in the parabolic bowl obtained by the LBH computation and the rate of energy loss over a period of 20 cycles of water sloshing. The total (kinetic plus potential) energy E is supposed to be constant but is reduced with time due to computation error. Therefore, the energy loss is a measure of the computation error. The analysis of this error against the block size shows that the convergence of the LBH simulation towards the exact

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

630

solution is following the first order (p ≈ 1.0) law of convergence of Celik et al.[21].

h (m)

h (m)

h (m)

h (m)

2.5 2 1.5 1 0.5 0 2.5 2 1.5 1 0.5 0 2.5 2 1.5 1 0.5 0 2.5 2 1.5 1 0.5 0 −4 −3 −2 −1 0 1 x (km)

2

3

3 2 u 1 (m/s) 0 −1 −2 −3 3 2 u 1 (m/s) 0 −1 −2 −3 3 2 u 1 (m/s) 0 −1 −2 −3 3 2 u 1 (m/s) 0 −1 −2 −3 4 −4 −3 −2 −1 0 1 x (km)

Energy (MJoule/(kg/m3 ))

10

Initial Energy Eo

2

3

4

ΔE

5

Kinetic Energy

0

Potential Energy 20000

5000

10000 Time (s)

15000

To demonstrate the existence of macro resistance, a series of LBH simulations have been conducted for the conceptual river-flow channel. Figures 7 and 8 show the dilation and vorticity profiles obtained by two of the LBH simulations. Besides the coherent turbulence, the breaking gravity waves are observed to play a significant role. (b)

Fig. 4 Water depth h and velocity u in the parabolic bowl of 2 km-radius as computed by the block advection using Δx = Δy = 5 m. The computed profiles on the plane of symmetry are compared with the exact solution of Thacker[22] denoted by the circle symbol (a)

also for the morphological changes to the channel.

Fig. 5 Kinetic plus potential energy of the water in the parabolic bowl computed by LBH using Δx = Δy = 5 m

6 CONCEPTUAL RIVER MODELS The advance of flood water into dry basins has been simulated by the LBH using a conceptual river-flow model. The simulation shown in Fig. 6 starts with initially dry basins on the side of a supercritical flow in the main channel. Despite the high speed of the supercritical flow, the LBH has been able to simulate the advance of the water into the dry basins with absolute computational stability (beside the usual Courant-Friedrichs-Lewy requirement). The resistance to the flow in the river model is dependent on the bottom friction and also on the macro resistance produced by the channel features such as the flood plains, spurs dikes and bends along the channel. The real river channels however are variable in space and time due to erosion and deposition of sediments along the channel[23]. A reliable computational model is necessary not only for the simulation of the flow but

Fig. 6 Advance of flood waters onto the dry basins at time t = 0.2 s to 0.8 s; (a) water depth h and (b) vorticity ω profiles of the flow of initial Froude number Fr = 1.30

Θ

(cms−1) 15

>1

0.5

0

ω

(s−1)

> 50 25

0

−25 <−50

(a) t = 0.6 s

y (cm) 7.5

0 y (cm) 7.5

(b) t = 0.8 s

0 y (cm) 7.5

(c) t = 1.0 s

0 y (cm) 7.5

0

−0.5 <−1

(d) t = 1.2 s

7.5

15 x (cm)

22.5

0

7.5

15 x (cm)

22.5

30

Fig. 7 Dilation Θ profile (left) and vorticity ω profile (right) for the river flow of water depth ho = 6 cm, slope So = 0.01, bottom friction coefficient cf = 0, and Froude number Fr = 1.06

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China Θ

(cms−1) 15

>2

1

0

−1 <−2

ω

(s−1)

> 100 50

0

−50 <−100

(a) t = 0.6 s

y (cm) 7.5 0 y (cm) 7.5

(b) t = 0.8 s

631

The steady state velocity u and Froude number Fr profiles at time t = 300 s are given in Fig. 11 for slopes of (a) So = 0.005, (b) So = 0.01, and (c) So = 0.02. 30

u 5 4 3 2 1 0 −1 −2 −3 −4 −5 (m/s)

Fr 2

1.5

1

0.5

0

(a) So = 0.005

y 15 (m)

0 y (cm) 7.5

(c) t = 1.0 s

0

(b) So = 0.01

y 15 (m)

0 y (cm) 7.5

(d) t = 1.2 s

0

(c) So = 0.02

y 15 (m)

7.5

0

15 x (cm)

22.5

0

7.5

15 x (cm)

22.5

30

0 30

45

60 x (m)

75

90 30

45

60 x (m)

75

90

Fig. 8 Dilation Θ profile (left) and vorticity ω profile (right) for river flow of water depth ho = 6 cm, slope So = 0.1, coefficient of friction cf = 0, and Froude number Fr = 2.68

Fig. 10 Steady-state velocity u profile (left) and Froude number Fr profile (right) at time t = 300 s for bed slopes (a) So = 0.005, (b) So = 0.01, and (c) So = 0.02

7 MEANDERING RIVER MODEL

8 SUPERCRITICAL JET

The meandering river as shown in Figs. 9 and 10 is another series of LBH simulations conducted for the macro resistance. The river has a 10 m wide sinusoidal meander in a flood plain of 30 m width. Initial water depth in the river is ho = 1 m. The elevation of flood plain is 2 m higher than the meander’s bed. Figure 9 shows the velocity u and the Froude number Fr for the flow on a channel of slope So = 0.02 at time t = 40 s, 46 s and 52 s. Fig. 10 shows the steady state flow for three channel slopes So = 0.005, 0.01 and 0.02.

The supercritical and subcritical jets results shown in Fig. 11 are examples to show the application of the LBH to turbulence simulation and the method’s application to the pollutant and sediment transport processes. Shear instabilities in the supercritical shear flow is affected by the wave radiation and is dependent on the convective Froude number[24−26].

30

u 5 4 3 2 1 0 −1 −2 −3 −4 −5 (m/s)

Fr 2

1.5

1

0.5

43 t= 2s

t = 2s

y (m)

ω (s−1) 1

37.5

0.8 0.6

0

43

(a) t = 40 s

t = 8s

y (m)

0.4

t= 8s

0.2 37.5

0

y 15 (m)

−0.2 −0.4

43

0

(b) t = 46 s

32

(c) t = 52 s

45

−0.6 −0.8

59

−1

65

71 x (m)

77

(a) ho = 0.305 m, h d = 0.1 m, Fro = 1.42

y 15 (m) 0 30

t = 16 s

37.5

y 15 (m) 0

t = 16 s

y (m)

60 x (m)

75

90 30

45

60 x (m)

75

90

Fig. 9 Velocity u profile (left) and Froude number Fr profile (right) at the instance when the meander overflows it’s bank. The slope So = 0.02, and bottom friction coefficient cf = 0

83

65

71 x (m)

77

83

(b) ho = 0.305 m, h d = 0.255 m, Fr o = 0.57

Fig. 11 Vorticity profiles are obtained by LBH to contrast the difference in the turbulence flow of (a) a super-critical jet and (b) a subcritical jet. The supercritical flow process in (a) are characterized by the formation of shocklets and the radiation of waves from the shocklets. The subcritical flow process in (b) are simply the formation of the eddies

632

9th International Conference on Hydrodynamics October 11-15, 2010 Shanghai, China

9 CONCLUSIONS The extension of the Lagrangian block hydrodynamics to primitive variable formulation has made the simulations of a number of previously difficult environmental fluid mechanics problems possible. These include the flood wave propagation on dry bed and the high speed supercritical turbulent shear flow. Sudden change in depth and velocity across the shock waves and the advancement of water on dry bed are computed by the LBH method with absolute computational stability. REFERENCES [1] Chu V H, Altai W. Simulation of shallow transverse shear flow by generalized second moment method. J Hydraul Res, 2001, 39(6): 575−582. [2] Chu V H, Altai W. Simulation of turbulence and gravity interfaces by Lagrangian block method. Computational Fluid Dynamics 2002, New York: Springer-Verlag, 2002, 299−304. [3] Tan L W, Chu V H. Simulation of wave fronts on dry beds using Lagrangian blocks. Engrg Comput Mech, 2009a, 162(EM2): 57−66. [4] Tan L W, Chu V H. Lauber and Hager’s dam-break wave data for numerical model validation. J Hydraul Res, 2009b, 47(4): 524−528. [5] Tan L W, Chu V H. Wave runup simulations using Lagrangian blocks on Eulerian mesh. J Hydro-environ Res, 2010a, 3(4): 193−200. [6] Tan L W, Chu V H. Wet-and-dry interface on steep slopes simulations using Lagrangian blocks. Proc 6th Int Symp Env Hydraul, 2010b, Athens, Greece. [7] Tan L W, Chu V H. Dam-break flood simulations using 2D Lagrangian blocks on Eulerian mesh method. Proc 6th Int Symp Env Hydraul, 2010c, Athens, Greece. [8] Harlow F H. The particle-in-cell computing method for fluid dynamics. Methods Comput Phys, New York: Academic Press, 1964, 3: 319−343. [9] Monaghan J J. Smoothed particle hydrodynamics. Reports Progress Phys, 2005, 68(8): 1703−1759. [10] Ritter A. The propagation of water waves. Zeitschrift des Vereines Deutsher Ingenieure, 1892, 36(33): 947−954. [11] Stoker J J. Water Waves: The Mathematical Theory with Applications. Wiley Interscience, 1957.

[12] Hogg A J. Lock-release gravity currents and dam-break flows. J Fluid Mech, 2006, 569(1): 61−87. [13] Ancey C, Iverson R M, Rentschler M, et al. An exact solution for ideal dam-break floods on steep slopes. Water Resources Res, 2008, 44: W01430. [14] Shen M C, Meyer R E. Climb of a bore on a beach: part 3, run-up. J Fluid Mech, 1963, 16(1): 113−125. [15] Peregrine D H, Williams S M. Swash overtopping a truncated plane beach. J Fluid Mech, 2001, 440: 391−399. [16] van Leer B. Towards the ultimate conservative difference scheme, II, Monotonicity and conservation combined in a second-order scheme. J Comput Phys, 1974, 14: 361−370. [17] Gaskell P H, Lau A K C. Curvature compensated convective transport: SMART, a new boundedness preserving transport algorithm. Int. J Num Meth in Fluids, 1988, 8: 617−641. [18] Leonard B P, Mokhtari S. Beyond first-order upwinding: the ULTRA-SHARP alternative for non-oscillatory steady-state simulation of convection. Int J Num Meth in Fluids, 1990, 30: 729−766. [19] Leonard B P, Drummond J E. Why you should not use ‘HYBRID’, ‘POWER-LAW’ or related exponential schemes for convective modelling – there are much better alternatives. Int J Num Meth in Fluids, 1995, 20: 421−442. [20] Pinilla C E, Bouhairie S, Tan L-W, et al. Minimal intervention to simulations of shallow-water equations. J Hydro-env Res, 2010, 3(4): 201−207. [21] Celik I B, Ghia U, Roache P J, et al. Procedure for estimation and reporting of uncertainty due to discretization in CFD applications. J Fluids Engrg, 2008, 130: 1−4. [22] Thacker C W. Some exact solutions to the nonlinear shallow-water wave equations. J Fluid Mech, 1981, 107: 499−508. [23] Webb R H, Schmidt J C, Marzolf G R, et al. The controlled flood in Grand Canyon. American Geophy. Union, 1999, Monograph 110. [24] Pinilla C, Chu V H. Waves and bed-friction effect on stability of transverse shear flow in shallow waters. J Coastal Res, 2008, 52: 207−214. [25] Pinilla C, Chu V H. The role of wave radiation on instability of coastal current. Proc. 33rd Congress IAHR, 2009a, Vancouver, Canada. [26] Pinilla C, Chu V H. Wave radiation and shear instability in rotating stratified flow. Proc 28th Int Conf on Ocean, Offshore and Arctic Engrg, 2009b, Honolulu, Hawaii.