Courant Condition Research Articles

A kinetic flux vector splitting scheme for shallow water flows

A kinetic flux vector splitting (KFVS) scheme for shallow water flows based on the collisionless Boltzmann equation is formulated and applied. The scheme is explicit and first order in space and time with stability governed by the Courant condition. The consistency of the KFVS scheme with the shallow water equations is proven using the equivalent differential equations approach. The accuracy and efficiency of the KFVS scheme in modeling complex flow features are compared to those of the Boltzmann Bhatnagar–Gross–Krook (BGK) scheme as well as a Riemann-based scheme. In particular, all schemes are applied to (i) strong shock waves, (ii) extreme expansion waves, (iii) a combination of strong shock waves and extreme expansion waves, and (iv) a one-dimensional dam break problem. Additionally, the KFVS, BGK and Riemann schemes are applied to a one-dimensional dam break problem for which laboratory data is available. These test cases reveal that all three schemes provide solutions of comparable accuracy, but the KFVS model is 1.5–2 times faster to execute than the BGK scheme and 2–3 times faster than the Riemann-based scheme. The absence of the collision term from the Boltzmann equation not only makes the mathematical formulation of KFVS easy but also helps elucidate this approach to the novice. The accuracy, efficiency, and simplicity of the KFVS scheme indicate its potential in modeling an array of water resources problems. Due to the scalar nature of the Boltzmann equation, the extension of the KFVS scheme to 2-D surface water flows is straightforward.

Fast transport simulation with an adaptive grid refinement.

One of the main difficulties in transport modeling and calibration is the extraordinarily long computing times necessary for simulation runs. Improved execution time is a prerequisite for calibration in transport modeling. In this paper we investigate the problem of code acceleration using an adaptive grid refinement, neglecting subdomains, and devising a method by which the Courant condition can be ignored while maintaining accurate solutions. Grid refinement is based on dividing selected cells into regular subcells and including the balance equations of subcells in the equation system. The connection of coarse and refined cells satisfies the mass balance with an interpolation scheme that is implicitly included in the equation system. The refined subdomain can move with the average transport velocity of the subdomain. Very small time steps are required on a fine or a refined grid, because of the combined effect of the Courant and Peclet conditions. Therefore, we have developed a special upwind technique in small grid cells with high velocities (velocity suppression). We have neglected grid subdomains with very small concentration gradients (zero suppression). The resulting software, MODCALIF, is a three-dimensional, modularly constructed FORTRAN code. For convenience, the package names used by the well-known MODFLOW and MT3D computer programs are adopted, and the same input file structure and format is used, but the program presented here is separate and independent. Also, MODCALIF includes algorithms for variable density modeling and model calibration. The method is tested by comparison with an analytical solution, and illustrated by means of a two-dimensional theoretical example and three-dimensional simulations of the variable-density Cape Cod and SALTPOOL experiments. Crossing from fine to coarse grid produces numerical dispersion when the whole subdomain of interest is refined; however, we show that accurate solutions can be obtained using a fraction of the execution time required by uniformly fine-grid solutions.

Vlasov model using kinetic phase point trajectories.

A method of solution of the collisionless Vlasov equation by following fixed collisionless phase point ("particle") trajectories (characteristics) in phase space is presented. It solves the coupled Vlasov Maxwell system self-consistently and employs the Leapfrog-Trapezoidal scheme to solve for the characteristics explicitly. It then uses the bilinear finite element interpolation scheme in phase space and maps vital instantaneous phase point information (distribution function) to a fixed background phase space mesh while retaining it at the phase point. The scheme is an enhanced second order one in time and fourth order in space. The code is then used to model a thermal plasma as well as two stream instability using mobile electrons and fixed background ions: the scheme being a momentum conserving one by construction allows energy conservation without assignment of particle shape functions; Langmuir waves are obtained with very good agreement with the Bohm-Gross dispersion relation; the two stream results do not show any numerically induced oscillations attributed to the initial well-ordered velocity distributions. Retention of the characteristics also minimized diffusion. Extensive numerical stability analysis deriving Courant condition for the scheme as well as behavior of computational modes are done in Appendix A, as well as estimating the impact of numerical diffusion in Appendix B. Two to five dimensional versions in phase space exist.

Reduced equations for electromagnetic turbulence in tokamaks

A set of reduced equations is derived for use in simulations of electromagnetic turbulence in tokamaks, in which the electron transit time scale is removed. The time step is not subject to a Courant condition depending on the electron thermal speed, but instead depends on the Alfvén speed. This allows significantly longer time steps in explicit numerical methods, for relevant values of the ratio of electron pressure to magnetic pressure.

Transient flow in pipe networks

