Proceedings of the Seventh International Conference on Hydroscience and Engineering, Philadelphia, PA, September 2006.  http://hdl.handle.net/1860/732The 2D shallow water equations (SWE) have been widely used in hydraulic science and
engineering. Numerical discretization based on midpoint rule (e.g. leapfrog method) is desirable for
2nd order temporal accuracy. However, the explicit 3-level Leapfrog approach must be applied with
appropriate strategies to decouple the known instability mode (Peyret and Taylor, 1983). In this
paper, we propose a “Time-Centered, Split-Implicit” (TCSI) method based on a 2-level midpoint
rule. This new method computationally linearizes the advection term such that individual solution
steps are implicit and linear, but their combination is discretely equivalent to nonlinear midpoint rule
discretization. To verify the capability of the new method, the 2D SWE are used to develop and test
the TCSI approach. Test cases of simulating a standing wave in a closed basin are performed to
examine a variety of model characteristics

Fu, Shipeng

Hodges, Ben R.

Drexel Libraries E-Repository and Archives

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
  
 
Proceedings of the 7th International Conference on HydroScience and Engineering 
Philadelphia, USA September 10-13, 2006 (ICHE 2006) 
ISBN: 0977447405 
 
Drexel University 
College of Engineering 
 
 
 
 
Drexel E-Repository and Archive (iDEA) 
http://idea.library.drexel.edu/   
 
 
Drexel University Libraries 
www.library.drexel.edu 
 
 
 
The following item is made available as a courtesy to scholars by the author(s) and Drexel University Library and may 
contain materials and content, including computer code and tags, artwork, text, graphics, images, and illustrations 
(Material) which may be protected by copyright law. Unless otherwise noted, the Material is made available for non 
profit and educational purposes, such as research, teaching and private study. For these limited purposes, you may 
reproduce (print, download or make copies) the Material without prior permission. All copies must include any 
copyright notice originally included with the Material. You must seek permission from the authors or copyright 
owners for all uses that are not allowed by fair use and other provisions of the U.S. Copyright Law. The 
responsibility for making an independent legal assessment and securing any necessary permission rests with persons 
desiring to reproduce or use the Material. 
 
 
Please direct questions to archives@drexel.edu 
 
The 7th Int.  Conf. on Hydroscience and Engineering (ICHE-2006),Sep 10 –Sep 13,Philadelphia,USA 1 
 
 
 
 
 
 
INITIAL TESTS OF A NEW COMUPATIONAL LINEARIZATION METHOD  
 
 
Shipeng Fu1 and Ben Hodges2 
 
 
ABSTRACT 
 
The 2D shallow water equations (SWE) have been widely used in hydraulic science and 
engineering. Numerical discretization based on midpoint rule (e.g. leapfrog method) is desirable for 
2nd order temporal accuracy. However, the explicit 3-level Leapfrog approach must be applied with 
appropriate strategies to decouple the known instability mode (Peyret and Taylor, 1983). In this 
paper, we propose a “Time-Centered, Split-Implicit” (TCSI) method based on a 2-level midpoint 
rule. This new method computationally linearizes the advection term such that individual solution 
steps are implicit and linear, but their combination is discretely equivalent to nonlinear midpoint rule 
discretization. To verify the capability of the new method, the 2D SWE are used to develop and test 
the TCSI approach. Test cases of simulating a standing wave in a closed basin are performed to 
examine a variety of model characteristics.  
 
1. INTRODUCTION 
The Leapfrog method, based on midpoint rule, is popular in meteorology and oceanography (Fujima 
and Shigemura 2000). The typical leap-frog approach is a three-level method, discretized as 
 
 ( )n 1 n 1 n n n 3f , t O( t )+ −φ = φ + φ φ φ Δ + Δ  (1) 
 
