i
CRWR Online Report 0804
A TimeCentered Split for Implicit Discretization of Unsteady
Advection Problems
by
Shipeng Fu, B.S.; M.S.
and
Ben R. Hodges, Ph.D.
The University of Texas at Austin
May 13, 2008
This document is available online via World Wide Web at
http://www.crwr.utexas.edu/online.html
ii
Copyright
by
Shipeng Fu
2008
iii
Acknowledgements
I would like to begin by expressing my gratitude to my advisor Dr. Ben Hodges. I
am thankful for his guidance and support, and I appreciate the freedom he has given me
to explore and pursue my own research interests. I would also like to thank Dr. Graham
Carey, Dr. Randall Charbeneau, Dr. Spyros Kinnas and Dr. Daene Mckinney for serving
on my committee.
The faculty members, facilities, and the world class program of EWRE have made
my research possible. I am grateful to Dr. Kinnas and Dr. Carey for their input in the
development of the numerical methods in this dissertation. I have gone to Dr. Liljestrand,
Dr. Lawler and Dr. Kinnas many times for help and advice, and they have always been
happy to help. I want to thank them for their encouragement and advice. I would like to
thank Dr. Charbeneau for giving me the opportunity to work with lab experiments (so I
could finally play with the real water!). I would like to thank Sharon Bernard, Marcy
Betak, and Susan SwansonCartwright at the CRWR Office for helping with all my
administrative paperwork and questions.
During my Ph.D study at UT, my friends in Austin have been extremely helpful
and supportive. I could not finish this without their friendship. I want to particularly
thank Paula Kulis, LiJung Chen and their wonderful husbands for feeding me meals and
listening to my whines. I also want to thank Paula for many discussions in fluid
mechanics and beyond. I want to thank Yihsiang Yu and Vimal Vinayan for their
insightful discussions on numerical methods and wave mechanics. I want to thank Becky
Teasley and John Allen, Rebekah and Nate Johnson, and Jeremy Seibert for graciously
reading and providing valuable feedbacks on my dissertation. I want to thank all my
friends in EWRE and CRWR. Our lunch topic at CRWR is always intellectually
stimulating, and sometimes wild. I also want to thank all the members in the BUSTED
(Best Underpaid Students Turning into Environmental Doctors) club. I really benefited
from the support from everyone. I will miss you girls. I want to thank my running buddy,
Susan De Long, for many interesting discussions about Yoga and philosophy.
I am thankful to my parents. I thank them for teaching me good values, and
encouraging me to think independently since I was very young. My diligent father and
curious mother have always been my role models and inspiration.
My deepest thanks go to my best friend and husband, Yunzhi Yang. I am really
fortunate to have such a wonderful life partner accompanying me for this journey. I want
to thank him for truly understanding me, encouraging me to strive for excellence, and
bringing life into perspective.
iv
Abstract
Environmental flows (e.g. river and atmospheric flows) governed by the shallow
water equations (SWE) are usually dominated by the advective mechanism over multiple
timescales. The combination of time dependency and nonlinear advection creates
difficulties in the numerical solution of the SWE. A fullyimplicit scheme is desirable
because a relatively large time step may be used in a simulation. However, nonlinearity in
a fully implicit method results in a system of nonlinear equations to be solved at each
time step. To address this difficulty, a new method for implicit solution of unsteady
nonlinear advection equations is developed in this research. This TimeCentered Split
(TCS) method uses a nested application of the midpoint rule to computationally decouple
advection terms in a temporally secondorder accurate timemarching discretization. The
method requires solution of only two sets of linear equations without an outer iteration,
and is theoretically applicable to quadraticallynonlinear coupled equations for any
number of variables.
To explore its characteristics, the TCS algorithm is first applied to one
dimensional problems and compared to the conventional nonlinear solution methods.
The temporal accuracy and practical stability of the method is confirmed using these 1D
examples. It is shown that TCS can computationally linearize unsteady nonlinear
advection problems without either 1) outer iteration or 2) calculation of the Jacobian. A
family of the TCS method is created in one general form by introducing weighting factors
to different terms. We prove both analytically and by examples that the value of the
weighting factors does not affect the order of accuracy of the scheme. In addition, the
TCS method can not only computationally linearize but also decouple an equation system
of coupled variables using special combinations of weighting factors. Hence, the TCS
method provides flexibilities and efficiency in applications.
v
Table of Contents
List of Tables ........................................................................................................ vii
List of Figures...................................................................................................... viii
Chapter 1 Introduction .............................................................................................1
1.1 Background...............................................................................................1
1.2 Challenge and Motivation of the Numerical model..................................3
1.3 Objective...................................................................................................4
1.4 Approach...................................................................................................4
1.5 Broad Application.....................................................................................5
Chapter 2 Review of Finite Difference Schemes for Solving the Unsteady Nonlinear
Advection in the Shallow Water Equations....................................................6
2.1 Explicit Method ........................................................................................6
2.2 SemiImplicit Method...............................................................................7
2.3 Implicit Method ........................................................................................7
2.4 Summary.................................................................................................11
Chapter 3 Theoretical Development of the TimeCentered Split (TCS) Method..12
3.1 The Theory of TimeCentered Splitting .................................................12
3.2 Derivation of the TCS method in a 1D AdvectionDiffusion Equation .16
3.3 TCS for Coupled Momentum and Scalar Transport...............................19
3.4 Summary.................................................................................................20
Chapter 4 Implementation of the TCS Method in 1D Problems ...........................22
4.1 Application of the TCS Method to the 1D Conservative Burgers? Equation
..............................................................................................................22
4.2 Applications of the TCS Method to the 1D Nonconservative Burgers?
Equation ...............................................................................................42
4.3 Application of the TCS Method to a 1D Nonlinear Ordinary Differential
Equation ...............................................................................................59
4.4 Summary.................................................................................................63
vi
Chapter 5 The TCS Family Method.......................................................................65
5.1 Derivation of the TCS Family Method ...................................................65
5.2 Application of the TCS Family Method to the 1D Nonconservative
Burgers? Equation ................................................................................68
5.3 Results and Discussion ...........................................................................70
5.4 Computational Decoupling .....................................................................78
5.5 Summary.................................................................................................81
Chapter 6 Application of the TCS Method to a 2D Depth Averaged Shallow Water
Equations (SWE) ..........................................................................................83
6.1. The 2D Depth Averaged SWE...............................................................83
6.2 The TCS Discretized SWE .....................................................................84
6.3 Decoupling the SWE...............................................................................87
6.4 Characteristics of the TCS Decoupled Equation System........................89
6.5 Numerical Tests ......................................................................................91
6.6 Summary...............................................................................................109
Chapter 7 Conclusions and Recommendations....................................................110
7.1 Summary of Discussion ........................................................................110
7.2 Conclusions...........................................................................................113
7.3 Recommendations for Future Work......................................................113
Appendix A Test of a Progressive Wave in an Open Boundary System.............115
A.1 Initial and Inlet Boundary Condition ...................................................115
A.2 Outlet Boundary Condition..................................................................116
A.3 Results and Discussion.........................................................................119
References............................................................................................................124
vii
List of Tables
Table 5.1 The mathematical meaning of the value of
1
? and
2
? ...................................... 70
Table 5.2 Weighting factors for computational decoupling the 2D Burgers? equation.... 80
Table 6.1 Weighting factors and the solution orders for decoupled Equations (6.11)
through (6.16). .......................................................................................................... 89
Table 6.2 Solution procedures of the three TCS discretizations....................................... 91
Table 6.3 Parameters of simulations of a 1D standing wave............................................ 93
Table 6.4 Parameters of simulations of the 2D standing wave......................................... 97
viii
List of Figures
Figure 1.1 Sulphur River (Texas) large woody debris at lowflow conditions .................. 2
Figure 4.1 Solution of 1D Burgers? equation evolving in time for t?{0, 0.1, 0.2, 0.3, 0.6,
1, 2, 3} using 0.05?= along with the initial and boundary conditions of Equations
(4.43) through (4.45)................................................................................................. 34
Figure 4.2 Absolute error time evolution for numerical solutions of 1D conservative
Burgers? equation for TCSF and TCSFD using { t? =0.01, x? =1/50, ?=0.05},
where
50
abs model i analytical i
i1
1
(t) u (x ,t) u (x ,t)
50
=
?= ?
?
and
analytical
u (x, t) is numerically
calculated from Equation (4.46) with K=30. ........................................................... 35
Figure 4.3
RMS
? vs. CFL number of various methods for solution of 1D conservative
Burgers? equation with three different x? , where 0.05?= and t0.3/? =?with
??{50, 30, 25, 20, 10, 8, 5, 4, 2}............................................................................ 37
Figure 4.4
2
L norm vs. CFL number of various methods for solution of 1D conservative
Burgers? equation with three different x? , where 0.05?= and t0.3/? =?with
??{50, 30, 25, 20, 10, 8, 5, 4, 2}............................................................................ 38
Figure 4.5 L
?
norm vs. CFL number of various methods for solution of 1D conservative
Burgers? equation with three different x? , where 0.05?= and t0.3/? =?with
??{50, 30, 25, 20, 10, 8, 5, 4, 2}............................................................................ 39
Figure 4.6
RMS
? vs. CFL number of TCSF, RK2 and RK4 methods for solution of 1D
conservative Burgers? equation, where 0.05?= x? =1/50, and t0.3/? =?with
??{500, 200, 100, 80, 50, 30, 25, 20, 10, 8, 5}...................................................... 40
Figure 4.7 Ideal operations per grid point for one time step using various nonlinear
solution methods for 1D conservative Burgers? equation. R is the number of outer
iterations taken by Newton method. It is assumed that the Picard method
convergences in R
2
outer iterations........................................................................... 42
Figure 4.8 Absolute error time evolution for numerical solutions of 1D nonconservative
Burgers? equation for the TCSF and TCSG methods using
{ }t 0.01, x 1/ 50, 0.05?= ?= ?= where
50
abs model i analytical i
i1
1
(t) u (x ,t) u (x ,t)
50
=
?= ?
?
and
analytical
u(x,t) is numerically calculated from Equation (4.46) with K=30........ 54
ix
Figure 4.9
RMS
? vs. CFL numbers of various methods for solution of 1D non
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50
and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}.................................................. 55
Figure 4.10
2
L norm vs. CFL numbers of various methods for solution of 1D non
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50
and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}.................................................. 56
Figure 4.11 L
?
norm vs. CFL numbers of various methods for solution of 1D non
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50
and t? =0.6/? with ? ?{50, 30, 25, 20, 10, 8, 5}.................................................. 57
Figure 4.12 Temporal accuracy of the TCS and RK methods for solution of 1D non
conservative Burgers? equation at different CFL numbers, where ?=0.05, x? =1/50
and t? =0.6/? with ? ?{500, 200, 100, 80, 50, 30, 25, 20, 10, 8, 5}.................... 58
Figure 4.13 Ideal operations per grid point for one time step using various nonlinear
solution methods for 1D Nonconservative Burgers equation. R is the number of
outer iterations taken by Newton method. It is assumed that the Picard method
convergences in
2
R outer iterations.......................................................................... 59
Figure 4.14 Time evolution of absolute errors for TCS methods applied to the ODE,
Equation (4.101), for t? =0.1, where
abs model analytical
(t) y (t) y (t)?= ? . .................... 62
Figure 4.15 RMS error from Equation (4.117) for discrete solutions of the ODE,
Equation (4.101), computed using a range of ? time steps. .................................... 63
Figure 5.1 Temporal accuracy of various combinations of
1
? and
2
? for solution of the
Burgers? equation at t=0.6, where x150/? = , 005.?= and t06./? =?
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 72
Figure 5.2 Temporal accuracy of various combinations of
2
? and
1
? for solution of the
Burgers? equation at t=0.6, where x150/? = , 005.?= and t06./? =?
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 73
Figure 5.3 Temporal accuracy of various combinations of
2
? and
1
? for solution of the
Burgers? equation at t=0.3, where x150/? = , 005.?= and t03? =?./
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 75
Figure 5.4 Temporal accuracy of various combinations of
2
? and
1
? for solution of the
Burgers? equation at t=1, where x150? = / , 005?= . and t1? =?/
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 76
x
Figure 5.5 Temporal accuracy of various combinations of
2
? and
1
? for solution of the
Burgers? equation at t=2, where x150/? = , 005.?= and t2? =?/
( { }5, 8, 10, 20, 25, 30, 50?? ).................................................................................. 77
Figure 6.1 Illustration of Arakawa C grid......................................................................... 85
Figure 6.2 A standing wave in a rectangular basin........................................................... 92
Figure 6.3 Simulations in the inviscid flow with different wave steepness...................... 94
Figure 6.4 Viscous damping effects in the shallow water model ..................................... 95
Figure 6.5 The initial 2D standing wave and its decomposed x and y direction 1D wave96
Figure 6.6 Computational domain (not to scale) .............................................................. 97
Figure 6.7 Normalized surface elevation contours in the computational domain at time
4T. ........................................................................................................................... 100
Figure 6.8 Normalized surface elevation contours in the computational domain at time
15T. ......................................................................................................................... 101
Figure 6.9 Normalized Uvelocity contours in the computational domain at time 4T... 102
Figure 6.10 Normalized Uvelocity contours in the computational domain at time 15T.
................................................................................................................................. 103
Figure 6.11 Normalized Vvelocity contours in the computational domain at time 4T. 104
Figure 6.12 Normalized Vvelocity contours in the computational domain at time 15T.
................................................................................................................................. 105
Figure 6.13 Surface displacement at points (0.5, 0) and (0, 0.5) simulated using three
TCS solutions.......................................................................................................... 106
Figure 6.14 Snapshot of the standing wave at 0.5 T. The surface displacements at the two
monitored points (0.5, 0) and (0, 0.5) are circled. .................................................. 107
Figure 6.15 Snapshot of the standing wave at 15 T. The surface displacements at the two
monitored points (0.5, 0) and (0, 0.5) are circled. .................................................. 108
Figure A.1 Initial water level, H (x, 0), in the rectangular open channel ....................... 115
Figure A.2 Inlet boundary condition H (0, t).................................................................. 116
Figure A.3 Shape function of viscosity........................................................................... 117
Figure A.4 An open boundary rectangular channel with a sponge layer........................ 117
Figure A.5 Schematic illustrations of grids .................................................................... 118
Figure A.6 Water level at point A with/without the sponge layer.................................. 118
xi
Figure A.7 Water level evolutions inside the test section and sponge layer................... 119
Figure A.8 Comparison of wave shape at time 7T and 14T for TCS solution1 and TCS
solution 2 respectively ............................................................................................ 120
Figure A.9 Comparison of wave shape simulated from TCS solution 1 and 2 at time 7T
and 14T respectively............................................................................................... 121
Figure A.10 Time evolution of a progressive wave at different wave periods from TCS
solution 1................................................................................................................. 122
Figure A.11 Time evolution of a progressive wave at different wave periods from TCS
solution 2................................................................................................................. 123
1
Chapter 1 Introduction
Models of river flow using mathematic tools have been developed since last half
century (Cao and Carling 2002). Various numerical river models have been created to
study flood control (Dutta et al. 2007; Liao et al. 2007), water allocation (Fleckenstein et
al. 2006; Luo 2007), sediment and morphological evolution (Giri and Shimizu 2006; Le
et al. 2006) and so on. The river flow is usually modeled at a coarsegrid scale due to the
limitations of the computational power (Dedong et al. 2007; Shaw et al. 2005; Thouvenin
et al. 2007). However, there exist some particular river systems such as a river with many
large woody debris (LWD), as shown in Figure 1.1. The flow information around LWD
is at subgridscale, which can?t satisfactorily be modeled at the conventional coarsegrid.
In addition, the river flow is often simulated using twodimensional (2D) depth
averaged shallow water equations (SWE) as the governing equations (Chen and Peng
2006; Nguyen et al. 2006; Weerakoon et al. 2003). However, the unsteadiness and the
dominant advective mechanism of a river flow create difficulties in SWE based
numerical simulation because of the combination of the ?stiffness? (i.e. the unsteady
term) and ?nonlinearity? (i.e. the nonlinear advection term) in the momentum equations.
Therefore, in this research, we develop a new numerical algorithm to address the
difficulties in simulating the unsteady nonlinear advection. This numerical scheme can
also serve as a platform to perform the new conceptual model which can integrate the
subgridscale physics into a coarsegrid scale model. The detailed development of the
conceptual model can be found in Fu and Hodges (2005). First, let us start with the
background information of the LWD.
1.1 BACKGROUND
Large woody debris (LWD) refers to woody material such as fallen tree trunks or
root balls that become lodged in stream channels. Figure 1.1 gives an example of a river
laden with large woody debris (LWD). The conventional size of LWD is diameters
larger than 0.1 m and lengths greater than 1.0 m (Keller and Swanson 1979; Andrus et al.
1988). LWD accumulations have been historically called ?snags?, a pejorative term
reflecting a perception of LWD as a nuisance to river navigation and efficient water use.
River ?improvement? schemes typically involve removing LWD (Shields & Nunnaly
1984; Gippel 1995), to improve water conveyance, rejuvenate channels, lessen the risk of
damage to bridges, improve recreational amenity, and remove barriers to fish migration
(Harmon et al. 1986). However, from an ecological perspective, LWD provides a stable
substratum for microorganisms, algae and invertebrates (Minshall 1984; Brown & May
2000; Statzner & Higler 1986). This, in turn, creates grater local productivity for higher
trophic levels such as invertebrates and fish (Price & Lovett 2002). Furthermore, LWD
enhances hydraulic diversity (i.e. a wide range of flow conditions), which in turn
enhances diversity in fish habitat, e.g. providing lowvelocity refuges sought by many
fish species (Benke et al. 1985; Matthews & Hill 1980), and also providing cover from
predators (Angermeier & Karr 1984, Everett & Ruiz, 1993), foraging habitat, and
2
spawning substratum  all of which vary dynamically with flow rate (Bao & Mathews
1991; Mathews & Tallent 1997; Marzolf 1978). In short, LWD creates greater habitat
complexity (O?Conner 1991), which should produce greater biodiversity (Gerhard &
Reich 2000).
Figure 1.1 Sulphur River (Texas) large woody debris at lowflow conditions
(Photo courtesy of Texas Water Development Board)
Over the past two decades, a wide variety of studies demonstrate how woody
debris may be a dominant player in aquatic habitat. Angermeier & Karr (1984)
selectively modified a stream by inducing and removing woody debris, subsequently
showing a greater abundance of fish and benthic invertebrates correlated with the
introduced debris. Benke et al. (1985) examined a lowgradient stream, showing that only
4% of habitat surface was in LWD and yet this supported 60% of the invertebrate
biomass and 16% of overall production. Similar results are noted in Benke et al. 1984 and
Jacobi & Benke1991. Effects of flooding on fish movements and distribution (Hill &
Grossman 1987; Harvey et al. 1999) have confirmed the importance of habitat
complexity (including LWD) for population persistence. More subtle and indirect effects
on habitat have been attributed to LWD?s modification of the flow field and water level
(Triska & Cromack 1980; Sedell et al. 1982; Wallace & Benke 1984). Near LWD,
reduced velocities allow retention of organic matter, building up ecologicallyimportant
debris dams, while water lever effects at high flow rates can influence the seed dispersal
from riparian vegetation (Merritt & Wohl 2002). Studies drawing similar conclusions
have been conducted over a variety of locales and river types, ranging from salmon
spawning streams in Alaska (Dolloff 1986), lowgradient rivers in West Virginia (Lobb
& Orth 1991) to an ephemeral river in Africa (Jacobson et al. 1999). The common theme
is that LWD plays a varied and critical role in aquatic habitat for rivers and streams ? a
role determined by the flow field near LWD. Moreover, water management agencies are
3
interested in including the effects of LWD in aquatic habitat assessments used for water
resource allocations.
The flow structure around LWD is very complicated and varied, which has been
revealed by several previous field studies (e.g. Mutz 2000; Beebe 2001; Daniels &
Rhoads 2003). Moreover, there is a large difference in scales between hydraulic
processes and the critical ecological processes. Because of this scale difference, the
smallscale flow heterogeneity (important to biota) is poorly represented in the channel
scale numerical models currently used for aquatic habitat analysis (e.g. Waddle 2001). To
address this problem, a new conceptual model is proposed to integrate the subgridscale
physics into a coarsegrid scale model (Fu and Hodges, 2005).
1.2 CHALLENGE AND MOTIVATION OF THE NUMERICAL MODEL
In addition to the scale difference, another challenge for solving the river flow
with LWD is to build a numerical model to simulate the unsteady river flow. River flow
is usually characterized as shallowwater flow. As shown in Figure 1.1, the flow around
LWD is indeed shallow as most of the LWD are exposed out of the water surface. This
kind of river system is usually modeled by 2D depth averaged shallow water equations
(SWE) (Arega et al. 2008; Burguete et al. 2006; Gejadze and Monnier 2007). 2D depth
averaged SWE can be obtained by integrating the NavierStokes equation over the water
depth with the underlying hydrostatic assumption. One of the difficulties in solving SWE
numerically is the combination of the ?stiffness? (i.e. the unsteady term) and
?nonlinearity? (i.e. the nonlinear advection term) in the momentum equations. Numerous
finitedifference numerical models for solving SWE have been developed over the past
several decades. Among these developed numerical models, explicit and semiimplicit
schemes are the most widely used. For example, in the atmospheric or weather research
community, explicit methods such as Leapfrog and LaxWendrofftype are the most
popular ones (MendezNunez and Carroll 1993). Casulli (Casulli 1990) proposed a semi
implicit scheme, on which a number of variations are based. A common property for
both the explicit and semiimplicit methods is that they explicitly discretize the advection
term. Hence, the nonlinearity does not arise in an explicit or a semiimplicit sachem.
However, their stability is restricted by the CourantFriedrichsLewy (CFL) condition
because of the explicitness. A fully implicit scheme is desirable because we can use a
relative large time step in the simulation. A fully implicit solution for an unsteady
nonlinear advection problem involves solving a system of nonlinear equation at each
step. The conventional techniques for solving a system of nonlinear equations require
either 1) iterations of rootfinding procedure or/and 2) calculations of the Jacobian. These
two processes make a fully implicit solution more complex and computationally
expensive. These motivate this research to develop a new implicit scheme to solve
unsteady nonlinear advection problems.
4
1.3 OBJECTIVE
The primary objective of this research is:
Develop a new implicit solution for unsteady nonlinear advection problems. The
new numerical algorithm should have the advantage in both stability and efficiency while
keeping the 2
nd
order temporal accuracy as in most existing SWE models. In addition,
this new method would provide a novel approach in decoupling a coupled equation
system.
1.4 APPROACH
To achieve the objective of this research, several steps need to be taken to guide
this study. The major steps are as follows:
Literature review
? Review current finitedifference models for the 2D depth averaged SWE.
? Examine how the nonlinear advection term is treated in different temporal
discretizations.
? Review the merits and limitations of the existing models.
This part of the dissertation work is presented in Chapter 2.
Analytical development of the new numerical algorithm
Illustrate the theoretical principle and concept of the new algorithm. This part of
the dissertation work is presented in Chapter 3.
Verify the new method using 1D test case
? Test the new algorithm using a 1D PDE with time dependent nonlinear
advection term.
? Test the new algorithm using 1D ordinary differential equation (ODE).
? Explore the characteristics of the method in accuracy, stability and
efficiency.
This part of the dissertation work can be found in Chapter 4.
Verify the new method using a multivariable and multidimensional problem
? Theoretically develop the new algorithm in multivariable and multi
dimensional problems.
? Demonstrate and analyze the decoupling process of the new method.
? Conduct initial test cases of the numerical algorithm in 2D depth averaged
SWE.
This work can be found in Chapter 5 and Chapter 6.
5
1.5 BROAD APPLICATION
The initial motivation and application for this research is simulating a river flow
with LWD. The newly developed numerical algorithm for solving unsteady nonlinear
advection problems is not limited only in simulating river flows. In most of the
environmental flows such as atmospheric flow, ocean circulation, storm surge and etc.,
advection is the dominant mechanism and we are often interested in simulating these
flows with a range of different time scales. Therefore, the difficulties arising from the
combination of the unsteadiness and nonlinear advection exist in modeling many types of
environmental flows. This new numerical method provides a new approach to address
this difficulty in these and many other research areas.
6
Chapter 2 Review of Finite Difference Schemes for Solving the
Unsteady Nonlinear Advection in the Shallow Water Equations
Environmental flows (e.g. river and atmospheric flows) governed by the SWE are
usually dominated by the advective mechanism over multiple timescales (Vreugdenhil
1994). The combination of time dependency and nonlinear advection creates difficulties
in the numerical solution of the shallow water equations (SWE) (Bourchtein and
Bourchtein 2006). To address this difficulty, numerous numerical schemes have been
developed over the past several decades. Iskandarani et al (2005) reviewed the finite
element/finite volume methods in solving advection equations. More reviews about finite
element methods can be found in the article by Thomee (2001). The research herein
focuses on developing a finite difference algorithm for solving unsteady nonlinear
advection problems. To review the existing finite difference techniques, we first classify
timemarching algorithms into three categories: explicit, semiimplicit and implicit
methods. In the explicit and semiimplicit methods, the nonlinear advection terms are
treated explicitly. Hence, the issue of solving a nonlinear matrix inversion does not arise
in these two temporal discretizations. However, to solve a fullyimplicit system, different
computational linearization methods are required treat the nonlinearity, so that the
inversion of a nonlinear matrix can be obtained. Therefore, we briefly review the explicit
and semiimplicit models for solving the SWE and examine different linearization
approaches in the implicit system. Comparisons among different existing methods are
summarized in the last section of this chapter.
2.1 EXPLICIT METHOD
Fully explicit methods are popular in atmospheric or weather research
community. The most widely used explicit methods in atmospheric research are Leapfrog
and Laxwendrofftype schemes (MendezNunez and Carroll 1993). The Leapfrog
method is a onestep method, requiring two time levels of known values to compute an
unknown level. Modifications of the Leapfrog methods are used in simulating ocean
circulation (Cho and Yoon 1998; Fujima and Shigemura 2000). The LaxWendrofftype
method is a twostep PredictorCorrector method. One of the popular variations of the
Laxwendroff scheme is developed by MacCormack (1969). The MacCormack scheme is
extended to many research areas such as aerospace engineering (Hirose et al. 1991;
Kurbatskii and Mankbadi 2004) and hydrological science in simulating free surface flow
(Fennema and Chaudhry 1990; Garcia and Kahawita 1986; Kazezy?lmazAlhana et al.
2005; Keshari and Koo 2007; Li and Jackson 2007), groundwater flow (Keshari and Koo
2007) and dambreak shock wave (Li and Jackson 2007). Another commonly used
explicit PredictorCorrector method is the RungeKutta method (Delis and Katsaounis
2005; Zhou et al. 2007). The Leapfrog and the twostep PredictorCorrector methods can
maintain 2
nd
order temporal accuracy. Although simpler explicit schemes such as
7
forwardtime scheme is also used in solving SWE (Murillo et al. 2008), it can only
maintain 1
st
order temporal accuracy.
Explicit schemes are straightforward and easy to implement. Nonlinear advective
terms in partial differential equations (PDE) are readilycomputed from known time
levels using fullyexplicit timemarching methods. However, the implementations of
explicit methods are restricted by stability requirements. RungeKutta and MacCormack
methods typically require an excessively small time step to satisfy the Courant
FriedrichsLewy (CFL) condition(MendezNunez and Carroll 1993); the Leapfrog
approach must be applied with appropriate numerical strategies such as alternate time
levels of the solution (Agoshkov et al. 1994; Peyret and Taylor 1983; Zhou 2002) or
sophisticated modesplitting (e.g.Blumberg and Mellor 1987).
2.2 SEMIIMPLICIT METHOD
Semiimplicit methods have been developed and widely used in temporal
discretized SWE. In most semiimplicit methods, the nonlinear advection terms are
discretized explicitly (Bonaventura and Rosatti 2002; Bourchtein and Bourchtein 2007;
Casulli 1990; Casulli and Cattani 1994; Casulli and Cheng 1992; Kar 2006; Spitaleri and
Corinaldesi 1997) so that the issue of solving a nonlinear matrix inversion does not arise.
For instance, Environmental Fluid Dynamics Computer code (EFDC) uses a threelevel
semiimplicit method, in which advection terms are explicitly discretized using upwind
scheme. As a result, the stability is constrained by the explicit advection terms (Hamrick
1992). Some other models (Kar 2006) use the explicit RungeKutta method to calculate
the advection term. It is obvious that the stability of this treatment is constrained by the
CFL number.
To overcome the stability constraint, different numerical treatments are
introduced to calculate the linearized advection term. One of the most widely used
techniques is the semiLagrangian method (Cullen 2001; Rosatti et al. 2005). Casulli
(Casulli 1990; Casulli and Cattani 1994; Casulli and Cheng 1992) used this method in
UnTRIM and the earlier version, TRIM. Some latest models including ELCIRC (Zhang
et al. 2004) and ELCOM(Hodges 2000) followed this idea. Although Casulli?s approach
has been successful and stable, its accuracy is relatively poor (Hodges, 2004).
The main advantage of the semiLagrangian method is that the stability is not
constrained by the CFL condition (Barros and Garcia 2007). However, a set of trajectory
equations need to be solved at each time step and the maximum allowable time step is
restricted by the convergence criteria of the iteration of the trajectory equations
(Bourchtein and Bourchtein 2007).
2.3 IMPLICIT METHOD
An obvious advantage of a fullyimplicit method is that the time step is not
restricted by the CFL number. Various fullyimplicit finite difference schemes have been
developed. Steppeler (2006) solved a fully implicit, SWE based meteorological model
using the Fourier transformation method at each grid point and each time step. Some
8
researchers derived a fullyimplicit scheme based on Alternating Direction Implicit (ADI)
method (Szymkiewicz 1992; Wilders et al. 1988). Yuan and Wu solved implicit Navier
Stokes equations using a staggered finite difference CrankNicolson scheme (Yuan and
Wu 2004). Burguete and GarciaNavarro (Burguete and GarciaNavarro 2004)
implemented a first order upwind implicit scheme to simulate river hydraulics problem.
Although the fullyimplicit methods are more accurate and robust, they are less
common, possibly due to computational complexity and expense (Turek 1996).
Nonlinearity in a fullyimplicit method results in a system of nonlinear algebraic
equations to be solved at each time step. Existing strategies for solving implicit nonlinear
equations include iterative and noniterative methods (Moin 2001). In the following, we
will review the conventional techniques for implicit nonlinear solutions.
The Newton method and Picard method (Ferziger and Peric 1999; Lehmann and
Ackerer 1998; Paniconi et al. 1991) are the most widely used twolevel iterative
techniques for implicit timemarching of nonlinear equations. For steadystate nonlinear
problems, Newton and Picard iterative algorithms have proven quite successful (Paniconi
and Putti 1994). These methods may be thought of as successive linear solutions of
=Ax b
%
, where A
%
is an approximation of A that is iteratively refined to solve the
nonlinear problem. However, for timeevolving CFD problems, each timemarching step
using the Newton or Picard method requires an outer iteration applied over an inner
solution of a linear equation set. When the inner problem also requires an iterative
solution (as is common in CFD), the time march for this doublyiterative approach is
computationally expensive. The principal differences between the Picard and Newton
methods are that the former is easier to implement and requires fewer computations per
outer iteration but has only 1storder convergence, whereas the latter is more difficult to
implement but provides 2ndorder convergence. Thus, which method is appropriate
depends upon the difficulty in implementing the Newton method compared to the slower
convergence of the Picard method. A further implicit technique, local linearization, has
not been widely used in CFD but does provide timemarching with 2ndorder accuracy
(Lomax et al. 1999) without an outer iteration. By providing a single linear
approximation of the nonlinear problem, local linearization sidesteps the convergence
issue of the outer iteration for Newton/Picard techniques. However, local linearization
requires discretization of the Jacobian, which is often difficult to derive and implement
for typical CFD applications.
To better illustrate these existing implicit linearization methods, we use a simple
scalar ODE with a quadratic nonlinearity as an example. The nonlinear ODE is written as
d
f( ,t)
dt
?
= ? (2.1)
Using a finitedifference CrankNicolson discretization (the simplest 2
nd
order implicit
method), the above can be approximated as
()(){}()
n1 n n1 n 3
t
ff t
2
++
?
?=?+ ? +?+?O (2.2)
9
where superscripts indicate the discrete time step and
( )
n1
f
+
? implies a quadratic
nonlinear relationship in
n1+
? (e.g.
n1+
?
n1+
? or
n1
x
+
?
n1+
? ). Equation (2.2) is our example
of an implicit nonlinear equation. In the following, we will demonstrate the principles of
existing techniques to linearize Equation (2.2) computationally.
2.3.1 Newton methods
The Newton method is one kind of rootfinding algorithm. It starts with an initial
guess and iteratively estimates the root of the function. The Newton method can be
derived from a Taylor series expansion or a geometric proof (Amat et al. 2003). One of
the advantages of the Newton method is that it can be applied to equation systems with
complex nonlinearities and keeps the quadratic convergence rate. In petroleum
engineering and other research areas concerning flow through porous media, the Newton
method is a standard approach because of the complex nonlinearities (Cao and Sun 2005;
Dettmer and Peric 2007; Kwok and Tchelepi 2007). Using the Newton method, the
linearized equation system of Equation (2.2) can be written as
()()
( )
k
n1
k1 k
n1 n1
k
n1
g
g
+
+
++
+
? ?
?
? ?
?=??
???
??
??
??
(2.3)
where the outer superscript indicates the iteration number and
() ()( ){}
n1 n1 n n n1
t
gff
2
++ +
?
?=???? ?+? (2.4)
In general, we can use the value at the previous time step as the value at the first iteration,
which means
( )
1
n1 n+
? =? (2.5)
Equation (2.3) is solved repeatedly until the difference between two successive iterations
satisfies the predefined convergence criteria. In Equation (2.3), a first order derivative
()
n1
g/
+
??? must be calculated and updated at each iteration. If ?is a vector and
discretized in space, the resulting
( )
n1
g/
+
??? is a matrix and called the Jacobian of the
linearized equation system. However, the Jacobian matrix may be computationally
expensive and difficult to calculate analytically (Ferziger and Peric 1999; Niet et al.
2007).
A number of approaches have been proposed to modify Newton method. For
example, it can be modified with a Chebyshev approximation to accelerate the
convergence (Bagatur 2007). To simplify the Newton method, much research has been
carried out to simplify the calculation of the Jacobian matrix. Niet et al. (2007) evaluated
the Jacobian matrix using a partialanalytical and partialnumerical technique in an
oceanclimate model. Li (1993) attempted to reduce the effort to compute the derivatives.
10
A Jacobianfree NewtonKrylov (JFNK) method has been developed recently to solve
nonlinear equation systems (Brown and Saad 1990; Chan and Jackson 1984; Knoll and
Keyes 2004; Mousseau et al. 2002; Reisner et al. 2005; Wubs et al. 2006). The key to
JFNK solver is approximating the Jacobianvector product iteratively instead of
evaluating each element of the Jacobian matrix. Although all these modifications reduce
the computation effort of Jacobian, inconvenient iteration still remains in the Newton
method.
2.3.2 Picard Iteration
The Picard iteration is more straightforward than the Newton method. In the
Picard iteration, the previous outer iteration value,
( )
k
n1+
? , is substituted for
n1+
? on the
RHS of Equation (2.2), resulting in
() () ( )
{ }
k1 k
n1 n n n1
t
ff
2
+
++
?
? ?
?=?+?+?
? ?
? ?
(2.6)
As in the Newton iteration, Equation (2.6) is solved repeatedly until the difference
between two successive iterations satisfies predefined convergence criteria. The same
initial condition, Equation (2.5) is used to start the solution. Due to the simplicity, Picard
iteration has been widely used in fully implicit nonlinear solver in Computational Fluid
Dynamics (CFD)(Clement et al. 1994; Kwag 2000; Webster 2007). Different
modifications of the Picard method have been reported in the literature. For example,
Celia et al.(1987) proposed a modified Picard iteration that is a combination of Picard
iteration and Newton iteration. However, this and other modified Picard methods remain
linearly convergent (Lehmann and Ackerer 1998). A detailed comparison of the Picard
and Newton methods is provided by Paniconi and Putti (1994). In general, Newton
methods require more computations per iteration, but converge more rapidly (quadratic)
than Picard methods (linear). However, for advection equations, the overall
computational cost of Newton methods is typically greater than Picard methods because
of the computational cost of the Jacobian (Ferziger and Peric 1999).
2.3.2 Local linearization
In contrast to the abovementioned Newton and Picard methods, local
linearization is a noniterative method. Let?s first review the derivation of the local
linearization techniques, which are based on a Taylorseries expansion for
n1
f
+
about
n
t in
Equation (2.2). The following derivation is based on Lomax et al.(1999),
()
n
n
n1 n 2
ff
ff t t
t
+
????
??
=+?? +? +?
????
?? ?
??
??
O (2.7)
11
where
n1 n+
??=? ?? . For advection equations, f is generally not an explicit function
of t, so ()
n
tf/t??? in Equation (2.7) is zero. Substituting Equation (2.7) into Equation
(2.2) provides
()
n
n1 n n n 3
tf
fft
2
O
+
??
??????
?=?+ +?? + +?
??
??
??
??
??
(2.8)
which is equivalent to
n
n
tf
1tf
2
??
????
? ??=???
??
??
??
??
(2.9)
Equation (2.9) is a linear secondorder discrete form of the original nonlinear ODE and
can be solved by any number of standard techniques once
()
n
f/? ?? is explicitly
computed. However, as in the Newton method, the evaluation of the Jacobian ()
n
f/???
complicates the overall computation.
2.4 SUMMARY
SWE is widely used in simulating environmental flows. Mathematically, the
stiffness and nonlinear advection increases the difficulty in the numerical solutions of
SWE. Explicit methods can easily solve the nonlinearity but are restricted by the CFL
stability criteria. Semiimplicit methods treat the advection term explicitly but need
additional numerical treatment such as the semiLagrangian method to calculate
advection terms for a higher CFL number. A fully implicit method is desirable but
requires further computational linearization to solve a nonlinear system of equations.
Existing implicit linearization techniques (including Newton, Picard and local
linearization methods) require either 1) additional explicit derivative evaluations to
provide an approximate linear problem, or 2) an outer iteration that converges an inner
approximate linear problem. These two characteristics make implicit systems complex
and expensive to solve. To address the nonlinearity in a fully implicit system, in the next
chapter we propose a new computational linearization method that can linearize the
nonlinear advection term without either 1) the outer iteration or 2) computation of the
Jacobian.
12
Chapter 3 Theoretical Development of the TimeCentered Split (TCS)
Method
In Chapter 2 we reviewed the current implicit linearization methods that require
either an outer iteration or a computation of the Jacobian. In this chapter, we develop a
new method that is similar to local linearization in that it allows noniterative
discretization of a temporally 2
nd
order approximation of the nonlinear equation set.
Instead of requiring the Jacobian, the new method splits the timemarching nonlinear
problem into two sets of linear problems that are solved in succession. The new method
has both the computational efficiency of a noniterative local linearization method and
the implementation simplicity of a Picard iterative method. We call this the Time
CenteredSplit (TCS) method.
The theory of the TCS method is first derived using a generic single variable
nonlinear equation in this chapter. The TCS method can generate different discrete forms
by introducing the timecentered split to different terms. This advantage is demonstrated
by applying the TCS method to a 1D advectiondiffusion equation (Burgers? equation).
Four different TCS formats such as TCSF (split the flux term), TCSG (split the gradient
term), TCSFD (split the flux and diffusion terms) and TCSGD (split the gradient and
the diffusion terms) are derived and presented for the 1D Burgers? equation. One of the
key advantages of the TCS method is that it provides noniterative coupling between
multiple variables. A 1D coupled advection diffusion equation and scalar transport
equation are used as an example to illustrate this advantage. A summary of the principles
of the TCS method and its advantages are provided in the last section.
3.1 THE THEORY OF TIMECENTERED SPLITTING
3.1.1 Computational Splitting of the Nonlinear Term
The TCS method is based on a nested application of the midpoint rule (i.e. a
centeredtime approximation). Midpoint rule discretizations are often used for time
marching to obtain secondorder temporal accuracy (Ferziger and Peric 1999). A
common approach for nonlinear timemarching in meteorology and oceanography is
application of the midpoint rule in an explicit formulation known as the 3level Leapfrog
method (Dubois et al. 2005; Fujima and Shigemura 2000). The explicit Leapfrog method
can be written as
( )
n1 n1 n n n 3
f2tt
+?
?=?+????+?,O() (3.1)
where superscripts represent time levels and t? is the model time step. Equation (3.1)
can be seen as a generic form of a nonlinear equation. The function f is a linear operator
of ?? and ?. Inside the function f, the product term ?? represents a generic quadratic
nonlinear term. This quadratic nonlinear term in a flow and transport problem is of the
13
interest in this research. To develop the TCS method, we note that in the same vein as
Equation (3.1), the midpoint rule can be written across only two time levels as
( )
n1 n n12 n12 n12 3
ftt
++++
?=?+? ? ? ?+?
// /
,O() (3.2)
Equation (3.2) 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 so that the advance is not leapfrogging
over alternating data. By introducing a timecentered linear approximation of
( ) ( )
n1/2 n n1 2
/2 tO
++
?=?+? +? (3.3)
in one part of the nonlinear product
n 1/2 n 1/2++
??in Equation (3.2), the discrete equation
becomes computationally linear in time (i.e. no products with the same time level):
()
nn1
n1 n 2 n12 n12 3
fOt tt
2
+
+++
??
???+?
?=?+ +? ? ? ?+?
??
??
????
//
,O() (3.4)
Reorganizing Equation (3.4) provides
()()
n1 n n n12 n12 n1 n12 n12 3
tt
ff t
22
+++++
? ?
?=?+?? ? +?? ? +?
// //
,,O() (3.5)
Equation (3.5) is computationally split into two steps and an intermediate variable
*
? is
defined as,
()
* n n n 1/2 n 1/2
t
f,
2
++
?
?=?+ ?? ? (3.6)
Subtracting Equation (3.6) from Equation (3.5) provides
()
n1 n1 n12 n12 3
t
ft
2
++++
?
?=?+?? ? +?
*//
,O() (3.7)
After introducing the above timecentered split, Equations (3.6) and (3.7) are both
computationally linear. Furthermore, Equation (3.6) is similar to the implicit Euler
approximation of
n1/2+
? written as:
()
2
n12 n n12 n12 n12
tt
f
22
////
,O
++++
? ?
??
?=?+??? +
??
??
(3.8)
Using the Taylor expansion,
n1/2+
? can be expanded as:
2
n12 n
t
tt
22
/
O
+
??
??
?=?+?+
??
??
(3.9)
Substituting Equation (3.9) into the nonlinear term of Equation (3.8), provides
2
n12 n n n12 n12
ttt
f
22
///
O,O
+++
?????
??
?? ??
?=?+?+ ?? +
??
?? ??
??
?? ??
????
(3.10)
Grouping the higher order terms,
14
()
2
n12 n n n12 n12
tt
f
22
///
,O
+++
? ?
??
?=?+??? +
??
??
(3.11)
Substituting Equation (3.6) into Equation (3.11), results in
2
n12
t
2
+
?
??
?=?+
??
??
/*
O (3.12)
Substitution of Equation (3.12) into Equations (3.6) and (3.7) provides a twostep
method
()
nn 3
t
ft
2
?
?=?+ ??? + ?
***
,O() (3.13)
()
n1 n1 3
t
ft
2
++
?
?=?+??? +?
***
,O() (3.14)
Equations (3.13) and (3.14) are a computationally linearized implicit equation system.
The summation of Equations (3.13) and (3.14) is a 2
nd
order computational equivalent of
Equation (3.2).
Remark:
Although the derivation above provides a twostep direct linearization method, an
outer iteration might also be combined for a big t? value. This process can be illustrated
as the following:
1. obtain
n1+
? from Equations (3.13) and (3.14).
2. calculate
n1/2+
? using
( )
n1 n
/2
+
?+? , where
n1+
? is calculated from step 1.
3. calculate a new
n1+
? using Equations (3.13) and (3.14) with the updated
n1/2+
? in step 2.
Repeat the procedure until the difference between the old and new
n1+
? is
acceptable.
3.1.2 Computationally Splitting the Linear Term
In the previous section, we introduced timecentered split only into the nonlinear
term of ??. In addition to splitting the nonlinear product term, one can also use the same
splitting idea in the linear term?, following
() ()
nn1 nn1
n1 n 2 n12 2 3
ft tt(t)
22
/
O,O O
++
++
??
?????+? ?+?
?=?+ +? ? +? ?+?
??
????
??????
(3.15)
Equation (3.15) is guaranteed 2
nd
order accurate in time because of Equation (3.3).
Reorganizing Equation (3.15) and grouping the higher order terms,
15
{}
n1 n n n12 n1 n12 n n1 3
t
f()
2
//
,,, O
++++
?
?=?+?? ?? ?? +? (3.16)
In Equation (3.16), not only the nonlinear term but also the linear term is split into two
parts. Therefore, a different intermediate variable
*
? is defined as
()
*n nn1/2n
t
f,
2
+
?
?=?+ ?? ? (3.17)
The second step of the splitting system then becomes,
()
n1 * n1 n1/2 n1
t
f,
2
++++
?
?=?+?? ? (3.18)
Equations (3.17) and (3.18) are different from the split system in the previous section. As
a result, the correspondence between
*
? and
n12/+
? is proved differently. Instead of the
implicit Euler equation, the explicit Euler approximation for
n12/+
? is introduced,
()
2
n12 n n n n
tt
f
22
/
,O
+
??
??
?=?+??? +
??
??
(3.19)
Substituting the Taylor expansion of
n12/+
? into Equation (3.17) and grouping the higher
order terms,
()
2
*n nnn
tt
f, O
22
??
??
?=?+ ??? +
??
??
(3.20)
Substituting Equation (3.20) into Equation (3.19), we obtain
2
n12
t
2
+
?
??
?=?+
??
??
/*
O (3.21)
Thus, we can use
*
? to replace
n12/+
? in Equations (3.17) and (3.18). A different set of
computationally linearized twostep equations can be written as:
()
*n n*n
t
f,
2
?
?=?+ ??? (3.22)
()
n1 * n1 * n1
t
f,
2
+++
?
?=?+??? (3.23)
The key to the TCS method is the secondorder timesplitting of quadratic
nonlinear terms to two different time levels, i.e.
n*
? ? and
n1 *+
? ? . The result is discretely
linear in any timelevel of information. A further application of the same splitting in the
linear terms will give another set of discrete linearized equations. The above derivation
for the single variable ? can be readily extended to a vector of variables,
12 N
[ , ,... ]?? ? , or
variables and linear operators, e.g.
112 2N 2
[ ,L( ), ,L( )... ,L( )]? ?? ? ? ? . However, as
additional variables (or operators) are introduced, the discretization method has multiple
implementations. For example, even with a 1D advection diffusion equation, we can
16
obtain at least four different TCS discrete formats. In the next section, we will present
these four different TCS formats by applying the TCS method to a 1D Burgers? equation.
3.2 DERIVATION OF THE TCS METHOD IN A 1D ADVECTIONDIFFUSION EQUATION
The Burgers? equation is the simplest advectiondiffusion test case and can be
viewed as a prototype of the NavierStokes equation or the SWE. The nonconservative
Burgers? equation can be written as:
2
2
uu u
u
txx
? ??
+=?
? ??
(3.24)
To develop the TCS method, we note that in the same vein as Equation (3.7), the
midpoint rule can be written across only two time levels as
()( )
n1 n n12 n12 2 n12 3
xx
uu u u u t t
// /
O
++++
=+? ? +?? ?+? (3.25)
where
x
? is the shorthand notation of the generic discretized spatial derivative of ?u?.
3.2.1 The TCSF Discrete Format
The nonlinear advection term in Equation (3.24) is constructed by the flux part ?u?
and the gradient part, ? ux/???. The timecentered splitting maybe applied to either part
of the advection term. We first introduce a timecentered linear approximation of
( ) ( )
n1/2 n n1 2
uuu/2tO
++
=+ +? to the flux term in the nonlinear product in Equation
(3.25). The discrete equation becomes computationally linear in time (i.e. no products
with the same time level):
()
n1 n
n1 n n12 2 n12 3
xx
uu
uu u u t t
2
//
O
+
+++
??
??+
=+? ? +?? ?+?
??
??
????
(3.26)
Multiplying the products in Equation (3.26) provides
()()()()
n1 n n1 n1/2 n n1/2 2 n1/2 3
xxx
tt
uu uu uu u tOt
22
+++ + +
? ?
=+?? +?? +?? ?+? (3.27)
A computationalsplitting technique is used to obtain a numericallysolvable set of
equations from the discrete form of Equation (3.27). Defining a generic intermediate
variable
*
u as
()()
n n n12 2 n12
xx
tt
uu uu u
22
*//++
? ?
=+?? +?? (3.28)
and subtracting Equation (3.28) from Equation (3.27) provides
()()()
n1 n1 n12 2 n12 3
xx
tt
uu uu u t
22
*/ /
O
+++ +
? ?
=+?? +?? +? (3.29)
17
Using an implicit Euler approximation of
n12
u
/+
and the same mathematical derivation as
in Equation (3.12), we can prove
2
n1/2 *
t
uu
4
O
+
???
=+
??
??
(3.30)
Substitution of Equation (3.30) into Equations (3.28) and (3.29) provides a twostep
method, which is called the TCSF method:
{}
nn2
xx
t
uu uu u
2
***
?
=+ ??+?? (3.31)
{}()
n1 n1 2 3
xx
t
uu uu u t
2
***
O
++
?
=+ ??+?? +? (3.32)
The summation of Equations (3.31) and (3.32) is equivalent to the original Equation
(3.25). The first step is an implicit solution for the variable u
*
involving only
n
u and with
an implicit discretization of the diffusion term.
n1
u
+
in the second step can be explicitly
calculated.
3.2.2 The TCSG Discrete Format
Splitting the gradient term instead of the flux term in Equation (3.25) results in a
second basic form of TCS. Thus, instead of Equation (3.26), we obtain
n1 n
n1 n n1/2 2 n1/2 3
xx
uu
uu u u t(t)
2
O
+
++ +
??
??+
=+? ? +?? ?+?
??
??
????
(3.33)
For this second form, we define a slightly different intermediate variable as
()()
* n n1/2 n 2 n1/2
xx
tt
uu u u u
22
++
? ?
=+? ? +?? (3.34)
Similarly in Equation (3.30), we have
2
n1/2 *
t
uu
4
O
+
???
=+
??
??
(3.35)
Subtracting Equation (3.34) from (3.33) and substituting Equation (3.35) in the same
manner as the transition from Equation (3.27) through (3.32), we obtain a second discrete
split form that we will call TCSG
{}
*n *n 2*
xx
t
uu uu u
2
?
=+ ??+?? (3.36)
{}
n1 * * n1 2 *
xx
t
uu uu u
2
++
?
=+ ?? +?? (3.37)
The summation of Equations (3.36) and (3.37) is equivalent to the original Equation
(3.25). The first step is an implicit solution for the variable u
*
involving
n
u and with an
18
implicit discretization of the diffusion term. The second step is an implicit solution for
the variable
n1
u
+
involving only u
*
and with an explicit discretization of the diffusion
term.
3.2.3 The TCSFD Discrete Format
The TCSF and TCSG forms are obtained by introducing the timecentered split
into the advection term. Further approximation of the diffusion term using the same time
splitting techniques will create different discretizations. For instance, in Equation (3.26),
we can substitute the time splitting into both the flux and diffusion terms,
()
n1 n n1 n
n1 n n12 2 3
xx
uu uu
uu u t t
22
/
O
++
++
??
?? ??
=+? ? +?? ?+?
??
?? ??
?? ????
(3.38)
Defining u
*
as,
()()
nnn12 2n
xx
tt
uu uu u
22
*/+
? ?
=+?? +?? (3.39)
the second step follows as,
()()()
n1 n1 n12 2 n1 3
xx
tt
uu uu u t
22
*/
O
+++ +
? ?
=+?? +?? +? (3.40)
Using an explicit Euler approximation of
n12
u
/+
and the same mathematical derivation as
in Equation (3.21), we can prove that
2
n1/2 *
t
uu
4
O
+
???
=+
??
??
(3.41)
Substituting Equation (3.30) into Equations (3.39) and (3.40), we obtained the TCSFD
format as
{}
nn2n
xx
t
uu uu u
2
**
?
=+ ??+?? (3.42)
{}()
n1 n1 2 n1 3
xx
t
uu uu u t
2
**
O
+++
?
=+ ??+?? +? (3.43)
We call Equations (3.42) and (3.43) the TCSFD method because we apply the time
centered splitting to both the flux term and the diffusion term. The summation of
Equations (3.42) and (3.43) is equivalent to the original Equation (3.25). The first step is
an implicit solution for the variable u
*
involving
n
u and with an explicit discretization of
the diffusion term. The second step is an implicit solution for the variable
n1
u
+
involving
only u
*
and with an implicit discretization of the diffusion term.
19
3.2.4 The TCSGD Discrete Format
Similar to the development of the TCSFD method, substitution of timecentered
splitting into the diffusion term in the TCSG provides another set of twostep equations
as
{}
*n *n 2n
xx
t
uu uu u
2
?
=+ ??+?? (3.44)
{}
n1 * * n1 2 n1
xx
t
uu uu u
2
+++
?
=+ ?? +?? (3.45)
Equations (3.44) and (3.45) are called the TCSGD discrete form. The summation of
Equations (3.44) and (3.45) is equivalent to the original Equation (3.25). u
*
in the first
step can be explicitly calculated and the diffusion term is explicitly discretized. The
second step is an implicit solution for the variable
n1
u
+
involving only u
*
and with an
implicit discretization of the diffusion term.
We have derived four discrete forms of the TCS method applied to a 1D Burgers?
equation in the previous sections. More possible discretizations will be generated when
we apply the TCS method to an equation system with multiple variables. We will discuss
this characteristic in Chapter 5.
3.3 TCS FOR COUPLED MOMENTUM AND SCALAR TRANSPORT
A key advantage of the TCS method is that it provides noniterative coupling
between multiple variables. The method is best understood through application to a
simple model problem. First, we will consider the 1D advectiondiffusion equation for a
scalar ?
2
2
u
txx
??????
+=?
? ??
(3.46)
where u is the velocity and ? is a diffusion coefficient. Equation (3.46) can be coupled to
the 1D Burgers? equation for momentum
2
2
uu u
u
txx
? ??
+=?
? ??
(3.47)
Applying TCSF to Equation (3.46) provides
2*
nn
2
t
u
2xx
?
?
? ?? ?? ? ?
?=?? ??
? ?
??
??
(3.48)
2*
n1 n1
2
t
u
2xx
?
+? +
? ?? ?? ? ?
?=?? ??
? ?
??
??
(3.49)
The first step, Equation (3.48), is an implicit equation for
*
? involving only time ?n?
values of u, whereas the second step, Equation (3.49), linearly couples solution of
n1+
? to
20
n1
u
+
. To complete the algorithm, Equation (3.47) can be similarly discretized so that the
coupled TCSF method for both scalar transport and momentum can be written as linear
operators,
() ()
2
nn
2
t
1u
2x x
?
??
????
? ??
+???=?
????
??
??
??
(3.50)
() ()
2
nn
2
t
1u uu
2x x
?
??
????
???
+??=
????
??
??
??
(3.51)
2*
n1
2
2*
n1
2
tt
1
2x 2 x
tu t u
u01 u
2x 2 x
?
+?
?
+ ?
????????? ? ??
??+?
??????
=
??
?? ? ?
++?
??
??????
(3.52)
where the parentheses indicate a spatial derivative operator on the term to the right of the
brackets. The first step of the TSCF method for N variables over Q grid cells results in N
independent linear problems of order Q. The second step of TCSF results in a single
linear problem of order NQ. This series of linear problems is a 2
nd
order temporal
equivalent of the original timemarching, coupled, nonlinear momentumadvection
problem. As the TCS uses an intermediate solution ( ,u )
? ?
? followed by final solution
n1 n1
(,u)
++
? , it resembles predictorcorrector schemes. However, classic predictor
corrector methods are formulated with an explicit predictor of the time n+1 values
followed by an implicit corrector to the time n+1 values, which makes the classic
predictorcorrector methods restricted by the CFL condition; furthermore, predictor
corrector methods do not provide an avenue for a simple nonlinear solution (see
discussions in Lomax et al, 1999 and Tannehill et al, 1997). In contrast, the new TCS
method provides an implicit linearized predictor of the time n+1/2 values that are used in
a coupled implicit solution of the time n+1 values. This new approach provides a simple
method of linearlycoupling equations that are nonlinearlycoupled in the original
problem. The method requires only linear matrix solutions and does not require the outer
iteration of Picard or Newton methods. As compared to the functional Jacobians required
for local linearization and Newton iteration, the TCS coefficient matrices are relatively
easy to derive and have forms very similar to the original model problem.
3.4 SUMMARY
A new computational linearization method, the TCS, is derived and analyzed in
this chapter. Without iteration and calculation of the Jacobian, the TCS method splits the
quadratic nonlinear term into two steps so that each step is computationally linear. That
is, for time marching from known state
n
x to unknown state
n1+
x we use a linear solution
of
n
()=Ax x b
% %
% followed by
n1
()
+
=Axx b
% %
% %
% that is a secondorder equivalent to
21
n1 n1
()
++
=Ax x b. In addition to the linearization, the TCS method can generate different
TCS discrete forms. All different TCS formats are mathematically equivalent to the
original implicit midpoint rule discretization. As a result, they share the same 2ndorder
temporal accuracy. Furthermore, we also demonstrated that the TCS method can couple
multiple variables without iterations. Characteristics of the TCS method such as
accuracy, stability and efficiency will be explored in the next chapter using examples
with analytical solutions.
22
Chapter 4 Implementation of the TCS Method in 1D Problems
The theory of the TCS method is illustrated in Chapter 3. To investigate the
characteristics of the new method, we apply the TCS method to three different test cases:
a 1D conservative Burgers? equation, a 1D nonconservative Burgers? equation and a 1D
nonlinear ordinary differential equation (ODE). For each test case, the TCS algorithm is
compared to the conventional implicit nonlinear solution methods (local linearization,
Picard iteration and Newton iteration) applied to CrankNicolson discretization. The
temporal accuracy of different TCS discretizations is verified by all three test cases. The
practical stability of the TCS method is confirmed using the unsteady flow test case with
an analytical solution in both conservative and nonconservative forms. The method is
shown to require computational effort similar to local linearization, but does not require
discrete computation of a functional Jacobian for solution.
4.1 APPLICATION OF THE TCS METHOD TO THE 1D CONSERVATIVE BURGERS?
EQUATION
4.1.1 Discrete formats of the 1D Conservative Burger?s equation using various
methods
Burgers? equation provides a useful model problem for comparing the TCS to
other nonlinear solution methods. We will begin from the conservative form of Burgers?
equation. The application of the TCS method to the nonconservative Burgers? equation
will be discussed in the next section. The conservative form of the 1D Burgers? equation
is:
()
2
2
2
u1 u
u
t2x x
? ??
+=?
?
(4.1)
TimeCentered Split Method
For an equation with one variable and a simple quadratic nonlinearity, TCSF and
TCSG collapse into a single form. For Equation (4.1), the TCSF and TCSG methods are
identical, and can be presented as linear operators
() ()
n2
n
2
u
t1
1uu
22 x x
?
??
????
? ??
+??=
????
??
??
??
(4.2)
( )
2*
n1
2
u
ttu
1uu
4x 2x
?
+?
???
????
+=+
??
??
??
(4.3)
23
For simplicity in the following discussion, we use the generic symbols A, x and b to
represent a coefficient matrix, the lefthandside (LHS) variable vector and the right
handside (RHS) known vector in a matrix equation of the form of Ax = b for various
solution methods. The number of grid points in space is Q. A central difference spatial
discretization is applied to derivatives in Equation (4.2) and(4.3). It is useful to define a
viscous scale,?,
2
t
x
??
??
?
(4.4)
and the ?i? grid cell CFL number for the
{ }
Ln,,n1? ?+time level for
{ }
i 1,2,...Q= as
L
L i
i
ut
C
x
?
?
?
(4.5)
Similarly, it will be useful to also define a diffusion operator at any time level as
( )
LLLL
ii1ii1
Du2uu
+?
?? ? + (4.6)
and a nonlinear adjective gradient operator
LLL LL
ii1i1i1i1
GCu Cu
+ +??
?? (4.7)
It follows that Eq, (4.2), the first step of TCSF can be written in the form Ax = b over Q
grid points where the A matrix is tridiagonal such that
n
2
nn
13
nn
Q2 Q
n
Q1
C
10
28
CC
1
28 28
CC
1
28 28
C
01
28
?
?
??
?
+? ? +
??
??
?? +? ?+
=
??
?? +? ?+
??
?
?? +?
??
A OOO (4.8)
T
**
12 Q
u ,u ,...u
? ?
=
? ?
x (4.9)
24
n
n
0
10
n
2
n
3
n
Q1
n
Q1n
QQ1
C
uu
28
u
u
u
C
uu
28
?
?
+ ?
+
? ?
???
++
? ?
??
??? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
??
?? ?
+?
??
? ?
??
? ?
b
M
(4.10)
with
nn
00Q1
C,u,C
?
+
and
Q1
u
?
+
implemented as Dirichlet boundary conditions. Neumann
boundary conditions can also be readily invoked, but are not presented here for brevity.
The second step of TCS for the 1D Burgers? equation can be as evaluated from
another Ax = b problem using
2
13
Q2 Q
Q1
C
10
8
CC
1
88
CC
1
88
C
01
8
?
??
??
?
?
?
??
+
??
?+
=
??
?+
?
??
A OOO (4.11)
T
n1 n1 n1
12 Q
u,u,..u
++ +
? ?
=
? ?
x (4.12)
25
0
11
22
33
Q1 Q1
Q1
QQ
1C
uD
28
1
uD
2
1
uD
2
1
uD
2
C
1
uD
28
?
??
??
??
??
??
?
+??
? ?
++
? ?
? ?
? ?
+
? ?
? ?
? ?
+
? ?
=
? ?
? ?
? ?
+
? ?
? ?
? ?
+?
? ?
? ?
b
M
(4.13)
Thus, the TCSF method for the 1D conservative Burgers? equation requires solution of
two tridiagonal linear problems at each time step.
If we further split the diffusion term, we obtain the TCSFD form for Equation
(4.1) as,
( )
n
2n
n
2
uttu
1uu
4x 2x
?
???
????
+=+
??
? ?
??
(4.14)
() ()
2
n1
2
u
t1
u
22 x x
?
+ ?
????
?
+??=
??
??
????
??
(4.15)
The first step is a matrix equation Ax=b with
n
2
nn
13
nn
Q2 Q
n
Q1
C
10
8
CC
1
88
CC
1
88
C
01
8
?
?
??
+
??
?+
=
??
?+
?
??
A OOO (4.16)
26
T
** *
12 Q
u,u,...u
? ?
=
? ?
x (4.17)
n
nn
0
11
nn
22
nn
33
nn
Q1 Q1
n
Q1nn
QQ
1C
uD
28
1
uD
2
1
uD
2
1
uD
2
C
1
uD
28
??
+
? ?
++
? ?
? ?
? ?
+
? ?
? ?
? ?
+
=? ?
? ?
? ?
? ?
+
? ?
? ?
? ?
+?
? ?
? ?
b
M
(4.18)
In the second step, we have
*
2
**
13
**
Q2 Q
*
Q1
C
10
28
CC
1
28 28
CC
1
28 28
C
01
28
?
?
??
?
+? ? +
??
??
?? +? ?+
=
??
??
?? +? ?+
?
?? +?
??
A OOO (4.19)
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?
=
? ?
x (4.20)
27
n
*n1
0
10
*
2
*
3
*
Q1
n
Q1*n1
QQ
C
uu
28
u
u
u
C
uu
28
+
?
+ +
+
? ?
???
++
? ?
??
??
? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ?
? ?
? ?
? ?
? ?
??
?
? ?
+?
??
? ?
??
? ?
b
M
(4.21)
Similar to TCSF, the TCSFD method for the 1D conservative Burgers? equation also
requires solution of two tridiagonal linear problems at each time step. However, TCSFD
reverses the solution process in TCSF. The first step in TCSFD is in the same structure
as the second step in TCSF; the second step in TCSFD is in the same structure as the
first step in TCSF.
To compare the above TCS methods to other temporally 2
nd
order accurate
implicit nonlinear solution methods, we apply CrankNicolson 2
nd
order temporal and
central difference spatial 2
nd
order discretization to Equation (4.1), resulting in
nn nn n n n
n1 n
i1 i1 i1 i1 i1 i i1
ii 2
n1 n1 n1 n1 n1 n1 n1
i1 i1 i1 i1 i1 i i1
2
tuu uu u 2u u
uu
24x x
uu uu u 2u u
4x x
+
++ ?? + ?
++ ++ + + +
++ ?? + ?
??? ?+
=? ??
?
??
?
?
+??
?
??
?
(4.22)
Conventional linearization methods such as Picard iteration, Newton iteration and local
linearization methods are used to solve Equation (4.22).
Picard Iteration
The simplest approach to implement for a nonlinear equation such as Equation
(4.22) is a laggedcoefficient iteration (Tannehill et al. 1997,pg 450) , which is a form of
Picard iteration. A linear inner equation is formed by estimating the time n+1 flux as
()
k
n1
i
u
+
, where the additional superscript ?k? is introduced as the outer iteration counter.
Equation (4.22) can then be represented as the Picard method
28
()
()()()()
() ()()
nn nn n n n
k1
n1 n
i1 i1 i1 i1 i1 i i1
ii 2
kk1 kk1
n1 n1 n1 n1
i1 i1 i1 i1
k1 k1 k1
n1 n1 n1
i1 i i1
2
tuu uu u 2u u
uu
24x x
uu uu
4x
u2u u
x
+
+
++ ?? + ?
++
++ ++
??
+++
+?
=?
??? ?+
??
?
??
?
?
+
?
?
?+
?
??
?
?
?
?
(4.23)
which is solved for k={1,2,3...} until an appropriate convergence criterion is reached. .
Thus, each outer iteration requires an inner linear solution of the form Ax = b. For
comparison with other methods, it is useful to define an
kk1 k+
= +Ax b c form in which
the b vector requires only a single computation in each time step whereas the c vector is
recomputed in each outer iteration. For Equation (4.23), the result is
k
3
kk
24
k
kk
Q2 Q
k
Q1
C
10
28
CC
1
28 28
CC
1
28 28
C
01
28
?
?
??
?
+? ? +
??
??
?? +? ?+
=
??
??
?? +? ?+
?
?? +?
??
A OOO (4.24)
()() ()
T
k1 k1 k1
k1 n1 n1 n1
12 Q
u , u ,... u
++ +
++ + +??
=
??
??
x (4.25)
29
nn n
11 1
nn n
22 2
nn n
33 3
nn n
Q1 Q1 Q1
nn n
QQ Q
11
uGD
82
11
uGD
82
11
uGD
82
11
uGD
82
11
uGD
82
???
? ?
?+
? ?
? ?
? ?
?+
? ?
? ?
? ??+
=
? ?
? ?
? ?
? ?
?+
? ?
? ?
? ?
?+
? ?
b
M
(4.26)
k
0
k
k
Q1
1
C
28
0
0
0
1
C
28
+
?
? ?
+
? ?
? ?
? ?
? ?
? ?
=
? ?
? ?
? ?
? ??
?
? ?
? ?
c
M
(4.27)
For the Picard iteration, the first outer iterative solution is started with
()
1
n1 n
ii
uu
+
= . The
outer iteration is stopped when
()()
{ }
k1 k
n1 n1
+
++
? 0.1.
However, the TCS methods remain stable in the whole range of CFL number.
10
?2
10
?1
10
0
10
1
10
?15
10
?14
10
?13
10
?12
10
?11
10
?10
10
?9
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
CFL
?
RMS
TCSF
TCSF?D
RK2
RK4
2nd order slope
unstable
unstable
Figure 4.12 Temporal accuracy of the TCS and RK methods for solution of 1D non
conservative Burgers? equation at different CFL numbers, where ?=0.05,
x? =1/50 and t? =0.6/? with ? ?{500, 200, 100, 80, 50, 30, 25, 20, 10, 8,
5}.
Computational requirement
The same assessment of operation counts is used in this nonconservative
Burgers? equation. All the methods perform similar to the nonconservative Burgers?
equation. The TCSF and TCSGD have the identical operation counts because they have
the same solution process only in a reverse mode. The same correspondence can be found
between the TCSG and TCSFD for the same reason. In this 1D nonconservative
Burgers? equation test case, the TCSF and TCSGD have the least operation counts.
59
1 2 3 4 5 6 7
0
100
200
300
400
500
600
R
operations per grid point
TCSF
TCSG
TCSF?D
TCSG?D
local linearization
Newton iteration
Picard iteration
Figure 4.13 Ideal operations per grid point for one time step using various nonlinear
solution methods for 1D Nonconservative Burgers equation. R is the
number of outer iterations taken by Newton method. It is assumed that the
Picard method convergences in
2
R outer iterations.
4.3 APPLICATION OF THE TCS METHOD TO A 1D NONLINEAR ORDINARY
DIFFERENTIAL EQUATION
The TCS method can be extended to any equation system with quadratic
nonlinearities. To test this general capability, we use a nonlinear ODE as our example:
dy
y(1 y) 0
dt
+ ?= (4.101)
For the initial condition y(0) 1/ 2= , Equation (4.101) has the analytical solution (Moin
2001)
[ ]
1
y 1 exp(t)
?
=+ (4.102)
4.3.1 The TCS discretizations
The TCSF Discretization
As in the 1D conservative Burgers? equation, the TCSF and TCSG are the same
for Equation (4.101). Applying the TCSF method to Equation (4.101) provides
60
()
nn**
t
yy yyy
2
?
?
=+ ? (4.103)
()
n1 * n1 * *
t
yyyyy
2
++
?
=+ ? (4.104)
The discrete form can be written from Equation (4.103) and (4.104) as
n
n
y
y
t
1y1
2
?
=
?
? ?? ?
? ?
(4.105)
**
n1
*
t
yy
2
y
t
1y
2
+
?
?
=
?
?
(4.106)
The TCSFD Discretization
Similarly, the TCSFD discrete format for Equation (4.101) is written as
nn
n
t
yy
2
y
t
1y
2
?
?
?
=
?
?
(4.107)
()
*
n1
*
y
y
t
1y1
2
+
=
?
? ?
(4.108)
It can be seen that Equation (4.108) has the same structure as Equation (4.105) and
Equation (4.107) has the same structure as Equation (4.106). Therefore, the TCSF and
TCSFD discrete forms for Equation (4.101) are similar but reversed.
4.3.2 Conventional linearization methods
CrankNicolson (CN) discretized Equation (4.101) is,
() ()
n1nnnn1n1
tt
yy y1y y1y
22
+++
? ?
=+ ? + ? (4.109)
To linearize the nonlinear term on right hand side of Equation (4.109), we use the
Newton iteration, Picard iteration and local linearization methods.
61
Newton method
Applying the Newton iteration approach, Equation (4.109) can be written as:
() ()
( )
k
n1
k1 k
n1 n1
n1
k
gy
yy
g
y
+
+
++
+
=?
???
??
?
??
(4.110)
where the outer superscript indicates the outer iteration number and
() () ()
n1 n1 n1 n1 n n n
tt
gy y y 1y y y 1y
22
++ + +
? ?
=? ? ?? ? (4.111)
Substituting Equation (4.111) into Equation (4.110), the original ODE in the Newton
linearized CN discretized form is
() ()
() () () ()
()
kkk
n1 n1 n1 n n n
k1 k
n1 n1
k
n1
tt
yy1yyy1y
22
yy
t
12y 1
2
+++
+
++
+
? ?
??
?????
??
??
=?
?
??
??
??
??
(4.112)
n1 n
0
yy
+
= can be used as the initial condition to start the solution.
Picard method
A second conventional approach, Picard iteration, in Equation (4.109), results in
() () () ()
k1 k k
n1 n n n n1 n1
tt
yyy1yy1y
22
+
+++
? ?
??
=+ ? + ?
??
??
(4.113)
The same initial condition,
n1 n
0
yy
+
= can be used to start the solution.
Local linearization
Applying the third conventional method, local linearization method, Equation
(4.109) can be written as
n
n1 n
n
tf
yy
tf
1
2y
+
?
=+
? ?
????
?
? ?
??
?
??? ?
? ?
(4.114)
where
( )fy1y=? (4.115)
62
Substituting Equation (4.115) into Equation (4.114) provides the locallylinearized
Equation (4.109) as
( )
()
nn
n1 n
n
ty 1y
yy
t
12y1
2
+
??
=+
?
? ?
??
? ?
? ?
(4.116)
In this specific case, the local linearization, TCSF and TCSFD are mathematically
equivalent. If we substitute Equation (4.105) to Equation (4.106) and substitute Equation
(4.107) to Equation (4.108), we will have two expressions which are exactly the same as
Equation (4.116).
4.3.3 Comparisons between the TCS and other methods
The evolution of the absolute errors for two TCS methods applied to the ODE is
shown in Figure 4.14. The maximum absolute errors for both TCS discrete formats occur
at time t=1.5.
0 2 4 6 8 10
0
1
x 10
?4
t
?
abs
TCSF
TCSF?D
Figure 4.14 Time evolution of absolute errors for TCS methods applied to the ODE,
Equation (4.101), for t? =0.1, where
abs model analytical
(t) y (t) y (t)?= ? .
Figure 4.15 provides a comparison of rootmeansquare (RMS) error over ? time
steps prior to maximum error, defined as
63
1/2
2
RMS mod el i analytical i
i1
1
y (t) y (t)
?
=
??
??
???? ?
??
??
?
??
?
(4.117)
where
max
t1.5= and { }10, 50, 100, 200, 500, 600, 800, 1000, 1200, 1500, 2000?? . All
tested methods provide secondorder accuracy. In this particular case, the iterative
methods have slightly smaller error than noniterative methods. The errors from the three
noniterative methods are identical since TCSF, TCSFD and local linearization are
mathematically equivalent.
10
1
10
2
10
3
10
?9
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
?
?
RMS
TCSF
TCSF?D
Picard
Newton
Local linearization
2nd order
Figure 4.15 RMS error from Equation (4.117) for discrete solutions of the ODE,
Equation (4.101), computed using a range of ? time steps.
4.4 SUMMARY
The TCS method has been tested on three different test cases: a 1D conservative
Burgers? equation, a 1D nonconservative Burgers? equation and an ODE. The TCS
method has multiple discrete forms for each test case. All the TCS discretizations have
been compared to conventional implicit linearization methods, including Picard iteration,
Newton iteration and local linearization methods. All the TCS discretizations
64
demonstrate 2
nd
order temporal accuracy and are shown stable up to a CFL O (10). All
the TCS discretizations are proved stable when the CFL value is in the order of 10. The
stability advantage of the TCS methods has been further demonstrated by the comparison
with the RK methods. The TCS methods require fewer operation counts than Picard and
Newton method, but are similar to local linearization. The principal advantage of the
TCS method over local linearization is the relative ease with which the TCS method can
be derived and implemented as it does not require discrete evaluation of a function
Jacobian. The TCS method is applicable to any quadraticallynonlinear problems, which
is verified by the ODE case.
For each test case, different TCS discretizations perform similarly in temporal
accuracy, stability and efficiency. However, difference can be observed among different
discretizations. In the 1D conservative Burgers? equation, the TCSF and TCSG collapse
into one single form. The TCSFD reverses the solution process of the TCSF.
Nevertheless, the TCSFD has a noticeably smaller relative error than the TCSF and other
methods. In the 1D nonconservative Burgers? equation, the TCSF and TCSG hold
distinct forms. The TCSFD reverses the solution process of the TCSG, and the TCSGD
reverses the solution process of the TCSF. Among the four discretizations, the TCSF has
the smallest relative error, which is about one order of magnitude smaller than the TCSG
D. In the 1D nonlinear ODE case, the TCSF and TCSG remain the same. The TCSFD
reverses the solution process of the TCSF, but they have the same relative error.
65
Chapter 5 The TCS Family Method
In Chapter 3, we revealed that the key to the TCS method is the 2ndorder
accurate timecentered splitting. Applying this splitting to different terms, different TCS
discretizations can be generated. In Chapter 4, we compared the TCS discretizations with
conventional linearization methods using various examples. In this chapter, we derive the
TCS family method in one general form by introducing weighting factors to different
terms. A significance of this TCS family method is that the value of the weighting factors
does not affect the order of accuracy of the scheme. Since weighting factors can be any
value between zero and one while still providing essentially the same advantages, it may
provide a flexible approach for a variety of problems. To examine the properties of
weighting factors, the 1D nonconservative Burgers? equation is solved using different
combinations of weighting factors. Furthermore, we demonstrate that the TCS method
can computationally decouple an equation system of coupled variables using special
combinations of weighting factors. This advantage makes the TCS method more efficient
in solving coupled multivariable equation systems.
5.1 DERIVATION OF THE TCS FAMILY METHOD
Let us start a time marching differential equation from a simplest form:
()() ()
{}
fFM L N
t
??
= ?? ???
??
?
, (5.1)
where ? s a generic variable, F, M, L and N represent generic linear operators. The
function f is a linear function of ( ) ( )FM L? ?? ?
? ?
and ( )N ? . Equation (5.1) then
represents a time marching differential equation with a quadratic nonlinearity of
()()FM L????
??
and a linear operator of the variable, ( )N ? . This quadratic nonlinearity
can be: the product of the variable itself,
2
? if F, M and L are equal to 1; the product of
the variable and the gradient of the variable, / x??? ? , if F=1, M=1 and L= ()/x??; the
gradient of the product of the variable, ( )/x? ?? ? if M=1, L=1 and F= ()/x??. Thus,
Equation (5.1) represents a generic time marching differential equation with any kind of
quadratic nonlinearity. Equation (5.1) can be discretized using the midpoint rule, written
as
( ) ( ) ( )
{ }
( )
n1 n n12 n12 n12 3
fFM L N t t
++++
??
? =? + ?? ??+ ?
??
// /
,O (5.2)
In the previous chapters, we demonstrated that the timecentered split can be introduced
to the quadratic nonlinear term and/or the linear term. Introducing two weighting factors,
66
1
? and
2
? , and combining with the timecentered linear approximation,
()
nn1
n12 2
t
2
/
O
+
+
?+?
?= +? , Equation (5.2) can be written as a general discretization
()()()
()() ()
nn1 nn1
n1 n n12 n12
11
nn1
n12 3
22
fFM L 1 FM L
N1N tt
2
++
+++
+
+
?
????
?? ???+? ?+?
?
?=?+? ? +?? ?
? ????
?? ??
?? ???
????
?
?
???+?
?? +?? ?+?
?
??
??
?
//
/
,
O
(5.3)
where the weighting factor
i
? ( { }i12,? ) is required that
i
01??? . Equation (5.3) is
mathematically equivalent to Equation (5.2). Similar to the derivation in Chapter 3, an
intermediate variable
*
? is introduced as
()( ) ()( )()
{
()()()}
nnn12 n12n
11
n12 n
22
fFM L 1 FM L
t
N1N
2
++
+
????
?=?+ ? ? ? + ?? ? ?
????
?
?? +?? ?
*//
/
,
(5.4)
Subtracting Equation (5.4) from Equation (5.3), the second step is written as:
()( ) ()( )( )
{
()()()} ()
n1 n1 n12 n12 n1
11
n12 n1 3
22
fFM L 1 FM L
t
N1N t
2
+++ ++
++
????
?=?+? ? ? +?? ? ?
????
?
?? +?? ? +?
*//
/
,
O
(5.5)
Expanding
n12/+
? in Equation (5.4) using the Taylor expansion,
()
() ()
()()
2
nnn
1t
2
nn
1t
2
2t 2
tt
fFM L
22
tt
1FM L
22
tt t
N1N
22 2
?
??
??
???
??
?=?+ ? ? ?+? +??
???
??
??
? ??
???
??
??
??
+?? ?+? + ?
??
??
??
??
??
?
??
? ???
??
??+?+ +?? ?
???
??
??
???
?
*
O
O,
O
(5.6)
Organizing Equation (5.6) and grouping the higher order terms,
()()
{
()}
2
nnnn
tt
fFM L N
22
??
??
??
?=?+ ? ? ? +
??
??
??
*
,O (5.7)
n12/+
? can be approximated using the explicit Euler approximation as
67
()() ()
{}
2
n12 n n n n
tt
fFM L N
22
+
? ?
??
??
?=?+ ?? ? +
??
??
??
/
,O (5.8)
Therefore, we obtain the correspondence between
*
? and
n12/+
? from Equations (5.7) and
(5.8) such that
( )
n12 2
t
*/
O
+
? =? + ? (5.9)
Substituting Equation (5.9) into Equations (5.4) and (5.5), the original equation set
becomes a computationally linearized twostep equation system:
( ) ( ) ()( ) ( )
{
()()()}
nn n
11
n
22
fFM L 1 FM L
t
N1N
2
????
?=?+ ? ? ? + ?? ? ?
????
?
??+?? ?
***
*
,
(5.10)
( ) ( ) ()( ) ( )
{
()()()} ()
n1 n1 n1
11
n1 3
22
fFM L 1 FM L
t
N1N t
2
++ +
+
????
?=?+? ? ?+?? ? ?
????
?
??+?? ? +?
***
*
,
O
(5.11)
The weighting factor,
i
? , distinguishes this derivation above from the original derivation
of TCS method in Chapter 3 because the value of
i
? will not affect the 2
nd
order accuracy
of the TCS methods. In addition, there is no dependent relationship between
1
? and
2
? . In
other words, as long as
i
01??? , Equations (5.10) and (5.11) are a 2
nd
order equivalent
to Equation (5.2). Although the derivation above uses a single variable, the same
mathematic principles apply to any quadratic nonlinearity. In ordinary differential
equations (ODE), the nonlinearity is simply
2
x or xy; in partial differential equations
(PDE), this quadratic nonlinearity can be UU/ x? ? or ( )HU / x? ? . Therefore, the
derivation above is applicable to an ODE or ODE system with quadratic nonlinearity and
typical flow transport equation such as 1D advectiondiffusion equation or shallow water
equations. In addition, the more variables are involved in an equation system, the more
weighting factors that can be introduced. As a result, the TCS method provides not a
single kind, but a whole family of TCS discretizations.
As we illustrated in the previous chapters, all the TCS family methods keep 2
nd