In order to apply the fixed-grid method of characteristics to transient flow in multi-pipe systems, some pipe reaches may not satisfy the Courant condition. Various numerical techniques can be applied for these reaches. In the present study, the finite difference or interpolation techniques are coupled with the method of characteristics. The location of the disproportionate pipe, which does not satisfy the Courant condition, is changed along the pipe in different schemes. The accuracy of the developed schemes has been checked by comparing their results to available exact analytical solutions. The comparison indicated that the scheme in which the finite difference method is employed for an intermediate location of the disproportionate pipe reach is the best. A general computer program based on this method has been developed. This program is capable of analyzing pipe networks including pumps, valves, surge tanks, air chambers, etc.

Enhanced grid scheme for one-dimensional integro-differential Volterra equation in the time domain

An enhanced numerical method for solving a one-dimensional (1-D) integro-differential Volterra equation for time- and spatially varying dispersive dielectric media is proposed. The scheme is based on spatial-time discretization (marching scheme) with an enhanced Courant stability condition. The algorithm performs refinement of the results by a prediction-correction scheme and imposes restrictions on neither the signal shape and duration nor on the temporal behavior of the medium and is therefore generic. Results show good correlation with existing theoretical solutions and considerable improvement in precision over previous techniques. The influence of electromagnetic transients in the lossless dielectric layer on reflection and transmission coefficients is investigated.

The phase velocity error and stability condition of the three-dimensional nonstandard FDTD method

The nonstandard finite-difference time-domain (NS-FDTD) method, using a rectangular parallelepipeds structured grid, has been proposed to overcome the dispersion and anisotropic errors of the FDTD method. However, the numerical dispersion and the stability condition have not been examined. Furthermore, the method has been defined only in the isotropic grids. This paper investigates the numerical dispersion and the stability condition of the three-dimensional NS-FDTD method for isotropic and nonisotropic grids. The method is compared with the FDTD method. As a result, this method demonstrates highly accurate characteristics and high Courant stability condition.

On the stability of the FDTD algorithm for elastic media at a material interface

In this paper, the stability behavior of the first-order finite-difference time-domain algorithm for elastodynamics at the interface between two different materials is investigated. A necessary condition for stability is established, which, dependent on the material properties of the two media, might be more restrictive than the well-known Courant condition. It is shown that this more restrictive stability condition can be avoided if the material properties are averaged on the boundary.

Dynamic Time Step for One-Dimensional Overland Flow Kinematic Wave Solution

Kinematic wave theory is widely used in modeling a variety of hydrologic processes. Results of applying the kinematic wave overland flow solution using different time steps showed that the conventionally used stability criterion known as the Courant condition fails to give a time step estimate that ensures stable and accurate numerical solutions. Accordingly, a new accuracy-based dynamic time step estimate for the one dimensional overland flow kinematic wave solution is developed. The newly developed dynamic time step estimates are functions of the mesh size, watershed slope, roughness, excess rainfall, and time of concentration. The new criteria were developed by comparing the consistent formulation of the Galerkin-Crank Nicholson numerical solution of the kinematic wave equation to the characteristic method-based analytical solution, using different time steps and meshes. For each simulation, characterized by boundary and initial conditions and mesh size, an optimal time step that integrates the problem within 5% error was determined. The series of mesh sizes and corresponding optimal time steps were used to develop the dynamic time step. The time step criteria were tested on a variety of problems, including a steady state and time varying rainfall scenarios, and proved to be adequate for accurate and stable results within an efficient computational time. The criteria can be easily integrated in flow routing models to select the optimal time step with minimal user input.

A fluid–kinetic hybrid electron model for electromagnetic simulations

A fluid–kinetic hybrid electron model for electromagnetic simulations of finite-β plasmas is developed based on an expansion of the electron response using the electron–ion mass ratio as a small parameter. (Here β is the ratio of plasma pressure to magnetic pressure.) The model accurately recovers low frequency plasma dielectric responses and faithfully preserves nonlinear kinetic effects (e.g., phase space trapping). Maximum numerical efficiency is achieved by overcoming the electron Courant condition and suppressing high frequency modes. This method is most useful for nonlinear kinetic (particle-in-cell or Vlasov) simulations of electromagnetic microturbulence and Alfvénic instabilities in magnetized plasmas.

The split-weight particle simulation scheme for plasmas

An efficient numerical method for treating electrons in magnetized plasmas has been developed. The scheme, which is based on the perturbative (δf) gyrokinetic particle simulation, splits the particle electron responses into adiabatic and nonadiabatic parts. The former is incorporated into the gyrokinetic Poisson’s equation, while the latter is calculated dynamically with the aid of the charge conservation equation. The new scheme affords us the possibility of suppressing unwanted high-frequency oscillations and, in the meantime, relaxing the Courant condition for the thermal particles moving in the parallel direction. It is most useful for studying low-frequency phenomena in plasmas. As an example, one-dimensional drift wave simulation has been carried out using the scheme and the results are presented in this paper. This methodology can easily be generalized to problems in three-dimensional toroidal geometry, as well as those in unmagnetized plasmas.