where f is a function of a nonlinear advection term ( n nφ φ ) and a diffusion term ( nφ ).  The midpoint 
rule is a desirable discretization as it has second-order temporal accuracy. However, the explicit 
Leapfrog approach must be applied with appropriate strategies to decouple the known instability 
mode (Peyret and Taylor, 1983).  We would prefer a two-level method based on the midpoint rule 
that can be written as  
 
 ( )n 1 n n 1/ 2 n 1/ 2 n 1/ 2 3f , t O( t )+ + + +φ = φ + φ φ φ Δ + Δ  (2) 
 
which has the desirable property that the time n+1/2 information is retained only within the 
computation of the n to n+1 time step, and does not directly appear in the discrete equation for 
n+3/2 time step.  If the n+1/2 values in eq. 2 are further discretized by the linear approximation 
                                                 
1 Graduate Student, Department of Civil, Architectural, and Environmental Engineering, University of Texas at Austin, 
Austin, TX 78712, USA (shipengfu@mail.utexas.edu) 
2 Associate Professor, Department of Civil, Architectural, and Environmental Engineering, University of Texas at 
Austin, Austin, TX 78712, USA (hodges@mail.utexas.edu) 
  2 
( )n 1/ 2 n n 1 / 2+ +φ = φ + φ , the resulting equation is both implicit and nonlinear, which is generally not a 
good starting point for a numerical method.  However, below we introduce a new splitting method 
based on exactly this form.  This new method computationally linearizes the advection term such 
that individual solution steps are implicit and linear, but their combination is discretely equivalent to 
nonlinear eq. 1.  We call this a “Time-Centered, Split-Implicit” (TCSI) method. 
  
1.1 Computational Linearization 
To solve a nonlinear equation system, computational linearization techniques are typically the most 
effective approach. One example is the local time-linearization method (e.g. Lomax, Pulliam and 
Zingg, 1996).  In local time-linearization, an implicit nonlinear term (e.g. n 1 n 1+ +φ φ ), is approximated 
with information from the known and unknown time steps (e.g. n n 1+φ φ ) by using a Taylor-series 
expansion.  The resulting equation retains the nonlinear term, but time-lagging one variable provides 
a linear equation set for time-integration. However, use of n n 1+φ φ  to represent n 1 n 1+ +φ φ  provides 
first-order temporal accuracy.  The new TCSI approach is similar to local-time linearization, but is 
formally second-order accurate in time by using a two-step approach. 
 The TCSI method applies the approximation 
 
 ( )n n 1n 1/ 2 2O t2
+
+ φ + φφ = + Δ  (3) 
 
Substituting eq. 3 into eq. 2 for one of the advective terms provides: 
 
 ( )n n 1n 1 n 2 n 1/ 2 n 1/ 2 3f O t , t O( t )2
+
+ + +⎧ ⎫⎡ ⎤φ + φ⎪ ⎪φ = φ + + Δ φ φ Δ + Δ⎨ ⎬⎢ ⎥⎪ ⎪⎣ ⎦⎩ ⎭
 (4) 
 
which can be written as: 
 
 ( )n 1 n n n 1/ 2 n 1 n 1/ 2 n 1/ 2 3tf , , O( t )2+ + + + + Δφ = φ + φ φ φ φ φ + Δ  (5) 
 
Thus the nonlinear advection term is discretely linearized without losing the original accuracy. 
 
1.2 Computational Splitting 
Observing that eq. 5 has two different nonlinear terms, a computational splitting approach can be 
applied: 
 
 ( )* n n n 1/ 2 n 1/ 2 tf , 2+ + Δφ = φ + φ φ φ  (6) 
 ( )n 1 * n 1 n 1/ 2 n 1/ 2 3tf , O( t )2+ + + + Δφ = φ + φ φ φ + Δ  (7) 
 
where *φ  is an intermediate solution and substituting eq. 6 into eq. 7 returns eq. 5. To create a 
solvable system, it can be shown that  
 
  3 
 
2
n 1/ 2 * tO
2
+ Δ⎛ ⎞φ = φ + ⎜ ⎟⎝ ⎠  (8) 
 