order temporal accuracy and they share the same simplicity advantages in computation.
However, with different combinations of
i
? , the numerical solutions perform differently.
In the next section, we will use the 1D nonconservative Burgers? equation as an example
to explore the different performance of different TCS discretizations.
68
5.2 APPLICATION OF THE TCS FAMILY METHOD TO THE 1D NONCONSERVATIVE
BURGERS? EQUATION
With the introduction of the weighting factors
1
? and
2
? , the general format of
the TCS discretized nonconservative Burgers? equation is written as,
()()()()(){}
nn n2 2n
1x 1 x 2x 2x
t
uu uu 1 uu u 1 u
2
****
?
=+ ??? +???? +??? +???? (5.12)
()()() ()
n1 n1 n1 2 2 n1
1x 1x 2x 2x
t
uu uu1 uu u1 u
2
****++ + +
?
=+ ??? +???? +???+???? (5.13)
The above equations are a twostep matrix equation system in the format of Ax=b. Using
a central difference in spatial derivatives, the LHS coefficient matrix A in the first step is
tridiagonal such that
1,1 1,2
2,1 2,2 2,3
3,2 3,3 3,4
Q2,N3 Q2,Q2 Q2,Q1
Q1,N 2 Q1,Q1 Q1,Q
Q,Q 1 Q,Q
aa 0
aaa
aaa
aaa
aa a
0a
?? ?? ??
?? ?? ?
?
??
??
=
??
A OOO (5.14)
where Q is the number of the grid points in space, and the nonzero components of A for
{ }i 1,2...Q= are of the form
n
12
i
nn
i,j 1 i 1 i 1 2
n
12
i
C:ji1
42
1
a1 CC :ji
4
C:1
42
+?
??
?
??? =?
?
?
?
??=+? ? +?? =
?
??
?
??
?
+?? =+
?
?
(5.15)
Also in the first step,
T
** *
12 Q
u ,u ,...u? ?=
? ?
x (5.16)
and
69
( )
()
()
()
n
n
n2111 2
1
n
n22
2
n
2Q1n
Q1
nn
2Q1Qn
2
Q
1D
C
u
24
1D
u
2
1D
u
2
1DC
u
24
?
?
???? ?
? ??
+++
??
?? ?
+
=
??
?? ?
+
?? ? ?
? ?
??
+?+
??
b M (5.17)
In the second step, we also have a tridiagonal coefficient matrix A with different
tridiagonal elements as
() ()
*
12
i
**
i,j 1 i 1 i 1 2
*
12
i
11
C:ji1
42
1
a11 CC 1 :ji
4
11
C:1
42
+?
?? ???
??? =?
?
?
?
??=+ ?? ? +??? =
?
??
?
?? ??
?
+?? =+
?
?
(5.18)
and
T
n1 n1 n1
12 Q
u ,u ,...u
++ +
? ?=
? ?
x (5.19)
( ) ( )
()()
*
*
*112
21
1
*
*
22
2
*
2Q1*
Q1
**
2Q 1Q*2
Q
1C1
D
u
24
D
u
2
D
u
2
D1 C1
u
24 2
?
?
???????
??
++ +
??
??
+
=
??
??
+
?? ?? ???
+? +
??
b M (5.20)
where C,? and D in these two steps have the same definitions as in Chapter 4.
70
5.3 RESULTS AND DISCUSSION
To explore the properties of different weighting factors, the 1D nonconservative
Burgers? equation is solved using Equations (5.12) and (5.13) with different
combinations of
1
? and
2
? . Before we analyze the results from various combinations of
1
? and
2
? , it is necessary to review their meanings. In the TCS scheme,
n12
u
/+
in the
Burgers? equation is either approximated using the timecentered splitting as
()()
n12 n n1 2
uuu2t
/
/O
++
=+ +? or u
*
as
( )
n12 2
uu t
/*
O
+
= +?. Different
i
? indicates
different approximations of
n12
u
/+
. When
i
? =0 or
i
? =1,
n1/2
u
+
is approximated only
using one kind of the approximations. When
i
? lies between 0 and 1,
n1/2
u
+
is
approximated by the combinations of these two approximations. Table 5.1 lists the
mathematical meaning of different values of
i
? .
Table 5.1 The mathematical meaning of the value of
1
? and
2
? .
value
1
?
2
?
0
split the entire gradient term and
approximate the entire flux term using
u
*
split the entire diffusion
term
0and1<<
split part of the gradient and flux terms
approximate part of the gradient and flux
terms using u
*
split part of the diffusion
term and approximate part
of the diffusion term using
u
*
1
split the entire flux term and
approximate the entire gradient term
using u
*
approximate the entire
diffusion term using u
*
5.3.1 Accuracy of Different Weighting Factors
The error analysis of different TCS discretizations is evaluated in the same vein as
in Chapter 4. In figure 5.1, the RMS error is plotted against a range of CFL numbers for
six different
1
? values ( { }
1
0131223341,/,/, /,/,?? ). For each
1
? value, four
different
2
? values ( { }
2
013121,/,/,?? ) are chosen for comparison, resulting
twenty four different TCS discretizations.
The total marching time t=0.6 is used as the basis for error comparison in Figure
5.1 because the maximum absolute error occurs at t=0.6 (as demonstrated in Chapter 4).
The CFL number is defined using
analytical
u, x150/? = and t? = 06./? , where ? is the
71
number of the time steps and { }5, 8, 10, 20, 25, 30, 50?? . Thus, the tested CFL number
covers the range of 04 CFL 42..??. It can be observed that every TCS discretization
has the same 2ndorder temporal accuracy and remains stable in the whole range of the
CFL numbers. The computation requirement of each TCS discretization is essentially the
same since they share the same general structure of matrix equations as in Equations
(5.14) to (5.20). Figure 5.1 clearly shows that for a given value of
1
? , any choice of
2
?
gives an error of the same magnitude. The reverse also holds true, as shown in Figure 5.2.
The RMS error is plotted for six values of
2
? , choosing four
1
? values for each
2
? value.
Again, given a value for
2
? , the choice of
1
? affects the RMS error less than or equal to
one order of magnitude.
Though
12
23 12/, /?= ?= gives the lowest RMS error in these two examples,
the best combination for the smallest error is problemspecific; error arises from
approximations of different terms. To obtain optimal combinations of
i
? for the lowest
error, one needs to examine the nature of the solution and each term at the designated
marching time. For this Burgers? equation, the solution u, the gradient term ux/??, and
the diffusion term
22
ux/?? ? all need to be taken into consideration because
approximations have been introduced to all of these terms.
72
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
1
=1
?
2
=0 ?
2
=1/3 ?
2
=1/2 ?
2
=1
2nd order slope
Figure 5.1 Temporal accuracy of various combinations of
1
? and
2
? for solution of the
Burgers? equation at t=0.6, where x150/? = , 005.?= and t06./?= ?
( { }5, 8, 10, 20, 25, 30, 50?? ).
73
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=0
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=1/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=1/2
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=2/3
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=3/4
10
?1
10
0
10
1
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
CFL
?
RM
S
?
2
=1
?
1
=0 ?
1
=1/3 ?
1
=1/2 ?
1
=1
2nd order slope
Figure 5.2 Temporal accuracy of various combinations of
2
? and
1
? for solution of the
Burgers? equation at t=0.6, where x150/? = , 005.?= and t06./?= ?
( { }5, 8, 10, 20, 25, 30, 50?? ).
74
5.3.2 Stability of Different Weighting Factors
Figure 5.1 and 5.2 display the temporal accuracy of the solution of the Burgers?
equation at t=0.6. All TCS discretizations produce stable solutions for the whole range of
tested CFL numbers. To further explore the stability characteristics, we investigate the
performance of the TCS family method at t=0.3, 1 and 2, and these results are shown in
Figure 5.3, 5.4 and 5.5, respectively. Similar to Figure 5.1, the RMS error is plotted
against a range of CFL numbers for several choices of
2
? given a
1
? value. Similar to
t=0.6, all 24 TCS discretizations stay stable at t=0.3.
At t=1and t=2, though, we find different results. At t=1 (Figure 5.4), for any given
1
? value, the solution of
2
1?=becomes unstable with the increase of the CFL number,
but the solutions for any
2
1?? remain stable. At t=2, the solutions of both
2
0?= and
2
1?=become unstable with the increase of the CFL number, but the solutions with any
2
01