Riemann problems are of prime importance in computational fluid dynamics simulations using finite elements or finite volumes discretizations. In some applications, billions of Riemann problems might need to be solved in a single simulation, therefore it is important to have reliable and computationally efficient algorithms to do so. Given the nonlinearity of the flux function in most systems considered in practice, to obtain an exact solution for the Riemann problem explicitly is often not possible, and iterative solvers are required. However, because of issues found with existing iterative solvers like lack of convergence and high computational cost, their use is avoided and approximate solvers are preferred. In this thesis work, motivated by the advances in computer hardware and algorithms in the last years, we revisit the possibility of using iterative solvers to compute the exact solution for Riemann problems. In particular, we focus on the development, implementation, and performance comparison of iterative Riemann solvers for the shallow water and Euler equations. In a one-dimensional homogeneous framework for these systems, we consider several initial guesses and iterative methods for the computation of the Riemann solution. We find that efficient and reliable iterative solvers can be obtained by using recent estimates on the Riemann solution to modify and combine well-known methods. Finally, we consider the application of these solvers in finite volume simulations using the wave propagation algorithms implemented in Clawpack.

We develop and analyse explicit in time schemes for the computation of compressible flows, based on staggered in space. Upwinding is performed equation by equation only with respect to the velocity. The pressure gradient is built as the transpose of the natural divergence. For the barotropic Euler equations, the velocity convection is built to obtain a discrete kinetic energy balance, with residual terms which are non-negative under a CFL condition. We then show that, in 1D, if a sequence of discrete solutions converges to some limit, then this limit is the weak entropy solution. For the full Euler equations, we choose to solve the internal energy balance since a discretization of the total energy is rather unnatural on staggered meshes. Under CFL-like conditions, the density and internal energy are kept positive, and the total energy cannot grow. To obtain correct weak solutions with shocks satisfying the Rankine-Hugoniot conditions, we establish a kinetic energy identity at the discrete level, then choose the source term of the internal energy equation to recover the total energy balance at the limit. More precisely speaking, we prove that in 1D, if we assume the L∞ and BV-stability and the convergence of the scheme, passing to the limit in the discrete kinetic and internal energy equations, we show that the limit of the sequence of solutions is a weak solution. Finally, we consider the computation of radial flows, governed by Euler equations in axisymetrical (2D) or spherical (3D) coordinates, and obtain similar results to the previous sections. In all chapters, we show numerical tests to illustrate for theoretical results. / We develop and analyse explicit in time schemes for the computation of compressible flows, based on staggered in space. Upwinding is performed equation by equation only with respect to the velocity. The pressure gradient is built as the transpose of the natural divergence. For the barotropic Euler equations, the velocity convection is built to obtain a discrete kinetic energy balance, with residual terms which are non-negative under a CFL condition. We then show that, in 1D, if a sequence of discrete solutions converges to some limit, then this limit is the weak entropy solution. For the full Euler equations, we choose to solve the internal energy balance since a discretization of the total energy is rather unnatural on staggered meshes. Under CFL-like conditions, the density and internal energy are kept positive, and the total energy cannot grow. To obtain correct weak solutions with shocks satisfying the Rankine-Hugoniot conditions, we establish a kinetic energy identity at the discrete level, then choose the source term of the internal energy equation to recover the total energy balance at the limit. More precisely speaking, we prove that in 1D, if we assume the L∞ and BV-stability and the convergence of the scheme, passing to the limit in the discrete kinetic and internal energy equations, we show that the limit of the sequence of solutions is a weak solution. Finally, we consider the computation of radial flows, governed by Euler equations in axisymetrical (2D) or spherical (3D) coordinates, and obtain similar results to the previous sections. In all chapters, we show numerical tests to illustrate for theoretical results.