so that we are free to substitute *φ  for n 1/ 2+φ  on the RHS of eq. 6 and eq. 7  without affecting the 
overall accuracy.  With this substitution, the original split system becomes: 
 
 ( )* n n * * 3tf , O( t )2Δφ = φ + φ φ φ + Δ  (9) 
 ( )n 1 * n 1 * * 3tf , O( t )2+ + Δφ = φ + φ φ φ + Δ  (10) 
 
Eq. 9 and eq. 10 are two implicit linear steps of a general time-centered split-implicit (TCSI) method 
that is formally second-order temporally accurate for nonlinear evolution equations.  Solution with 
second-order spatial differences will result in each step being a simple tridiagonal inversion. 
 
2. IMPLEMENTATION AND TESTING FOR  2D SHALLOW WATER EQUATIONS  
2.1 Discretized Equation 
We use the 2D shallow water equation to develop and test the TCSI approach. The 2D shallow 
water equations can be written as: 
 
 
2 2
2 2
U U U U UU V g
t x y x x y
⎛ ⎞∂ ∂ ∂ ∂ζ ∂ ∂+ + = − + ν +⎜ ⎟∂ ∂ ∂ ∂ ∂ ∂⎝ ⎠
 (11) 
 
2 2
2 2
V V V V VU V g
t x y y x y
⎛ ⎞∂ ∂ ∂ ∂ζ ∂ ∂+ + = − + ν +⎜ ⎟∂ ∂ ∂ ∂ ∂ ∂⎝ ⎠
 (12) 
 H HU HV 0
t x y
∂ ∂ ∂+ + =∂ ∂ ∂  (13) 
 bH Zζ = +  (14) 
 
where U  and V  are depth-averaged velocities in x  and y  directions; g is the gravitational 
acceleration; H  is the water depth and ζ  is surface elevation; bZ is the bottom elevation.  The 
above statement of the 2D SWE neglects the turbulence closure and bottom stress as a test case to 
our numerical method,   
 Applying the midpoint rule between time ‘n’ and n+1’ along with central-differencing of 
spatial derivatives on an Arakawa C grid (Arakawa and Lamb, 1977), eqs. 11, 12 and 13 can be 
discretized as a set of coupled nonlinear equations: 
 ( ) ( )
( )
n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2n 1 n
i, j 3/ 2 i, j 1/ 2 i 1, j 1/ 2 i 1, j 1/ 2i, j 1/ 2 i, j 1/ 2 n 1/ 2 n 1/ 2
i, j 1/ 2 i, j 1/ 2
n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2
i, j 1/ 2 i, j 1/ 2 i, j 3/ 2 i 1, j 1/ 2 i, j 1/ 22 2
U U U UU U
U V
t 2 x 2 y
U 2U U U 2U U
x y
+ + + ++
+ − + + − ++ + + +
+ +
+ + + + +
− + + − + +
− −− = − −Δ Δ Δ
ν ν+ − + + − +Δ Δ ( )
( )
n 1/ 2
i 1, j 1/ 2
n 1/ 2 n 1/ 2
i, j 1 i, j
g
x
+
+ +
+ +
+− ζ − ζΔ
 (15) 
 
  4 
( ) ( )
( )
n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2n 1 n
i 1/ 2, j 1 i 1/ 2, j 1 i 3/ 2, j i 1/ 2, ji 1/ 2, j i 1/ 2, j n 1/ 2 n 1/ 2
i 1/ 2, j i 1/ 2, j
n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2
i 1/ 2, j 1 i 1/ 2, j i 1/ 2, j 1 i 1/ 2, j i 1/ 2, j2 2
V V V VV V
U V
t 2 x 2 y
v vV 2V V V 2V
x y
+ + + ++
+ + + − + −+ + + +
+ +
+ + + + +
+ − + + + − +
− −− = − −Δ Δ Δ
+ − + + −Δ Δ ( )
( )
n 1/ 2
i 1/ 2, j
n 1/ 2 n 1/ 2
i 1, j i, j
V
g
y
+
+
+ +
+
+
− ζ −ζΔ
 (16) 
 
 ( ) ( )n 1 ni, j i, j n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2 n 1/ 2i, j 1/ 2 i, j 1/ 2 i, j 1/ 2 i, j 1/ 2 i 1/ 2, j i 1/ 2, j i 1/ 2, j i 1/ 2, jH H 1 1H U H U H V H Vt x y