FARGO: A fast eulerian transport algorithm for differentially rotating disks

We present an efficient and simple modification of the standard transport algorithm used in explicit eulerian fixed polar grid codes, aimed at getting rid of the average azimuthal velocity when applying the Courant condition. This results in a much larger timestep than the usual procedure, and it is particularly well-suited to the description of a Keplerian disk where one is traditionally limited by the very demanding Courant condition on the fast orbital motion at the inner boundary. In this modified algorithm, the timestep is limited by the perturbed velocity and by the shear arising from the differential rotation. FARGO stands for “Fast Advection in Rotating Gaseous Objects”. The speed-up resulting from the use of the FARGO algorithm is problem dependent. In the example presented here, which shows the evolution of a Jupiter sized protoplanet embedded in a minimum mass protoplanetary nebula, the FARGO algorithm is about an order of magnitude faster than a traditional transport scheme, with a much smaller numerical diffusivity.

Multi-dimensional hydrodynamic code SOVA for interfacial flows: Application to the thermal layer effect

Eulerian, three-dimensional, numerical code, which conserves mass, momentum and energy simultaneously both in the Lagrangian and remap steps, is developed. The use of special form of linear viscosity provides a weaker time step restriction as compared with the Courant condition. The code is designed to investigate the multi-material problems, including dusty flows. The performance of the code is illustrated by the modeling of shock wave interaction with a dusty thermal layer.

A comparison of models for characterizing the distribution of radionuclides with depth in soils

SummaryIt is common practice to fit mathematical models to radionuclide activity–depth profiles in soils in order to quantify rates of vertical migration through the soil profile. We have fitted six such models to 21 different activity–depth profiles of radiocaesium (137Cs) derived from Chernobyl and determined relations between the models and the values of their parameters. The advection and dispersion parameters obtained using four solutions to the advection–dispersion equation (each based on different initial and boundary conditions or different simplifications) are in good agreement. We further develop a relation between parameter values obtained using the advection–dispersion models and those determined by a simpler exponential function of the form Aexp(–Bt) where t is the time and A and B are parameters to be estimated. One of the advection–dispersion models proved to be significantly better than the others in terms of goodness‐of‐fit, versatility and ease of use. A simple model, using calculations based on measured characteristics of the activity–depth profile, was shown to accord well with parameters derived from more complex models based on statistical curve fitting. We have also evaluated the ‘residence time’ or ‘compartmental’ model approach to characterizing radionuclide activity–depth profiles. We relate such models to a numerical solution of a simple advection equation, and we show that apparent dispersion in compartmental models is an artefact of numerical dispersion, which can be quantified by the Courant condition. For activity profiles that have a significant advection component, using solutions to the advection–dispersion equation, we have observed a strong positive correlation between advection and dispersion in the profile.

An SPH Projection Method

A new formulation is introduced for enforcing incompressibility in Smoothed Particle Hydrodynamics (SPH). The method uses a fractional step with the velocity field integrated forward in time without enforcing incompressibility. The resulting intermediate velocity field is then projected onto a divergence-free space by solving a pressure Poisson equation derived from an approximate pressure projection. Unlike earlier approaches used to simulate incompressible flows with SPH, the pressure is not a thermodynamic variable and the Courant condition is based only on fluid velocities and not on the speed of sound. Although larger time-steps can be used, the solution of the resulting elliptic pressure Poisson equation increases the total work per time-step. Efficiency comparisons show that the projection method has a significant potential to reduce the overall computational expense compared to weakly compressible SPH, particularly as the Reynolds number, Re, is increased. Simulations using this SPH projection technique show good agreement with finite-difference solutions for a vortex spin-down and Rayleigh–Taylor instability. The results, however, indicate that the use of an approximate projection to enforce incompressibility leads to error accumulation in the density field.

THREE-DIMENSIONAL RESIN TRANSFER MOLDING: ISOTHERMAL PROCESS MODELING AND IMPLICIT TRACKING OF MOVING FRONTS FOR THICK, GEOMETRICALLY COMPLEX COMPOSITES MANUFACTURING APPLICATIONS - PART 2

In the manufacture of complex thick composite structures, isothermal process modeling simulation tools analyzing the three-dimensional flow of resin impregnating the thick fiber preform are instrumental in initial process optimizations. We present a new implicit formulation based on accurately taking into account the transient nature of the time-dependent conservation of resin mass and employing the pure finite element method to implicitly solve for the pressure field and to track the flow front progression of the resin inside the mold cavity. Our emphasis is on thick composites. Our previous efforts involved two-dimensional thin sections [1] and the implicit type approach. The present method is an extension of that effort to three-dimensional thick composites. The numerical developments involving the present implicit methodology and tracking of the resin progression flow fronts provide an improved and physically accurate representation of the physical problem. They do not involve the time step restrictions based on the Courant condition as in the traditional explicit finite element-control volume (FE-CV) associated formulation discussed in Part 1 in this issue. Instead, we treat the transient process as a series of quasi-steady-state processes. For illustrative purposes, the present developments are first validated with simple geometries, then extended to practical applications for isothermal conditions, and finally, contrasted to part 1 in this issue of this study.