Many applications in engineering practice can be described by thehyperbolic partial differential equations (PDEs). Numerical modeling of this typeof equations often involves large gradients or shocks, which makes it achallenging task for conventional numerical methods to accurately simulate suchsystems. Thus developing accurate and efficient shock capturing numericalschemes becomes important for the study of hyperbolic equations.In this dissertation, a detailed study of the numerical methods for linearand nonlinear unsteady hyperbolic equations was carried out. A new finitedifference shock capturing scheme of finite volume style was developed. Thisscheme is based on the high order Pad?? type compact central finite differencemethod with the weighted essentially non-oscillatory (WENO) reconstruction toeliminate non-physical oscillations near the discontinuities while maintain stablesolution in the smooth areas. The unconditionally stable semi-implicit Crank-Nicolson (CN) scheme is used for time integration.The theoretical development was conducted based on one-dimensionalhomogeneous scalar equation and system equations. Discussions were alsoextended to include source terms and to deal with problems of higher dimension.For the treatment of source terms, Strang splitting was used. For multidimensionalequations, the ?? -form Douglas-Gunn alternating direction implicit(ADI) method was employed. To compare the performance of the scheme withENO type interpolation, the current numerical framework was also applied usingENO reconstruction.The numerical schemes were tested on 1-D and 2-D benchmark problems,as well as published experimental results. The simulated results show thecapability of the proposed scheme to resolve discontinuities while maintainingaccuracy in smooth regions. Comparisons with the experimental results validatethe method for dam break problems. It is concluded that the proposed scheme isa useful tool for solving hyperbolic equations in general, and from engineeringapplication perspective it provides a new way of modeling open channel flows.

A finite volume, shock-capturing scheme is used to solve the shallow water equations on unstructured triangular meshes. The conditions are characterised by: slow flow velocities (up to 1m/s), long time scale (around 10 days), and large domains (50-100km across). Systematic verification is carried out by comparing numerical with analytical results, and by comparing parameter variation in the numerical scheme with perturbation analysis, and good agreement is found. It is the first time a shock-capturing scheme has been applied to slow flows in Moreton Bay. The scheme is used to simulate transport of a pollutant in Moreton Bay, to the east of the city of Brisbane, Australia. Tidal effects are simulated using a sinusoidal time-dependent boundary condition. An advection equation is solved to model the path of a contaminant that is released in the bay, and the effect of tide and wind on the contaminant is studied. Calibration is done by comparing numerical results with measurements made at a study site in Moreton Bay. It is found that variation in the wind speed and bed friction coefficients changes the solution in the way predicted by the asymptotics. These results vary according to the shape of the bathymetry of the domain: in shallower areas, flow is more subject to shear and hence changes in wind speed or bed friction had a greater effect in adding energy to the system. The results also show that the time-dependent boundary condition reproduces the tidal effects that are found on the Queensland coast, i.e. semi-diurnal with amplitude of about 1 metre, to a reasonable degree. It is also found that the simulated path of a pollutant agrees with field measurements. The computer model means different wind speeds and directions can be tested which allows management decisions to be made about which conditions have the least damaging effect on the area.

This thesis examines chaotic mixing in wavy-type channels and two-layer shallow water flow. For wavy-type channels, the equations of motion for vortices and fluid particles are derived assuming two-dimensional irrotational, incompressible flow. Instantaneous positions of the vortices and particles are determined using Lagrangian tracking, and are conformally mapped to the physical domain. Unsteady vortex motion is analysed, and vortex-induced chaotic mixing in the channels studied. The dynamics of mixing associated with the evolution of the separation bubble, and the invariant manifolds are examined. Mixing efficiencies of the different channel configurations are compared statistically. Fractal enhancement of productivity is identified in the study of auto-catalytic reaction in the wavy channel. For the two-layer shallow water model, an entropy-correction free Roe type two-layer shallow water solver is developed for a hyperbolic system with non-conservative products and source terms. The scheme is well balanced and satisfies the C-property such that smooth steady solutions are second order accurate. Numerical treatment of the wet-dry front of both layers and the loss of hyperbolicity are incorporated. The solver is tested rigorously on a number of 1D and 2D benchmark test cases. For 2D implementation, a dynamically adaptive quadtree grid generation system is adopted, giving results which are in excellent agreement with those on regular grids at a much lower cost. It is also shown that algebraic balancing cannot be applied directly to a two-layer shallow water flow due to the lack of simultaneous referencing for the still water position for both layers. The adaptive two-layer shallow water solver is applied successfully to flow in an idealised tidal channel and to tidal-driven flow in Tampa Bay, Florida. Finally, chaotic advection and particle mixing is studied for wind-induced recirculation in two-layer shallow water basins, as well as Tampa Bay, Florida.