+
+ + + + + + + +
+ + − − + + − −
− = − − − −Δ Δ Δ  (17) 
 
 
2.2 Discussion and Results of Simulating a Standing Wave in a Closed Basin 
As a test of the TCSI method, a free surface standing wave in a rectangular basin has been simulated 
for the basin shown in Figure 1. An equally-spaced 20 5×  mesh is used.  It is recognized that the 
model mesh and initial conditions provide essentially a 1D exercise that does not fully test the x-y 
decoupling of the method; however, this mesh does allow testing of how the method treats the 
nonlinear u u / dx∂  term. Free-slip boundary conditions are enforced on all walls and the bottom.  
L
H
2a
x
yz
 
Figure 1  A standing wave in a rectangular basin 
  
The period of an inviscid free surface wave is (Dean and Dalrymple 1998): 
 
 ( ) 12T tanh kH
g
−πλ=  (18) 
 
where ‘T’ is the wave period; ‘λ ’ is the wave length; ‘k’ is the wave number and ‘H’ is the water 
depth. According to linear wave theory, the wave function is a sinusoid for small amplitude waves. 
The surface height (h) above the still water level is: 
 
 ( ) ( )h(x, t) a sin kx sin t= σ  (19) 
 
  5 
where the frequency (σ) is 2 / Tπ ; ‘a’ is the free surface wave amplitude. For a viscous unconfined 
wave in deep water, the evolution of the wave amplitude can be approximated as (Lamb 1932), 
 
 
22 k ta(t) a(0)e− ν=  (20) 
 
where ‘ν ’ is kinematic viscosity and ‘ k ’ is the wave number.  However, wave damping based on 
eq. 20 is unlikely to be well-represented in a 2D depth-averaged model as vertical velocity and 
vertical velocity shears are part of the closure term. The Reynolds number is defined as: 
 
 LuRe = ν  (21) 
 
where L  is the length of the basin, ν  is the kinematic viscosity and u  is the characteristic Cartesian 
fluid velocity based on the wave amplitude (a) and the wave frequency (σ ) : 
 
 u a= σ  (22) 
 
We have run simulations for different cases to examine a variety of model characteristics. To 
find out how out model perform in different ratios of water depth to basin length (i.e. H / L ), three 
different cases are tested and compared in Figure 2. Furthermore, three simulations with different 
initial wave steepness (i.e. ratio of wave amplitude to basin length) are presented in Figure 3. All the 
six simulations are run in an inviscid flow to exclude viscous damping effects. Results in Figure 2 
and Figure3 illustrate the nondimensionalized water surface elevation (i.e. ( )H / aζ − ) at x x / 2= Δ , 
which is the grid cell center closest to the left wall. The important parameters for these six test cases 
are listed in Table 1.   
 
Table 1 Parameters of simulations of a standing wave 
 
case H/L a/L Re L /λ  
3a 0.05 5 .0x 10-4 ∞  0.5 
3b 0.1 5 .0x 10-4 ∞  0.5 
3c 0.2 5 .0x 10-4 ∞  0.5 
4a 0.05 5 .0x 10-4 ∞  0.5 
4b 0.05 2 .0x 10-3 ∞  0.5 
4c 0.05 5 .0x 10-3 ∞  0.5 
 
In Figure 2, three cases with different H / L  ratio are compared. As the water depth increases, 
the simulation wave period is shortened.  The wave dispersion behavior of this method requires 
further research. 
In Figure 3, we compare waves with different initial steepness. Wave steepening is inherently 
a nonlinear effect, which increases with increasing initial amplitude, thus we do not expect 
agreement between linear theory and. As the SWE governing equations retain nonlinearity but 
neglect non-hydrostatic dispersion, the wave steepening will be greater. Moreover, we keep the 
nonlinearity while the hydrostatic approximation is inherent in SWE. This two facts will introduce 
additional steepening (Hodges et al. 2006). In case 4a, a very small initial wave steepness is 
introduced and the simulation show little nonlinear effects. With the initial steepness increased in 
case 4b and 4c, the nonlinear effects are more and more obvious. 
 
  6 
