The thermal radiative transfer (TRT) equations form an integrodifferential system that describes the propagation and collisional interactions of photons. Computing accurate and efficient numerical solutions TRT are challenging for several reasons, the first of which is that TRT is defined on a highdimensional phase space that includes the independent variables of time, space, and velocity. In order to reduce the dimensionality of the phase space, classical approaches such as the P$_N$ (spherical harmonics) or the S$_N$ (discrete ordinates) ansatz are often used in the literature. In this work, we introduce a novel approach: the hybrid discrete (H$^T_N$) approximation to the radiative thermal transfer equations. This approach acquires desirable properties of both P$_N$ and S$_N$, and indeed reduces to each of these approximations in various limits: H$^1_N$ $\equiv$ P$_N$ and H$^T_0$ $\equiv$ S$_T$. We prove that H$^T_N$ results in a system of hyperbolic partial differential equations for all $T\ge 1$ and $N\ge 0$. Another challenge in solving the TRT system is the inherent stiffness due to the large timescale separation between propagation and collisions, especially in the diffusive (i.e., highly collisional) regime. This stiffness challenge can be partially overcome via implicit time integration, although fully implicit methods may become computationally expensivemore »
A selfadaptive theta scheme using discontinuity aware quadrature for solving conservation laws
Abstract We present a discontinuity aware quadrature (DAQ) rule and use it to develop implicit selfadaptive theta (SATh) schemes for the approximation of scalar hyperbolic conservation laws. Our SATh schemes require the solution of a system of two equations, one controlling the cell averages of the solution at the time levels, and the other controlling the spacetime averages of the solution. These quantities are used within the DAQ rule to approximate the time integral of the hyperbolic flux function accurately, even when the solution may be discontinuous somewhere over the time interval. The result is a finite volume scheme using the theta time stepping method, with theta defined implicitly (or selfadaptively). Two schemes are developed, selfadaptive theta upstream weighted (SAThup) for a monotone flux function using simple upstream stabilization, and selfadaptive theta Lax–Friedrichs (SAThLF) using the Lax–Friedrichs numerical flux. We prove that DAQ is accurate to second order when there is a discontinuity in the solution and third order when it is smooth. We prove that SAThup is unconditionally stable, provided that theta is set to be at least 1/2 (which means that SATh can be only first order accurate in general). We also prove that SAThup satisfies the maximum more »
 Award ID(s):
 1912735
 Publication Date:
 NSFPAR ID:
 10356572
 Journal Name:
 IMA Journal of Numerical Analysis
 ISSN:
 02724979
 Sponsoring Org:
 National Science Foundation
More Like this


Simulation of flow and transport in petroleum reservoirs involves solving coupled systems of advectiondiffusionreaction equations with nonlinear flux functions, diffusion coefficients, and reactions/wells. It is important to develop numerical schemes that can approximate all three processes at once, and to high order, so that the physics can be well resolved. In this paper, we propose an approach based on high order, finite volume, implicit, Weighted Essentially NonOscillatory (iWENO) schemes. The resulting schemes are locally mass conservative and, being implicit, suited to systems of advectiondiffusionreaction equations. Moreover, our approach gives unconditionally Lstable schemes for smooth solutions to the linear advectiondiffusionreaction equation in the sense of a von Neumann stability analysis. To illustrate our approach, we develop a third order iWENO scheme for the saturation equation of twophase flow in porous media in two space dimensions. The keys to high order accuracy are to use WENO reconstruction in space (which handles shocks and steep fronts) combined with a twostage RadauIIA RungeKutta time integrator. The saturation is approximated by its averages over the mesh elements at the current time level and at two future time levels; therefore, the scheme uses two unknowns per grid block per variable, independent of the spatial dimension. Thismore »

The Langevin Dynamics (LD) method (also known in the literature as Brownian Dynamics) is routinely used to simulate aerosol particle trajectories for transport rate constant calculations as well as to understand aerosol particle transport in internal and external fluid flows. This tutorial intends to explain the methodological details of setting up a LD simulation of a population of aerosol particles and to deduce rate constants from an ensemble of classical trajectories. We discuss the applicability and limitations of the translational Langevin equation to model the combined stochastic and deterministic motion of particles in fields of force or fluid flow. The drag force and stochastic “diffusion” force terms that appear in the Langevin equation are discussed elaborately, along with a summary of common forces relevant to aerosol systems (electrostatic, gravity, van der Waals, …); a commonly used first order and a fourth order RungeKutta time stepping schemes for linear stochastic ordinary differential equations are presented. A MATLAB® implementation of a LD code for simulating particle settling under gravity using the first order scheme is included for illustration. Scaling analysis of aerosol transport processes and the selection of timestep and domain size for trajectory simulations are demonstrated through two specific aerosol processes:more »

Seismic waves in earth media usually undergo attenuation, causing energy losses and phase distortions. In the regime of highfrequency asymptotics, a complexvalued eikonal is an essential ingredient for describing wave propagation in attenuating media, where the real and imaginary parts of the eikonal function capture dispersion effects and amplitude attenuation of seismic waves, respectively. Conventionally, such a complexvalued eikonal is mainly computed either by tracing rays exactly in complex space or by tracing rays approximately in real space so that the resulting eikonal is distributed irregularly in real space. However, seismic data processing methods, such as prestack depth migration and tomography, usually require uniformly distributed complexvalued eikonals. Therefore, we have developed a unified framework to Eulerianize several popular approximate realspace raytracing methods for complexvalued eikonals so that the real and imaginary parts of the eikonal function satisfy the classic realspace eikonal equation and a novel realspace advection equation, respectively, and we dub the resulting method the Eulerian partialdifferentialequation method. We further develop highly efficient highorder methods to solve these two equations by using the factorization idea and the LaxFriedrichs weighted essentially nonoscillatory schemes. Numerical examples demonstrate that our method yields highly accurate complexvalued eikonals, analogous to those from raytracing methods.more »

The YBJ equation (Young & Ben Jelloul, J. Marine Res. , vol. 55, 1997, pp. 735–766) provides a phaseaveraged description of the propagation of nearinertial waves (NIWs) through a geostrophic flow. YBJ is obtained via an asymptotic expansion based on the limit $\mathit{Bu}\rightarrow 0$ , where $\mathit{Bu}$ is the Burger number of the NIWs. Here we develop an improved version, the YBJ + equation. In common with an earlier improvement proposed by Thomas, Smith & Bühler ( J. Fluid Mech. , vol. 817, 2017, pp. 406–438), YBJ + has a dispersion relation that is secondorder accurate in $\mathit{Bu}$ . (YBJ is firstorder accurate.) Thus both improvements have the same formal justification. But the dispersion relation of YBJ + is a Padé approximant to the exact dispersion relation and with $\mathit{Bu}$ of order unity this is significantly more accurate than the powerseries approximation of Thomas et al. (2017). Moreover, in the limit of high horizontal wavenumber $k\rightarrow \infty$ , the wave frequency of YBJ + asymptotes to twice the inertial frequency $2f$ . This enables solution of YBJ + with explicit timestepping schemes using a time step determined by stable integration of oscillations with frequency $2f$ . Other phaseaveraged equations have dispersion relations with frequency increasing as $k^{2}$more »