Hydraulic structures such as dykes and groynes are commonly used to help control river flows and reduce flood risk. The present research aims to develop an idealized model of the hydrodynamics in the vicinity of a large river bend, and the advection and mixing processes where groynes are located. In this study a curvilinear model of shallow water equations is applied to investigate chaotic advection of particles in a river bend similar in dimensions to a typical bend in the River Danube, Hungary. First, a curvilinear grid generator is developed based on Poisson-type elliptic partial differential equations. The grid generator is verified for benchmark tests concerning a circular domain and for distorted grids in a rectangular domain. It is found that multi-grid (MG) and conjugate gradient (CG) methods performed better computationally than successive over-relaxation (SOR) in generating the curvilinear grids. The open channel hydrodynamics are modelled using the shallow water equations (SWEs) derived by depth-averaging the continuity and Navier-Stokes momentum equations. Both Cartesian and curvilinear forms of the shallow water equations are presented. Both sets of equations are discretized spatially using finite differences and the solution marched forward in time using fourth-order Runge-Kutta scheme. The shallow water solvers are verified and validated for uniform flow in the rectangular channel, wind-induced set up in rectangular and circular basins, flow past a sidewall expansion, and Shallow flow in a rectangular channel with single groyne. A Lagrangian particle tracking model is used to predict the trajectories of tracer particles, and bilinear interpolation is used to provide a representation of the continuous flow field from discrete results. The particle tracking model is verified for trajectories in the flow field of a single free vortex and in the alternating flow field of a pair of blinking vortices. Excellent agreement is obtained with analytical solutions, previously published results in the literature. The combined shallow flow and Lagrangian particle tracking model is then used to simulate particle advection in the flow past a side-wall cavity containing a groyne and reasonable agreement is obtained with published experimental and alternative numerical data. Finally, the combined model is applied to simulate the shallow flow hydrodynamics, advection and mixing processes in the vicinity of groynes in river bend, the dimensions representative of a typical bend in the Danube River, Hungary.

Flux gradient terms and source terms are two fundamental components of hyperbolic systems of balance law. Though having distinct mathematical natures, they form and maintain an exact balance in a special class of solutions, which are called steady-state solutions. In this dissertation, we are interested in the construction of well-balanced schemes, which are the numerical methods for hyperbolic systems of balance laws that are capable of exactly preserving steady-state solutions on the discrete level. We first introduce a well-balanced scheme for the Euler equations of gas dynamics with gravitation. The well-balanced property of the designed scheme hinges on a reconstruction process applied to equilibrium variables---the quantities that stay constant at steady states. In addition, the amount of numerical viscosity is reduced in the areas where the flow is in (near) steady-state regime, so that the numerical solutions under consideration can be evolved in a well-balanced manner. We then consider the shallow water equations with friction terms, which become very stiff when the water height is close to zero. The stiffness in the friction terms introduces additional difficulty for designing an efficient well-balanced scheme. If treated explicitly, the stiff friction terms impose a severe restriction on the time step. On the other hand, a straightforward (semi-) implicit treatment of the stiff friction terms can greatly enhance the efficiency, but will break the well-balanced property of the resulting scheme. To this end, we develop a new semi-implicit Runge-Kutta time integration method that is capable of maintaining the well-balanced property under the time step restriction determined exclusively by non-stiff components in the underlying equations. The well-balanced property of our schemes are tested and verified by extensive numerical simulations, and notably, the obtained numerical results clearly indicate that the well-balanced property plays an important role in achieving high resolutions when a coarse grid is used. / acase@tulane.edu