0 1 2 3
-1
0
1
t/T
( ζ
-H
)/a
case a: H/L=0.05
0 1 2 3
-1
0
1
t/T
( ζ
-H
)/a
case b: H/L=0.1
0 1 2 3
-1
0
1
t/T
( ζ
-H
)/a
case c: H/L=0.2
 
 
Figure 2  Modeled water surface elevation at x=0 for different ratios of H/L with small amplitude 
wave (a/L = 5 x 10-4) and inviscid model.  analytical solution, simulation. 
 
 
 
  7 
0 1 2 3
-1
0 
1 
t/T
( ζ
-H
)/a
case a: a/L=0.05%
0 1 2 3
-1
0 
1 
t/T
( ζ
-H
)/a
case b: a/L=0.2%
0 1 2 3
-1
0 
1 
t/T
( ζ
-H
)/a
case c: a/L=0.5%
 
 
Figure 3 Modeled water surface elevation at x=0 for different ratios of a/L in shallow water (H/L = 
0.05) with inviscid model.  analytical solution, simulation. 
 
Figure 4 (below) displays results for a viscous simulation.  Walls and boundaries are free slip, 
so the only viscous terms are wave-induced (rather than wave-boundary induced). The Reynolds 
number is equal to 0.1. Our simulation result damped slower than the analytical solution, eq. 20. 
This result is not unexpected because that 2D SWE neglect the velocity gradients in the vertical 
(i.e. u / z 0∂ ∂ = ).  As a result, the viscous damping due to 2 2u / zν∂ ∂  is not reflected in our 
simulation. Thus, 2D SWE cannot correctly represent the damping of a standing wave without an 
additional modeling term to represent the viscous term which is partially lost in the depth averaging 
process.  
 
0 1 2 3
-1
0
1
t/T
( ζ
-H
)/a
Re=0.1,  H/L=0.05,  a/L=0.05%
 
 
Figure 4 Viscous simulation results.  analytical solution, simulation. 
 
  8 
 
 
REFERENCES 
Arakawa, A., and Lamb, V.R. (1977). “Computational design of the basic dynamical processes of 
the UCLA general circulation model”, Methods in Computational Physics, vol.17, Academic 
Press, New York, 173-265. 
Dean, R.G. and Dalrymple, R.A. (1998). Water Wave Mechanics for Engineering and Scientists, 
World Scientific Publishing Co. Pte. Ltd., New Jersey. 
Fujima, K. and Shigemura, T. (2000). “Determination of Grid Size for Leap-Frog Finite Difference 
Model to Simulate Tsunamis around A Conical Island”, Costal Engineering, vol42, No.2, pp. 
197-210. 
Hodges, B.R., Laval, B., and Wadzuk, B.M. (2006).  “Numerical error assessment and temporal 
horizon for internal waves in a hydrostatic model”, Ocean Modeling Vol. 13, No. 1, pp. 44-64. 
Lamb, H. (1932). Hydrodynamics. Cambridge University Press. 
Lomax, H., Pulliam, T.H. and Zingg, D.W. (1996) Fundamentals of Computational Fluid Dynamics. 
Springer, New York Inc. 
Peyret, R. and Taylor, T.D. (1983). Computational Method s for Fluid Flow. Springer-Verlag, New 
York Inc. 
 


Initial tests of a new comupational linearization method

https://idea.library.drexel.edu/islandora/object/idea%3A1436/datastream/OBJ/download/Initial_tests_of_a_new_comupational_linearization_method.pdf

Initial tests of a new comupational linearization method

Abstract

Similar works

Full text

Available Versions

Drexel Libraries E-Repository and Archives