A finite-difference time-domain method without the Courant stability conditions

In this paper, a finite-difference time-domain method that is free of the constraint of the Courant stability condition is presented for solving electromagnetic problems. In it, the alternating direction implicit (ADI) technique is applied in formulating the finite-difference time-domain (FDTD) algorithm. Although the resulting formulations are computationally more complicated than the conventional FDTD, the proposed FDTD is very appealing since the time step used in the simulation is no longer restricted by stability but by accuracy. As a result, computation speed can be improved. It is found that the number of iterations with the proposed FDTD can be at least three times less than that with the conventional FDTD with the same numerical accuracy.

The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures

AbstractWe present the spectral element method to simulate elastic-wave propagation in realistic geological structures involving complieated free-surface topography and material interfaces for two- and three-dimensional geometries. The spectral element method introduced here is a high-order variational method for the spatial approximation of elastic-wave equations. The mass matrix is diagonal by construction in this method, which drastically reduces the computational cost and allows an efficient parallel implementation. Absorbing boundary conditions are introduced in variational form to simulate unbounded physical domains. The time discretization is based on an energy-momentum conserving scheme that can be put into a classical explicit-implicit predictor/multi-corrector format. Long-term energy conservation and stability properties are illustrated as well as the efficiency of the absorbing conditions. The associated Courant condition behaves as ΔtC < O (nel−1/ndN−2), with nel the number of elements, nd the spatial dimension, and N the polynomial order. In practice, a spatial sampling of approximately 5 points per wavelength is found to be very accurate when working with a polynomial degree of N = 8. The accuracy of the method is shown by comparing the spectral element solution to analytical solutions of the classical two-dimensional (2D) problems of Lamb and Garvin. The flexibility of the method is then illustrated by studying more realistic 2D models involving realistic geometries and complex free-boundary conditions. Very accurate modeling of Rayleigh-wave propagation, surface diffraction, and Rayleigh-to-body-wave mode conversion associated with the free-surface curvature are obtained at low computational cost. The method is shown to provide an efficient tool to study the diffraction of elastic waves by three-dimensional (3D) surface topographies and the associated local effects on strong ground motion. Complex amplification patterns, both in space and time, are shown to occur even for a gentle hill topography. Extension to a heterogeneous hill structure is considered. The efficient implementation on parallel distributed memory architectures will allow to perform real-time visualization and interactive physical investigations of 3D amplification phenomena for seismic risk assessment.

Implicit SPH Drag and Dusty Gas Dynamics

In a multiphase flow such as gas and dust, the time step for numerical integration is controlled by the Courant condition, by the drag between the fluids, and by gravity. If the dust particles are sufficiently small the time step is determined by the drag. Since the primary effect of very large drag is to keep the fluids moving together, and the details of how this is achieved are not important when the drag is very large, it is desirable to treat the drag terms implicitly. For standard grid-based methods this is straightforward, but for particle methods like SPH it appears difficult because the drag appears in the form of pair interactions between particles, and any given particle can interact with ∼40 neighbours. The idea exploited in this paper is to treat each pair interaction separately. The velocities of the two particles involved are updated implicitly, according to their pair interaction, and the initial velocities of these particles are then replaced by the new velocities. This process is repeated for each pair interaction. Two sweeps over the particles gives satisfactory convergence. In this paper the method is tested first by keeping just the drag terms, then including pressure and drag to test wave propagation, and finally including gravity to consider the fall of a layer of gas in an isothermal atmosphere. The results are in good agreement with theory. The basic idea of working with implicit pair interactions can be extended to other SPH problems.

Adaptive Local Refinement with Octree Load Balancing for the Parallel Solution of Three-Dimensional Conservation Laws

Conservation laws are solved by a local Galerkin finite element procedure with adaptive space-time mesh refinement and explicit time integration. The Courant stability condition is used to select smaller time steps on smaller elements of the mesh, thereby greatly increasing efficiency relative to methods having a single global time step. Processor load imbalances, introduced at adaptive enrichment steps, are corrected by using traversals of an octree representing a spatial decomposition of the domain. To accommodate the variable time steps, octree partitioning is extended to use weights derived from element size. Partition boundary smoothing reduces the communications volume of partitioning procedures for a modest cost. Computational results comparing parallel octree and inertial partitioning procedures are presented for the three-dimensional Euler equations of compressible flow solved on an IBM SP2 computer.