Waves and circulation processes interact in daily wind and tide driven flows as well as in more extreme events such as hurricanes. Currents and water levels affect wave propagation and the location of wave-breaking zones, while wave forces induce setup and currents. Despite this interaction, waves and circulation processes are modeled separately using different approaches. Circulation processes are represented by the shallow water equations, which conserve mass and momentum. This approach for wind-generated waves is impractical for large geographic scales due to the fine resolution that would be required. Therefore, wind-waves are instead represented in a spectral sense, governed by the action balance equation, which propagates action density through both geographic and spectral space. Even though wind-waves and circulation are modeled separately, it is important to account for their interactions by coupling their respective models. In this dissertation we use discontinuous-Galerkin (DG) methods to couple spectral wave and circulation models to model wave-current interactions. We first develop, implement, verify and validate a DG spectral wave model, which allows for the implementation of unstructured meshes in geographic space and the utility of adaptive, higher-order approximations in both geographic and spectral space. We then couple the DG spectral wave model to an existing DG circulation model, which is run on the same geographic mesh and allows for higher order information to be passed between the two models. We verify and validate coupled wave/circulation model as well as analyzing the error of the coupled wave/circulation model. / text

Storm surge, the pileup of seawater occurring as a result of high surface stresses and strong currents generated by extreme storm events such as hurricanes, is known to cause greater loss of life than these storms' associated winds. For example, inland flooding from the storm surge along the Gulf Coast during Hurricane Katrina killed hundreds of people. Previous storms produced even larger death tolls. Simultaneously, dune, barrier island, and channel erosion taking place during a hurricane leads to the removal of major flow controls, which significantly affects inland inundation. Also, excessive sea bed scouring around pilings can compromise the structural integrity of bridges, levees, piers, and buildings. Modeling these processes requires tightly coupling a bed morphology equation to the shallow water equations (SWE). Discontinuous Galerkin finite element methods (DGFEMs) are a natural choice for modeling this coupled system, given the need to solve these problems on large, complicated, unstructured computational meshes, as well as the desire to implement hp-adaptivity for capturing the dynamic features of the solution. Comprehensive modeling of these processes in the coastal zone presents several challenges and open questions. Most existing hydrodynamic models use a fixed-bed approach; the bottom is not allowed to evolve in response to the fluid motion. With respect to movable-bed models, there is no single, generally accepted mathematical model in use. Numerical challenges include coupling models of processes that exhibit disparate time scales during fair weather, but possibly similar time scales during intense storms. The main goals of this dissertation include implementing a robust, efficient, tightly-coupled morphological model using the local discontinuous Galerkin (LDG) method within the existing Advanced Circulation (ADCIRC) modeling framework, performing systematic code and model verification (using test cases with known solutions, proven convergence rates, or well-documented physical behavior), analyzing the stability and accuracy of the implemented numerical scheme by way of a priori error estimates, and ultimately laying some of the necessary groundwork needed to simultaneously model storm surges and bed morphodynamics during extreme storm events. / text

Este trabalho aborda a solução numérica das equações lineares de águas rasas. O método dos elementos finitos e utilizado para a discretização espacial das equações que modelam o problema, e para a discretização temporal, o esquema semi-implícito de Crank-Nicolson é empregado. Além de alguns conceitos comuns quando se trabalha com escoamentos geofísicos, são descritas também a formulação das equações de águas rasas, sua linearização e uma solução analítica para um caso onde o parâmetro de Coriolis é nulo. A escolha adequada de pares de elementos finitos é a principal dificuldade quando se trabalha com esse método para a resolução da equação de águas rasas. Assim, é discutido o uso de quatro pares de elementos finitos e técnicas de estabilização para contornar o surgimento de modos espúrios na solução discreta. Os resultados numéricos são realizados com auxílio do software FreeFem++, onde se pode notar a capacidade dos pares de elementos de reproduzirem o escoamento, através da solução discreta, além das propriedades de conservação de massa e energia de cada discretização. / This work is about the numerical solution of the linear shallow water equations. The finite element method is used for spatial discretization of the equations that model the problem and for the time discretization the semi-implicit Crank-Nicolson scheme is used. Besides the concepts related to geophysical flows, the formulation of the shallow water equations, their linearization and an analytical solution for a case where the Coriolis parameter is zero are also described. The appropriate choice of a pair of finite elements is the main difficulty when working with this method for solving the shallow water equations. The use of four pairs of finite elements and stabilization techniques to circumvent the appearance of spurious modes in the discrete solution are discussed. The numerical results are performed using the software FreeFem++, where one can notice the ability of the elements to represent the discrete solution and mass and energy conservation of each discretization.