Resumen de High resolution schemes for hyperbolic conservation laws with source terms
The laws establishing the conservation of mass, momentum and energy in a physical system translate into a well defined system of partial differential equations, In these equations, the effect of sources, sinks, chemical reactions and other phenomena of interest are modeled by the inclusion of additional terms, which are generically referred to as source terms.
This memoir is devoted to the study of the numerical treatment of source terms in hyperbolic conservation laws and systems. In particular, we study two types of situations that are particularly delicate from the point of view of their numerical approximation: The case of balance laws, with the shallow water system as the main example, and the case of hyperbolic equations with stiff source terms.
There are nowadays many techniques that can produce accurate numerical solutions of homogeneous conservation laws and systems. It is well known that these solutions can be discontinuous, even when the initial data in a Cauchy problem is perfectly smooth. Standard (linear, data-independent) finite-difference, finite-volume or finite-element techniques tend to produce oscillatory numerical approximations when the order of accuracy of the scheme is larger than one. A large amount of research in the last decades has resulted in a well established technology to construct High Resolution Shock Capturing (HRSC henceforth) schemes. These schemes lead to accurate results away from discontinuities, as well as well defined, very steep, monotone profiles at those locations where discontinuities in the solution do occur.
On the other hand, the numerical solution of Ordinary Differential Equations (ODE) is a well established discipline that has produced a variety of numerical techniques that have proved to be useful in many applications.
Because of these facts, a standard approach to solve hyperbolic conservation laws with source terms is to apply the so-called fractional step approach. This technique alternates between solving a homogeneous conservation law and solving an ODE that contains only the source term. However, there are situations where the fractional-step approach does not lead to acceptable numerical approximations.
When computing numerical approximations to balance laws, such as the shallow water equations, in steady-state or quasi-steady-state situations, the numerical solutions of the homogeneous PDE and the ODE are required to balance exactly. This exact balance is not likely to be respected by the fractional splitting procedure, and parasitic waves of a purely numerical nature can occur.
For stiff source terms, the usage of a stiff ODE solver combined with a HRSC method can lead to numerical solutions that look reasonable but are completely wrong. This phenomenon was observed as early as 1986 by Colella, Majda and Roytburd in [19] on a model combustion problem that involved the Euler equations of Gas Dynamics coupled with a single chemistry variable representing the mass fraction of unburnt gas in a detonation wave. The structure of the detonation waves obtained was well understood, and it was observed that the numerical solution obtained was qualitatively incorrect when computing on coarse grids.
Stiff source terms could describe also the effect of relaxation as in the kinetic theory of rarefied gases, hydrodynamical models for semiconductors elasticity with memory, water waves, traffic flows, etc.
There has been a large amount of literature in recent years devoted to the numerical problems that occur in these two types of situations.
For shallow water equations, several authors extended the classical Riemann solver of Roe to nonhomogeneous problems related to balance laws [5], [7], [17], [58], [104]. In these works, the discrete form of the source terms is constructed in a way similar to that employed for the construction of the numerical fluxes, seeking an equilibria that exists in a steady-state conservation law with source terms. The idea of source-term upwinding lead Bermúdez and Vázquez-Cendón [5] to formulate the so-called C-property (for Conservation property) for a numerical scheme, which prevents the propagation of parasitic waves in steady and quasi-steady flows. Independently, Greenberg and Leroux [46] coined the term well-balanced for schemes that preserve steady states at the discrete level. These ideas have been explored and developed for shallow water flows in the recent literature [4], [24], [42], [61], [71], [78], [106] ...
In this work, we follow the strategy described by Gascón and Corberán in [38] and Donat, Caselles and Haro in [10]. In [38], the authors propose to formally write the source term in divergence form so that the nonhomogeneous problem can be 'transformed' into a 'homogeneous form' through the definition of a new flux function. This change seeks to preserve the balance of the source and flux terms at steady states in an almost automatic manner, and suggests a way to apply well known schemes for homogeneous conservation laws to the non-homogenous case. However, as they readily observe, the application of the numerical methods for the homogeneous case is not immediate and adequate formalizations are required.
In [10], the idea of flux gradient and source term balancing in [38] was incorporated into the numerical scheme developed by Donat and Marquina in [26], thus effectively extending this scheme to balance laws.
In this work, we concentrate on the theoretical foundations of high-resolution total variation diminishing (TVD) schemes for homogeneous scalar conservation laws, firmly established through the work of Harten [50], Sweby [95], and Roe [80] and analyzed the properties of a second order, flux-limited version of the Lax-Wendroff scheme which avoids oscillations around discontinuities, while preserving steady states [38]. When applied to homogeneous conservation laws, TVD schemes prevent an increase in the total variation of the numerical solution, hence guaranteeing the absence of numerically generated oscillations. They are successfully implemented in the form of flux-limiters or slope limiters for scalar conservation laws and systems. Our technique is based on a flux limiting procedure applied only to those terms related to the physical flow derivative/Jacobian.
With respect to the numerical treatment of stiff source terms, we follow Leveque and Yee in [73]. Taking the simple model problem considered in [73], we study the properties of the numerical solution obtained with different numerical techniques. We are able to identify the delay factor, which is responsible for the anomalous speed of propagation of the numerical solution on coarse grids. The delay is due to the introduction of non-equilibrium values through numerical dissipation, and can only be controlled by adequately reducing the spatial resolution of the simulation. Explicit schemes suffer from the same numerical pathology, even after reducing the time step so that the stability requirements imposed by the fastest scales are satisfied. We study the behavior of Implicit-Explicit (IMEX) numerical techniques, as a tool to obtain high resolution simulations that incorporate the stiff source term in an implicit, systematic, manner. The IMEX framework has been also successfully applied to hyperbolic systems with relaxation (see [9], [73], [110]).
Usually, when using very fine uniform grids, we find that the computational time becomes the main drawback in the numerical simulation. For some high resolution shock capturing schemes (HRSC), fine mesh simulations in two dimensions are out of reach simply because they cost too much. The numerical flux evaluations are too expensive. However, the flux computations are needed only because nonsmooth structures may develop spontaneously in the solution of a hyperbolic system of conservation laws and evolve in time, which lead to develop techniques that reduce the computational effort associated to these simulations. Harten in [49] proposed a scheme based on reducing the computational cost using the smoothness information of the data, and replacing the expensive numerical flux with a cheap polynomial interpolation in the smooth regions. The key is the use of different multilevel strategy to reduce the computational effort associated to HRSC scheme. In smooth regions, Harten in [49] proposes to evaluate the numerical flux function of the HRSC only on a coarse grid and to use these values to compute the fluxes on the finest grid using an inexpensive polynomial interpolation process in a multilevel fashion. Here, we extend the technique developed in [16] to hyperbolic conservation laws with source terms and apply the multilevel technique to the shallow water system.
