The Maxwell's equations are written as normal Hamilton equations using functional variation method. We discretize Maxwell's equations using sympletic propagation technique combined with fourth-order finite difference approximations to construct symplectic finite difference time domain (SFDTD) scheme...

In this paper we study numerically the KdV-top equation and compare it with the Boussinesq equations over uneven bottom. We use here a finite-difference scheme that conserves a discrete energy for the fully discrete scheme. We also compare this approach with the discontinuous Galerkin method. For th...

The method of equivariant moving frames on multi-space is used to construct symmetry preserving finite difference schemes of partial differential equations invariant under finite-dimensional symmetry groups. Invariant numerical schemes for a heat equation with a logarithmic source and the spherical Burgers equation are obtained. Numerical tests show how invariant schemes can be more accurate than standard discretizations on uniform rectangular meshes.

This paper defines discretizations of the divergence and flux operators that produce symmetric, positive-definite, and accurate approximations to steady-state diffusion problems. Because discontinuous material properties and highly distorted grids are allowed, the flux operator, rather than the gradient, is used as a fundamental operator to be discretized. Resulting finite-difference scheme is similar to those obtained from the mixed finite-element method.

The impetus of this paper is the comparative applications of two numerical schemes for supersonic flows using computational algorithms tailored for a supercomputer. The mathematical model is the conservation form of Navier-Stokes equations with the effect of turbulence being modeled algebraically. The first scheme is an implicit, unfactored, upwind-biased, line-Gauss-Seidel relaxation scheme based on finite-volume discretization. The second scheme is the explicit-implicit MacCormack scheme based on finite-difference discretization. The best overall efficiencies are obtained using the upwind relaxation scheme. The integrity of the solutions obtained for the example cases is shown by comparisons with experimental and other computational results.

We present a new non-conforming space-time mesh refinement method for the symmetric first order hyperbolic system. This method is based on the one hand on the use of a conservative higher order discontinuous Galerkin approximation for space discretization and a finite difference scheme in time, on the other hand on appropriate discrete transmission conditions between the grids. We use a discrete energy technique to drive the construction of the matching procedure between the grids and guarantee the stability of the method.

We consider the discretization problem for U(1)-invariant nonlinear wave equations in any dimension. We show that the classical finite-difference scheme used by Strauss and Vazquez (in J. Comput. Phys. 28, 271?278 (1978)) conserves the positive-definite discrete analog of the energy if the grid ratio satisfies Formula Not Shown , where dt and dx are the mesh sizes of the time and space variables and n is the spatial dimension. We also show that, if the grid ratio is Formula Not Shown , then there is a discrete analog of charge, and this discrete analog is conserved. We prove the existence and uniqueness of solutions to the discrete Cauchy problem. We use energy conservation to obtain a priori bounds for finite energy solutions, thus showing that the Strauss-Vazquez finite-difference scheme...

A new benchmarking concept is presented for verifying the PEBBED 3D multigroup finite difference/nodal diffusion code with application to pebble bed modular reactors (PBMRs). The key idea is to perform convergence acceleration, also called extrapolation to zero discretization, of a basic finite difference numerical algorithm to give extremely high accuracy. The method is first demonstrated on a 1D cylindrical shell and then on an r,8 wedge where the order of the second order finite difference scheme is confirmed to four places.

A mixed finite-difference scheme is presented for the free-vibration analysis of simply supported closed noncircular cylindrical shells. The problem is formulated in terms of eight first-order differential equations in the circumferential coordinate which possess a symmetric coefficient matrix and are free of the derivatives of the elastic and geometric characteristics of the shell. In the finite-difference discretization, two interlacing grids are used for the different fundamental unknowns in such a way as to avoid averaging in the difference-quotient expressions used for the first derivative. The resulting finite-difference equations are symmetric. The inverse-power method is used for obtaining the eigenvalues and eigenvectors.

Numerical schemes for Einstein's vacuum equation are developed. Einstein's equation in harmonic gauge is second-order symmetric hyperbolic. It is discretized in four-dimensional spacetime by finite differences, finite elements and interior penalty discontinuous Galerkin methods, the latter being related to Regge calculus. The schemes are split into space and time and new time-stepping schemes for wave equations are derived. The methods are evaluated for linear and nonlinear test problems of the Apples-with-Apples collection.

We study the differential-difference regularization for a two-dimensional inverse heat conduction problem, i.e. the heat equation is semi-discretized by a differential-difference equation, where the time derivative and a spatial second-order derivative have been replaced by the finite differences, while the other spatial second-order derivative is preserved. We analyze the properties of the discretized approximation using Fourier transform techniques. Some error estimates, which give the information about how to choose the step lengths in the discretization, show that the semi-discretized form has a 'regularized effect'. We also proved the unconditional convergence of a discretization scheme involving spatial marching.

We present a new non conforming space-time mesh refinement method for symmetric first order hyperbolic system. This method is based on the one hand on the use of a conservative higher order discontinuous Galerkin approximation for space discretization and a finite difference scheme in time, on the o...

Most of thermal hydraulic processes in nuclear engineering can be described by general convection-diffusion equations that are often can be simulated numerically with finite-difference method (FDM). An effective scheme for finite-difference discretization of such equations is presented in this report. The derivation of this scheme is based on analytical solutions of a simplified one-dimensional equation written for every control volume of the finite-difference mesh. These analytical solutions are constructed using linearized representations of both diffusion coefficient and source term. As a result, the Efficient Finite-Differencing (EFD) scheme makes it possible to significantly improve the accuracy of numerical method even using mesh systems with fewer grid nodes that, in turn, allows to speed-up numerical simulation. EFD has been carefully verified on the series of sample problems for which either analytical or very precise numerical solutions can be found. EFD has been compared with other popular FDM schemes including novel, accurate (as well as sophisticated) methods. Among the methods compared were well-known central difference scheme, upwind scheme, exponential differencing and hybrid schemes of Spalding. Also, newly developed finite-difference schemes, such as the the quadratic upstream (QUICK) scheme of Leonard, the locally analytic differencing (LOAD) scheme of Wong and Raithby, the flux-spline scheme proposed by Varejago and Patankar as well as the latest LENS discretization of Sakai have been compared. Detailed results of this comparison are given in this report. These tests have shown a high efficiency of the EFD scheme. For most of sample problems considered EFD has demonstrated the numerical error that appeared to be in orders of magnitude lower than that of other discretization methods. Or, in other words, EFD has predicted numerical solution with the same given numerical error but using much fewer grid nodes. In this report, the detailed description of EFD is given. It includes basic assumptions, the detailed derivation, the verification procedure, as well verification and comparisons. Conclusion summarizes results and highlights the problems to be solved. (author)

We investigate the connections between several recent methods for the discretization of anisotropic heterogeneous diffusion operators on general grids. We prove that the Mimetic Finite Difference scheme, the Hybrid Finite Volume scheme and the Mixed Finite Volume scheme are in fact identical up to some slight generalizations. As a consequence, some of the mathematical results obtained for each of the method (such as convergence properties or error estimates) may be extended to the unified common framework. We then focus on the relationships between this unified method and nonconforming Finite Element schemes or Mixed Finite Element schemes, obtaining as a by-product an explicit lifting operator close to the ones used in some theoretical studies of the Mimetic Finite Difference scheme. We also show that for isotropic operators, on particular meshes such as triangular meshes with acute angles, the unified method boils down to the well-known efficient two-point flux Finite Volume scheme.

The relationship is analyzed between layer-resolving transformations and mesh-generating functions for numerical solution of singularly perturbed boundary-value problems. The analysis is carried out for one-dimensional quasilinear problems without turning points, which are discretized by first-order finite-difference schemes. It is proved that if a general layer-resolving function is used to generate the discretization mesh, then the numerical solution converges uniformly in the perturbation parameter.

In the past several years, with the development of the coarse-mesh schemes for the spatial discretization of the diffusion equations, the computing time for the static neutron flux calculations has been reduced considerably over the conventional finite difference scheme. For the temporal discretization in the time-dependent case, a fully implicit method (FIM) was found to be more efficient than either the explicit or the alternating-direction explicit method. In this paper, a new approach using the collocation method is investigated for the temporal discretization in an effort to further improve the computing time over that of the FIM.

The NEWT (NEW Transport algorithm) code is a multi-group discrete ordinates neutral-particle transport code with flexible meshing capabilities. This code employs the Extended Step Characteristic spatial discretization approach using arbitrary polygonal mesh cells. Until recently, the coarse mesh finite difference acceleration scheme in NEWT for fission source iteration has been available only for rectangular domain boundaries because of the limitation to rectangular coarse meshes. Therefore no acceleration scheme has been available for triangular or hexagonal problem boundaries. A conventional and a new partial-current based coarse mesh finite difference acceleration schemes with unstructured coarse meshes have been implemented within NEWT to support any form of domain boundaries. The computational results show that the new acceleration schemes works well, with performance often improved over the earlier two-level rectangular approach.

In this paper, a high-order technique to accurately predict flow-generated noise is introduced. The technique consists of solving the viscous incompressible flow equations and inviscid acoustic equations using a incompressible/compressible splitting technique. The incompressible flow equations are solved using the in-house flow solver EllipSys2D/3D which is a second-order finite volume code. The acoustic solution is found by solving the acoustic equations using high-order finite difference schemes. The incompressible flow equations and the acoustic equations are solved at the same time levels where the pressure and the velocities obtained from the incompressible equations form the input to the acoustic equations. To achieve low dissipation and dispersion errors, either Dispersion-Relation-Preserving (DRP) schemes or optimized compact finite difference schemes are used for spatial discretizations of the acoustic equations. The classical fourth-order Runge-Kutta time scheme is applied to the acoustic equations for time discretization.

The utility of the finite-element Galerkin technique in advection-diffusion flow problems is examined by comparison with several finite-difference schemes in one dimension. The calculations show that for relatively coarse grids, finite-element solutions are either comparable to or significantly better than those obtained from the finite-difference schemes considered. For advection-dominated flows, the superiority of the finite-element technique is attributed to spatial coupling of time-derivative terms inherent in the Galerkin discretization. This procedure, absent from conventional finite-difference schemes, leads to very accurate phase properties for the approximate solution even when coarse grids are used. A two-dimensional analogue of the advection-diffusion problem further illustrates the advantages and accuracy of the finite-element method in conjunction with the use of isoparametric elements.

We propose a finite-difference algorithm for solving the time-dependent Ginzburg-Landau (TDGL) equation coupled to the appropriate Maxwell equation. The time derivatives are discretized using a second order semi-implicit scheme which, for intermediate values of the Ginzburg-Landau parameter $\\kappa$, allows time-steps two orders of magnitude larger than commonly used in explicit schemes. We demonstrate the use of the method by solving a fully three-dimensional problem of a current-carrying wire with longitudinal and transverse magnetic fields.

In the present paper, we gauge the performance of finite-difference schemes with Runge-Kutta time integration for wave propagation problems by rigorously defining appropriate cost and error metrics in a simple setting represented by the linear advection equation. Optimal values of the grid spacing and of the time step are obtained as a result of a cost minimization (for given error level) procedure. The theory suggests superior performance of high-order schemes when highly accurate solutions are sought for, and in several space dimensions even more. The analysis of the global discretization error shows the occurrence of two (approximately independent) sources of error, associated with the space and time discretizations. The improvement of the performance of finite-difference schemes can th...

Surface topography has been considered a difficult task for seismic wave numerical modelling by the finite-difference method (FDM) because the most popular staggered finite-difference scheme requires a rectilinear grid. Even though there are numerous collocated grid schemes in other computational fields that could be used to solve the first-order hyperbolic equations, the lack of a stable free-surface boundary condition implementation for curvilinear grids also obstructs the adoption of curvilinear grids in seismic wave FDM modelling. In this study, we use generalized curvilinear grids that can fit the surface topography to discretize the computational domain and describe the implementation of a collocated grid finite-difference scheme, a higher order MacCormack scheme, to solve the first-order hyperbolic velocity-stress equations on the curvilinear grid. To achieve a sufficiently accurate and stable free-surface boundary condition implementation on the curvilinear grids, we propose the traction image method that antisymmetrically images the traction components instead of the stress components to the ghost points above the free surface. Since the velocity derivatives at the free surface are provided by the free-surface condition, we use a compact scheme to compute the velocity derivatives near the free surface and avoid the use of velocity values on the ghost points. Numerical tests verify that using the curvilinear grid, the collocated finite-difference scheme and the traction image technique can simulate seismic wave propagation in the presence of surface topography with sufficient accuracy.

In this work, some ordering schemes for mesh points are presented which enable algorithms such as the Gauss-Seidel or SOR iteration to be performed efficiently for the nine-point operator finite difference method on computers consisting of a two-dimensional grid of processors. Convergence results are presented for the discretization of u /SUB xx/ + u /SUB yy/ on a uniform mesh over a square, showing that the spectral radius of the iteration for these orderings is no worse than that for the standard row by row ordering of mesh points. Further applications of these mesh point orderings to network problems, more general finite difference operators, and picture processing problems are noted.

The numerical solution for a class of fractional sub-diffusion equations is studied. For the time discretization, we use a generalized Crank?Nicolson method combined with the second central finite difference (FD) for the spatial discretization which will then define a fully discrete implicit FD scheme. An error of order O(h 2 max (1, log?k ???1)?+?k 2?+?? ) has been shown where h and k denote the maximum space and time steps, respectively. A non-uniform time step is employed to compensate for the singular behaviour of the exact solution at?t?=?0. Our theoretical results are numerically validated in a series of test problems.

We consider solving time-independent (steady-state) flow problems in 2D or 3D governed by hyperbolic or {open_quotes}almost hyperbolic{close_quotes} systems of partial differential equations (PDE). Examples of such PDE are the Euler and the Navier-Stokes equations. The PDE is discretized using a finite difference or finite volume scheme with arbitrary order of accuracy. If the matrix B describes the discretized differential operator and u denotes the approximate solution, the discrete problem is given by a large system of equations.

A mixed finite-difference scheme is presented for the stress and free vibration analysis of simply supported nonhomogeneous and layered orthotropic thick plates. The analytical formulation is based on the linear, three-dimensional theory of orthotropic elasticity and a Fourier approach is used to reduce the governing equations to six first-order ordinary differential equations in the thickness coordinate. The governing equations possess a symmetric coefficient matrix and are free of derivatives of the elastic characteristics of the plate. In the finite difference discretization two interlacing grids are used for the different fundamental unknowns in such a way as to reduce both the local discretization error and the bandwidth of the resulting finite-difference field equations. Numerical studies are presented for the effects of reducing the interior and boundary discretization errors and of mesh refinement on the accuracy and convergence of solutions. It is shown that the proposed scheme, in addition to a number of other advantages, leads to highly accurate results, even when a small number of finite difference intervals is used.

In this paper we present and compare two unconditionally energy stable finite-difference schemes for the phase field crystal equation. The first is a one-step scheme based on a convex splitting of a discrete energy by Wise et al. [S.M. Wise, C. Wang, J.S. Lowengrub, An energy stable and convergent finite-difference scheme for the phase field crystal equation, SIAM J. Numer. Anal., in press]. In this scheme, which is first order in time and second order in space, the discrete energy is non-increasing for any time step. The second scheme we consider is a new, fully second-order two-step algorithm. In the new scheme, the discrete energy is bounded by its initial value for any time step. In both methods, the equations at the implicit time level are nonlinear but represent the gradients of strictly convex functions and are thus uniquely solvable, regardless of time step-size. We solve the nonlinear equations using an efficient nonlinear multigrid method. Numerical simulations are presented and confirm the stability, efficiency and accuracy of the schemes.

The energy of solutions of the wave equation with a suitable boundary dissipation decays exponentially to zero as time goes to infinity. We consider the finite-difference space semi-discretization scheme and we analyze whether the decay rate is independent of the mesh size. We focus on the one-dimensional case. First we show that the decay rate of the energy of the classical semi-discrete system in which the 1?d Laplacian is replaced by a three-point finite difference scheme is not uniform with respect to the net-spacing size h. Actually, the decay rate tends to zero as h goes to zero. Then we prove that adding a suitable vanishing numerical viscosity term leads to a uniform (with respect to the mesh size) exponential decay of the energy of solutions. This numerical viscosity term damps ou...

We present a new finite-difference numerical method to solve the incompressible Navier-Stokes equations using a collocated discretization in space on a logically Cartesian grid. The method shares some common aspects with, and it was inspired by, the Box scheme. It uses centered second-order-accurate finite-difference approximations for the spatial derivatives combined with semi-implicit time integration. The proposed method is constructed to ensure discrete conservation of mass and momentum by discretizing the primitive velocity-pressure form of the equations. The continuity equation is enforced exactly (to machine accuracy) at the collocated locations, whereas the momentum equations are evaluated in a staggered manner. This formulation preempts the appearance of spurious pressure modes in...

A one dimensional fractional diffusion model with the Riemann-Liouville fractional derivative is studied. First, a second order discretization for this derivative is presented and then an unconditionally stable weighted average finite difference method is derived. The stability of this scheme is established by von Neumann analysis. Some numerical results are shown, which demonstrate the efficiency and convergence of the method. Additionally, some physical properties of this fractional diffusion system are simulated, which further confirm the effectiveness of our method.

In this paper we present and analyze an implicit fully discrete local discontinuous Galerkin (LDG) finite element method for solving the time-fractional Schrödinger equation, where the fractional derivative is described in the Caputo sense. The scheme is based on a finite difference method in time and local discontinuous Galerkin methods in space. A stability and error analysis is performed on the numerical methods. Numerical results confirm the expected convergence rates and illustrate the effectiveness of the method.

In this article, we study the dynamics of Mackey-Glass system applied a nonstandard finite-difference scheme. For the discrete system we show that a sequence of Hopf bifurcations occur at the positive fixed point as the delay increasing, analyze the stability of the fixed point and calculate the direction of the Hopf bifurcations. At last, by giving some numerical experiments, we illustrate the relation between the time of producing blood cells and symptom.

Within the framework of the fractal mobile-immobile medium model describing non-Fickian effects occurring in admixture seepage due to particle adhesion to the solid matrix, an expression for the admixture flux is derived. Flow discretization intended for finite-difference calculations is proposed and used as a basis for a conservation-law scheme for solving the model equations with account for admixture sources. Several one-dimensional test problems of admixture propagation in an imposed seepage flow are solved using the approach developed.

In this paper, a fully discrete local discontinuous Galerkin (LDG) finite element method is considered for solving the time-fractional KdV-Burgers-Kuramoto (KBK) equation. The scheme is based on a finite difference method in time and local discontinuous Galerkin methods in space. We prove that our scheme is unconditional stable and $L^2$ error estimate for the linear case with the convergence rate $O(h^{k+1}+(\\Delta t)^2+(\\Delta t)^\\frac{\\alpha}{2}h^{k+1/2})$. Numerical examples are presented to show the efficiency and accuracy of our scheme.

The aim of this paper is to present the class of high order compact schemes in the context of numerical simulation of stratified flow. The numerical schemes presented here are based on the approach outlined in Lele [1]. The numerical model presented in this contribution is based on the solution of the Boussinesq approximation by a finite-difference scheme. The numerical scheme itself follows the principle of semi-discretization, with high order compact discretization in space, while the time integration is carried out by suitable Runge-Kutta time-stepping scheme. In the case presented here the steady flow was considered and thus the artificial compressibility method was used to resolve the pressure from the modified continuity equation. The test case used to demonstrate the capabilities of the selected model consists of the flow of stably stratified fluid over low, smooth hill.

The major thrust of this proposal was to continue our investigations of so-called non-standard finite-difference schemes as formulated by other authors. These schemes do not follow the standard rules used to model continuous differential equations by discrete difference equations. The two major aspects of this procedure consist of generalizing the definition of the discrete derivative and using a nonlocal model (on the computational grid or lattice) for nonlinear terms that may occur in the differential equations. Our aim was to investigate the construction of nonstandard finite-difference schemes for several classes of ordinary and partial differential equations. These equations are simple enough to be tractable, yet, have enough complexity to be both mathematically and scientifically interesting. It should be noted that all of these equations differential equations model some physical phenomena under an appropriate set of experimental conditions. The major goal of the project was to better understand the process of constructing finite-difference models for differential equations. In particular, it demonstrates the value of using nonstandard finite-difference procedures. A secondary goal was to construct and study a variety of analytical techniques that can be used to investigate the mathematical properties of the obtained difference equations. These mathematical procedures are of interest in their own right and should be a valuable contribution to the mathematics research literature in difference equations. All of the results obtained from the research done under this project have been published in the relevant research/technical journals or submitted for publication. Our expectation is that these results will lead to improved finite difference schemes for the numerical integration of both ordinary and partial differential equations. Section G of the Appendix gives a concise summary of the major results obtained under funding by the grant.

For a Spacelab flight, a model experiment of the earth's atmospheric circulation has been proposed. This experiment is known as the Atmospheric General Circulation Experiment (AGCE). In the experiment concentric spheres will rotate as a solid body, while a dielectric fluid is confined in a portion of the gap between the spheres. A zero gravity environment will be required in the context of the simulation of the gravitational body force on the atmosphere. The present study is concerned with the development of pseudospectral/finite difference (PS/FD) model and its subsequent application to physical cases relevant to the AGCE. The model is based on a hybrid scheme involving a pseudospectral latitudinal formulation, and finite difference radial and time discretization. The advantages of the use of the hybrid PS/FD method compared to a pure second-order accurate finite difference (FD) method are discussed, taking into account the higher accuracy and efficiency of the PS/FD method.

A higher order accurate numerical procedure has been developed for solving incompressible Navier-Stokes equations for 2D or 3D fluid flow problems. It is based on low-storage Runge-Kutta schemes for temporal discretization and fourth and sixth order compact finite-difference schemes for spatial discretization. The particular difficulty of satisfying the divergence-free velocity field required in incompressible fluid flow is resolved by solving a Poisson equation for pressure. It is demonstrated that for consistent global accuracy, it is necessary to employ the same order of accuracy in the discretization of the Poisson equation. Special care is also required to achieve the formal temporal accuracy of the Runge-Kutta schemes. The accuracy of the present procedure is demonstrated by application to several pertinent benchmark problems.

The three-dimensional, time-dependent convection-diffusion equation (CDE) is considered. An exponential transformation is used to collectively transform the CDE. The idea of global discretization is used, where attention is paid to the whole transformed CDE, but not to the individual spatial and temporal derivatives in the equation. Four finite difference schemes for both CDE and transformed CDE are established. The modified partial differential equations of these schemes are obtained, which indicate that the trunction errors of the schemes can be of second and fourth order, depending on the prescription of the time step length. Some characteristic physical parameters, i.e., local Reynolds number, local Strouhal number, and viscous diffusive length, are introduced into the schemes and the viscous diffusive length is found to be a significant parameter in relating temporal discretisation with spatial discretisation. A series of benchmark analytical solutions of Navier-Stokes and Burgers equations, as well s the numerical solutions using the well-known discretisation schemes, are used to investigate the properties of the derived schemes. The high-order schemes achieve higher resolutions over the conventional schemes without decreasing much the sparsity of the matrix structures. Grid refinement studies reveal that the inverse exponential transformation of the finite difference schemes tends to destroy some resolution of the schemes, especially for large local Reynolds number. 18 refs., 9 figs., 4 tabs.

This paper examines a class of explicit finite-difference advection schemes derived along the method of lines. An important application field is large-scale atmospheric transport. The paper therefore focuses on the demand of positivity. For the spatial discretization, attention is confined to conservative schemes using five points per direction. The fourth-order central scheme and the family of {kappa}-schemes, comprising the second-order central, the second-order upwind, and the third-order upwind biased, are studied. Positivity is enforced through flux limiting. It is concluded that the limited third-order upwind discretization is the best candidate from the four examined. For the time integration attention is confined to a number of explicit Runge-Kutta methods of orders two to four. With regard to the demand of positivity, these integration methods turn out to behave almost equally and no best method could be identified. 16 refs., 4 figs., 4 tabs.

Von Neumann established that discretized algebraic equations must be consistent with the differential equations, and must be stable in order to obtain convergent numerical solutions for the given differential equations. The "stability" is required to satisfactorily approximate a differential derivative by its discretized form, such as a finite-difference scheme, in order to compute in computers. His criterion is the necessary and sufficient condition only for steady or equilibrium problems. It is also a necessary condition, but not a sufficient condition for unsteady transient problems; additional care is required to ensure the accuracy of unsteady solutions.

A three-dimensional numerical model for large-eddy simulation (LES) of oceanic turbulent processes is described. The numerical formulation comprises a spectral discretization in the horizontal directions and a high-order compact finite-difference discretization in the vertical direction. Time-stepping is accomplished via a second-order accurate fractional-step scheme. LES subgrid-scale (SGS) closure is given by a traditional Smagorinsky eddy-viscosity parametrization for which the model coefficient is derived following similarity theory in the near-surface region. Alternatively, LES closure is given by the dynamic Smagorinsky parametrization for which the model coefficient is computed dynamically as a function of the flow. Validation studies are presented demonstrating the temporal and spa...

Hyperbolic systems of conversation laws are considered which are discretized in space by spectral collocation methods and advanced in time by finite difference schemes. At any time-level a domain deposition method based on an iteration by subdomain procedure was introduced yielding at each step a sequence of independent subproblems (one for each subdomain) that can be solved simultaneously. The method is set for a general nonlinear problem in several space variables. The convergence analysis, however, is carried out only for a linear one-dimensional system with continuous solutions. A precise form of the error reduction factor at each iteration is derived. Although the method is applied here to the case of spectral collocation approximation only, the idea is fairly general and can be used in a different context as well. For instance, its application to space discretization by finite differences is straight forward.

A computational tool called "Directional Diffusion Regulator (DDR)" is proposed to bring forth real multidimensional physics into the upwind discretization in some numerical schemes of hyperbolic conservation laws. The direction based regulator when used with dimension splitting solvers, is set to moderate the excess multidimensional diffusion and hence cause genuine multidimensional upwinding like effect. The basic idea of this regulator driven method is to retain a full upwind scheme across local discontinuities, with the upwind bias decreasing smoothly to a minimum in the farthest direction. The discontinuous solutions are quantified as gradients and the regulator parameter across a typical finite volume interface or a finite difference interpolation point is formulated based on fractio...

Finite difference schemes for transient convection-diffusion problems on grids with local refinement in time and space are constructed and studied. The construction utilizes a modified upwind approximation and linear interpolation at the slave nodes. The proposed schemes are implicit of backward Euler type and unconditionally stable. Error analysis is presented in the maximum norm, and convergence estimates are derived for smooth solutions. Optimal approximation results for ratios between the spatial and time discretization parameters away from the CFL condition are shown. Finally, numerical examples illustrating the theory are given.

The method of difference potentials was originally proposed by Ryaben?kii and can be interpreted as a generalized discrete version of the method of Calderon?s operators in the theory of partial differential equations. It has a number of important advantages; it easily handles curvilinear boundaries, variable coefficients, and non-standard boundary conditions while keeping the complexity at the level of a finite-difference scheme on a regular structured grid. The method of difference potentials assembles the overall solution of the original boundary value problem by repeatedly solving an auxiliary problem. This auxiliary problem allows a considerable degree of flexibility in its formulation and can be chosen so that it is very efficient to solve. Compact finite difference schemes enable hig...

In this paper a finite-difference-based lattice Boltzmann method for curvilinear coordinates is proposed in order to improve the computational efficiency and numerical stability of a recent method [R. Mei and W. Shyy, J. Comput. Phys. 143, 426 (1998)] in which the collision term of the Boltzmann Bhatnagar-Gross-Krook equation for discrete velocities is treated implicitly. In the present method, the implicitness of the numerical scheme is removed by introducing a distribution function different from that being used currently. As a result, an explicit finite-difference lattice Boltzmann method for curvilinear coordinates is obtained. The scheme is applied to a two-dimensional Poiseuille flow, an unsteady Couette flow, a lid-driven cavity flow, and a steady flow around a circular cylinder. The numerical results are in good agreement with the results of previous studies. Extensions to other lattice Boltzmann models based on nonuniform meshes are also discussed.

Previous applications of QUICK for the discretization of convective transport terms in finite-volume calculation procedures have failed to employ a rigorous and systematic approach for consistently deriving this finite difference scheme. Instead, earlier formulations have been established numerically, by trial and error. The new formulation for QUICK presented here is obtained by requiring that it satisfy four rules that guarantee physically realistic numerical solutions having overall balance. Careful testing performed for the wall-driven square enclosure flow configuration shows that the consistently derived version of QUICK is more stable and converges faster than any of the formulations previously employed. This testing includes the relative evaluation of boundary conditions approximated by second- and third-order finite-difference schemes as well as calculations performed at higher Reynolds numbers than previously reported. 14 refs., 10 figs., 1 tab.

We describe a general optimization procedure for both maximizing the resolution characteristics of existing finite differencing schemes as well as designing finite difference schemes that will meet the error tolerance requirements of numerical solutions. The procedure is based on an optimization process. This is a generalization of the compact scheme introduced by Lele in which the resolution is improved for single, one-dimensional spatial derivative, whereas in the present approach the complete scheme, after spatial and temporal discretizations, is optimized on a range of parameters of the scheme and the governing equations. The approach is to linearize and Fourier analyze the discretized equations to check the resolving power of the scheme for various wave number ranges in the solution and optimize the resolution to satisfy the requirements of the problem. This represents a constrained nonlinear optimization problem which can be solved to obtain the nodal weights of discretization. An objective function is defined in the parametric space of wave numbers, Courant number, Mach number and other quantities of interest. Typical criterion for defining the objective function include the maximization of the resolution of high wave numbers for acoustic and electromagnetic wave propagations and turbulence calculations. The procedure is being tested on off-design conditions of non-uniform mesh, non-periodic boundary conditions, and non-constant wave speeds for scalar and system of equations. This includes the solution of wave equations and Euler equations using a conventional scheme with and without optimization and the design of an optimum scheme for the specified error tolerance.

In the present paper, we gauge the performance of finite-difference schemes with Runge Kutta time integration for wave propagation problems by rigorously defining appropriate cost and error metrics in a simple setting represented by the linear advection equation. Optimal values of the grid spacing and of the time step are obtained as a result of a cost minimization (for given error level) procedure. The theory suggests superior performance of high-order schemes when highly accurate solutions are sought for, and in several space dimensions even more. The analysis of the global discretization error shows the occurrence of two (approximately independent) sources of error, associated with the space and time discretizations. The improvement of the performance of finite-difference schemes can then be achieved by trying to separately minimize the two contributions. General guidelines for the design of problem-tailored, optimized schemes are provided, suggesting that significant reductions of the computational cost are in principle possible. The application of the analysis to wave propagation problems in a two-dimensional environment demonstrates that the analysis carried out for the scalar case directly applies to the propagation of monochromatic sound waves. For problems of sound propagation involving disparate length-scales the analysis still provides useful insight for the optimal exploitation of computational resources; however, the actual advantage provided by optimized schemes is not as evident as in the single-scale, scalar case.

A method is proposed which allows to efficiently treat elliptic problems on unbounded domains in two and three spatial dimensions in which one is only interested in obtaining accurate solutions at the domain boundary. The method is an extension of the optimal grid approach for elliptic problems, based on optimal rational approximation of the associated Neumann-to-Dirichlet map in Fourier space. It is shown that, using certain types of boundary discretization, one can go from second-order accurate schemes to essentially spectrally accurate schemes in two-dimensional problems, and to fourth-order accurate schemes in three-dimensional problems without any increase in the computational complexity. The main idea of the method is to modify the impedance function being approximated to compensate for the numerical dispersion introduced by a small finite-difference stencil discretizing the differential operator on the boundary. We illustrate how the method can be efficiently applied to nonlinear problems arising in modeling of cell communication.

We study the numerical approximation to the solution of the steady convection-diffusion equation. The diffusion term is discretized by using the hybrid mimetic method (HMM), which is the unified formulation for the hybrid finite-volume (FV) method, the mixed FV method and the mimetic finite-difference method recently proposed in Droniou et al. (2010, Math. Models Methods Appl. Sci., 20, 265-295). In such a setting we discuss several techniques to discretize the convection term that are mainly adapted from the literature on FV or FV schemes. For this family of schemes we provide a full proof of convergence under very general regularity conditions of the solution field and derive an error estimate when the scalar solution is in H2(). Finally, we compare the performance of these schemes on a ...

A detailed finite-difference approach is presented for the simulation of transient radial flow in multi-layer systems. The proposed discretization scheme simulates drawdown within the well more accurately than commonly applied schemes. The solution is compared to existing (semi) analytical models for the simulation of slug tests and pumping tests with constant discharge in single- and multi-layer systems. For all cases, it is concluded that the finite-difference model approximates drawdown to acceptable accuracy. The main advantage of finite-difference approaches is the ability to account for the varying saturated thickness in unconfined top layers. Additionally, it is straightforward to include radial variation of hydraulic parameters, which is useful to simulate the effect of a finite-thickness well skin. Aquifer tests with variable pumping rate and/or multiple wells may be simulated by superposition. The finite-difference solution is implemented in MAxSym, a MATLAB tool which is designed specifically to simulate axi-symmetric flow. Alternatively, the presented equations can be solved using a standard finite-difference model. A procedure is outlined to apply the same approach with MODFLOW. The required modifications to the input parameters are much larger for MODFLOW than for MAxSym, but the results are virtually identical. The presented finite-difference solution may be used, for example, as a forward model in parameter estimation algorithms. Since it is applicable to multi-layer systems, its use is not limited to the simulation of traditional pumping and slug tests, but also includes advanced aquifer tests, such as multiple pumping tests or multi-level slug tests.

A systematic analysis shows how results from the finite difference code SEAWAT are sensitive to choice of grid dimension, time step, and numerical scheme for unstable flow problems. Guidelines to assist in selecting appropriate combinations of these factors are suggested. While the SEAWAT code has been tested for a wide range of problems, the sensitivity of results to spatial and temporal discretization levels and numerical schemes has not been studied in detail for unstable flow problems. Here, the Elder-Voss-Souza benchmark problem has been used to systematically explore the sensitivity of SEAWAT output to spatio-temporal resolution and numerical solver choice. A grid size of 0.38 and 0.60% of the total domain length and depth respectively is found to be fine enough to deliver results with acceptable accuracy for most of the numerical schemes when Courant number (Cr) is 0.1. All numerical solvers produced similar results for extremely fine meshes; however, some schemes converged faster than others. For instance, the 3rd-order total variation-diminishing method (TVD3) scheme converged at a much coarser mesh than the standard finite difference methods (SFDM) upstream weighting (UW) scheme. The sensitivity of the results to Cr number depends on the numerical scheme as expected.

We analyse 13 3-D numerical time-domain explicit schemes for modelling seismic wave propagation and earthquake motion for their behaviour with a varying P-wave to S-wave speed ratio (VP/VS). The second-order schemes include three finite-difference, three finite-element and one discontinuous-Galerkin schemes. The fourth-order schemes include three finite-difference and two spectral-element schemes. All schemes are second-order in time. We assume a uniform cubic grid/mesh and present all schemes in a unified form. We assume plane S-wave propagation in an unbounded homogeneous isotropic elastic medium. We define relative local errors of the schemes in amplitude and the vector difference in one time step and normalize them for a unit time. We also define the equivalent spatial sampling ratio as a ratio at which the maximum relative error is equal to the reference maximum error. We present results of the extensive numerical analysis. We theoretically (i) show how a numerical scheme sees the P and S waves if the VP/VS ratio increases, (ii) show the structure of the errors in amplitude and the vector difference and (iii) compare the schemes in terms of the truncation errors of the discrete approximations to the second mixed and non-mixed spatial derivatives. We find that four of the tested schemes have errors in amplitude almost independent on the VP/VS ratio. The homogeneity of the approximations to the second mixed and non-mixed spatial derivatives in terms of the coefficients of the leading terms of their truncation errors as well as the absolute values of the coefficients are key factors for the behaviour of the schemes with increasing VP/VS ratio. The dependence of the errors in the vector difference on the VP/VS ratio should be accounted for by a proper (sufficiently dense) spatial sampling.

In the present paper we consider the Navier-Stokes equations for the two-dimensional viscous incompressible fluid flows and apply to these equations our earlier designed general algorithmic approach to generation of finite-difference schemes. In doing so, we complete first the Navier-Stokes equations to involution by computing their Janet basis and discretize this basis by its conversion into the integral conservation law form. Then we again complete the obtained difference system to involution with eliminating the partial derivatives and extracting the minimal Gröbner basis from the Janet basis. The elements in the obtained difference Gröbner basis that do not contain partial derivatives of the dependent variables compose a conservative difference scheme. By exploiting arbitrariness in the numerical integration approximation we derive two finite-difference schemes that are similar to the classical scheme by Harlow and Welch. Each of the two schemes is characterized by a 5×5 stencil on an orthogonal and uniform grid. We also demonstrate how an inconsistent difference scheme with a 3×3 stencil is generated by an inappropriate numerical approximation of the underlying integrals.

In this paper, a new 27-point finite difference method is presented for solving the 3D Helmholtz equation with perfectly matched layer (PML), which is a second order scheme and pointwise consistent with the equation. An error analysis is made between the numerical wavenumber and the exact wavenumber, and a refined choice strategy based on minimizing the numerical dispersion is proposed for choosing weight parameters. A full-coarsening multigrid-based preconditioned Bi-CGSTAB method is developed for solving the linear system stemming from the Helmholtz equation with PML by the finite difference scheme. The shifted-Laplacian is extended to precondition the 3D Helmholtz equation, and a spectral analysis is given. The discrete preconditioned system is solved by the Bi-CGSTAB method, with a multigrid method used to invert the preconditioner approximately. Full-coarsening multigrid is employed, and a new matrix-based prolongation operator is constructed accordingly. Numerical results are presented to demonstrate the efficiency of both the new 27-point finite difference scheme with refined parameters, and the preconditioned Bi-CGSTAB method with the 3D full-coarsening multigrid.

To simulate seismic wave propagation in the spherical Earth, the Earth's curvature has to be taken into account. This can be done by solving the seismic wave equation in spherical coordinates by numerical methods. In this paper, we use an optimized, collocated-grid finite-difference scheme to solve the anisotropic velocity-stress equation in spherical coordinates. To increase the efficiency of the finite-difference algorithm, we use a non-uniform grid to discretize the computational domain. The grid varies continuously with smaller spacing in low velocity layers and thin layer regions and with larger spacing otherwise. We use stress-image setting to implement the free surface boundary condition on the stress components. To implement the free surface boundary condition on the velocity components, we use a compact scheme near the surface. If strong velocity gradient exists near the surface, a lower-order scheme is used to calculate velocity difference to stabilize the calculation. The computational domain is surrounded by complex-frequency shifted perfectly matched layers implemented through auxiliary differential equations (ADE CFS-PML) in a local Cartesian coordinate. We compare the simulation results with the results from the normal mode method in the isotropic and anisotropic models and verify the accuracy of the finite-difference method.

The Artificial Compressibility Method (ACM) for the incompressible Navier-Stokes equations is (link-wise) reformulated (referred to as LW-ACM) by a finite set of discrete directions (links) on a regular Cartesian mesh, in analogy with the Lattice Boltzmann Method (LBM). The main advantage is the possibility of exploiting well established technologies originally developed for LBM and classical computational fluid dynamics, with special emphasis on finite differences (at least in the present paper), at the cost of minor changes. For instance, wall boundaries not aligned with the background Cartesian mesh can be taken into account by tracing the intersections of each link with the wall (analogously to LBM technology). LW-ACM requires no high-order moments beyond hydrodynamics (often referred to as ghost moments) and no kinetic expansion. Like finite difference schemes, only standard Taylor expansion is needed for analyzing consistency. Preliminary efforts towards optimal implementations have shown that LW-ACM is...

A coupled fluid-structure method to study the coupled nonlinear problem of flapping aero-elastic structures was developed. The fluid-dynamic solver is a finite-difference solution to the Navier-Stokes equations solved on a grid fitted to a moving thin body. Direct solution of the NavierStokes equations is carried out using a previously developed high-order 2D immersed boundary method on a moving curvilinear grid. Fluid-dynamic forcing on the body surface is calculated and used as input to a finite difference structural-dynamic solver. The structural solver is geometrically nonlinear and able to take on arbitrary configurations. Both solvers are explicit and use a Runge-Kutta 4th order time discretization scheme. Case studies that the coupled fluid/structure method was applied include flag flapping and harmonically pitching flexible membranes with different densities and rigidities.

In this paper, we introduce a general family of variational, multisymplectic numerical methods for solving Maxwell's equations, using discrete differential forms in spacetime. In doing so, we demonstrate several new results, which apply both to some well-established numerical methods and to new methods introduced here. First, we show that Yee's finite-difference time-domain (FDTD) scheme, along with a number of related methods, are multisymplectic and derive from a discrete Lagrangian variational principle. Second, we generalize the Yee scheme to unstructured meshes, not just in space, but in 4-dimensional spacetime. This relaxes the need to take uniform time steps, or even to have a preferred time coordinate at all. Finally, as an example of the type of methods that can be developed within this general framework, we introduce a new asynchronous variational integrator (AVI) for solving Maxwell's equations. These results are illustrated with some prototype simulations that show excellent energy and conservatio...

An explicit finite-difference scheme is presented for solving the two-dimensional Biot equations of poroelasticity across the full range of frequencies. The key difficulty is to discretize the Johnson-Koplik-Dashen (JKD) model which describes the viscous dissipations in the pores. Indeed, the time-domain version of Biot-JKD model involves order 1/2 shifted fractional derivatives which amounts to a time convolution product. To avoid storing the past values of the solution, a diffusive representation of fractional derivatives is used: the convolution kernel is replaced by a finite number of memory variables that satisfy local-in-time ordinary differential equations. The coefficients of the diffusive representation follow from an optimization procedure of the dispersion relation. Then, various methods of scientific computing are applied: the propagative part of the equations is discretized using a fourth-order ADER scheme, whereas the diffusive part is solved exactly. An immersed interface method is implemented ...

In this article, a new finite element method, discontinuous finite difference streamline diffusion method (DFDSD), is constructed and studied for first-order linear hyperbolic problems. This method combines the benefit of the discontinuous Galerkin method and the streamline diffusion finite element method. Two fully discrete DFDSD schemes (Euler DFDSD and Crank-Nicolson (CN) DFDSD) are constructed by making use of the difference discrete method for time variables and the discontinuous streamline diffusion method for space variables. The stability and optimal L^2 norm error estimates are established for the constructed schemes. This method makes contributions to the discontinuous methods. Finally, a numerical example is provided to show the benefit of high efficiency and simple implementati...

The numerical solution for a class of sub-diffusion equations involving a parameter in the range - 1 a < 0 is studied. For the time discretization, we use an implicit finite-difference Crank-Nicolson method and show that the error is of order k2+a, where k denotes the maximum time step. A nonuniform time step is employed to compensate for the singular behaviour of the exact solution at t = 0. We also consider a fully discrete scheme obtained by applying linear finite elements in space to the proposed time-stepping scheme. We prove that the additional error is of order h2 max(1, log k-1), where h is the parameter for the space mesh. Numerical experiments on some sample problems demonstrate our theoretical result.

Difficulties for the conventional computational fluid dynamics and the standard lattice Boltzmann method (LBM) to study the gas oscillating patterns in a resonator have been discussed. In light of the recent progresses in the LBM world, we are now able to deal with the compressibility and non-linear shock wave effects in the resonator. A lattice Boltzmann model for viscid compressible flows is introduced firstly. Then, the Boltzmann equation with the Bhatnagar-Gross-Krook approximation is solved by the finite-difference method with a third-order implicit-explicit (IMEX) Runge-Kutta scheme for time discretization, and a fifth-order weighted essentially non-oscillatory (WENO) scheme for space discretization. Numerical results obtained in this study agree quantitatively with both experimental...

The sensitivity of a numeric steam flooding model with respect to mobility weighting is examined in depth. Three numeric discretization procedures are used in this investigation: a new numeric scheme, a 5-point finite difference method, and a procedure which, under certain assumptions, is equivalent to that introduced by McCracken and Yanosik. Three mobility weighting schemes also are investigated: (1) upstream mobility weighting; (2) harmonic total mobility weighting; and (3) upstream weighting of fractional flow terms. The approach introduced uses the kinematic viscosity in the total mobility and the fractional flow terms. The steam displacement model formed from the combination of this mobility weighting approach and the McCracken and Yanosik discretization procedure is shown to produce realistic simulations of an inverted 7-spot pattern under a continuous steam drive. 20 references.

In recent years, there has been a growing interest in the development and use of mathematical models for the simulation of fluid flow, heat transfer and combustion processes in engineering equipment. The equations representing the multi-dimensional transport of mass, momenta and species are numerically solved by finite-difference or finite-element techniques. However despite the multiude of differencing schemes and solution algorithms, and the advancement of computing power, the calculation of multi-dimensional flows, especially three-dimensional flows, remains a mammoth task. The following discussion is concerned with the author's recent work on the construction of accurate discretization schemes for the partial derivatives, and the efficient solution of the set of nonlinear algebraic equations resulting after discretization. The present work has been jointly supported by the Ramjet Engine Division of the Wright Patterson Air Force Base, Ohio, and the NASA Lewis Research Center.

In this paper we consider the stability and convergence of finite-difference discretizations of a reaction-diffusion equation on a one-dimensional domain that is growing in time. We consider discretizations of conservative and nonconservative formulations of the governing equation and highlight the different stability characteristics of each. Although nonconservative formulations are the most popular to date, we find that discretizations of the conservative formulation inherit greater stability properties. Furthermore, we present a novel adaptive time integration scheme based on the well-known method and describe how the parameter should be chosen to ensure unconditional stability, independently of the rate of domain growth. This work is a preliminary step towards an analysis of numerical ...

Numerical techniques are developed to solve the Navier-Stokes equations for unsteady incompressible flow. The extension of the finite-difference Galerkin (FDG) method of Stephens et al. (1984) to the continuous-time case in two or three space dimensions is explained, and the numerical implementation of the method is discussed with particular attention to the staggered-MAC-grid primitive-variable discretization, the application of discrete mass balance to avoid problems inherent in FDG schemes, the direct interpretation of the FDG expansion variables as a discrete streamfunction, and a mass-balance approach to two-dimensional problems with throughflow or obstacles. Numerical results are presented graphically for the evolution of asymptotic steady flow in a driven cavity at Reynolds number 400, 1000, or 3200; good agreement with published experimental data is demonstrated, with accurate predictions of secondary-vortex formation from wall bubble recirculations at Reynolds number 1000.

We build Gaussian wave packets for the linear Schr\\"odinger equation and its finite difference space semi-discretization and illustrate the lack of uniform dispersive properties of the numerical solutions as established in Ignat, Zuazua, Numerical dispersive schemes for the nonlinear Schr\\"odinger equation, SIAM. J. Numer. Anal., 47(2) (2009), 1366-1390. It is by now well known that bigrid algorithms provide filtering mechanisms allowing to recover the uniformity of the dispersive properties as the mesh size goes to zero. We analyze and illustrate numerically how these high frequency wave packets split and propagate under these bigrid filtering mechanisms, depending on how the fine grid/coarse grid filtering is implemented.

We present a finite-difference scheme which solves the Stokes problem in the presence of curvilinear non-conforming interfaces and provides second-order accuracy on physical field (velocity, vorticity) and especially on pressure. The gist of our method is to rely on the Helmholtz decomposition of the Stokes equation: the pressure problem is then written in an integral form devoid of the spurious sources known to be the cause of numerical boundary layer error in most implementations, leading to a discretization which guarantees a strict enforcement of mass conservation. The ghost method is furthermore used to implement the boundary values of pressure and vorticity near curved interfaces.

In this short note, we show how to use a highly accurate finite-difference scheme to compute second derivatives in the Navier-Stokes equations while ensuring targeted numerical dissipation. This approach, essentially non conservative, is shown to be close to an upwind method and is straightforward to implement with a negligible computational extra cost. The benefit offered by the resulting discrete operator is illustrated for the direct computation of sound in aeroacoustics and in the more general context of large-eddy simulation through connections with hyperviscosity and spectral vanishing viscosity.

Initial-boundary value problems for self-adjoint parabolic equations on a semiaxis and a semibounded strip are considered. For finite-difference ?-schemes, an alternative method for stating approximate transparent boundary conditions is suggested and conditions ensuring unconditional stability in the energy norm with respect to the initial data and free terms for a weight ? ? 1/2 are presented. The validity of these stability conditions in the case of discrete transparent boundary conditions is proved (by several methods), and the derivation of the latter conditions is revisited.

In this research, we present the first MPI-CUDA implementation of Finite-Difference Time-Domain (FDTD) discretization of Maxwell's equations in dispersive media that uses the MPI API to assign each CPU node its share of the computational domain and GPUs to their corresponding CPU threads. By taking advantage of the CUDA programming model, we present a unique implementation of the FDTD scheme that exploits the memory hierarchy of GPUs, including the global, texture, and shared memory. This enables us to tackle problems that are otherwise computationally prohibitive. Practical results will be presented along with a measure of speedup factors achieved when using multiple GPU processors.

A set of governing partial differential equations (PDE) derived from fundamental physical principles can describe the behavior of a fluid, but due to nonlinearity they cannot be analytically solved. Instead, they must be approximated. This work is a survey of well established finite-difference methods applied to two sets of equations-- the linear advection equation (LAE) and the linearized shallow-water equations (LSWE). Finite difference schemes replace the partial derivatives of a variable with the differences between discrete points in space and time. The resulting equations only approximate the original PDE leading to unwanted behavior, such as computational instability, damping, dispersion, and unphysical solutions. The following numerical schemes were applied to both the LAE and the LSWE: forward-in-time-and-space, Euler, backward, leapfrog, Lax-Wendroff. While there are many types of numerical schemes used for solving PDE, the schemes we employed show the necessity of weighing the characteristics of a particular scheme against the physical behavior being simulated and the resources available for computation.

Reservoir simulators have traditionally used centered finite differences with upstream weighting to discretize the spatial terms in the partial differential equations describing fluid flow in porous media. These techniques are relatively simple to implement in huge simulation codes, and they produce highly stable results which are sufficiently accurate for black-oil problems. However, as more compositional processes are being studied, it appears that the upstream weighting may degrade the accuracy of the usual discretization schemes to the point that they cannot represent the sharp fronts and rapid velocity fluctuations associated with these more complex problems. A weighting technique based on mixed finite elements which reduces numerical dispersion by lowering the phase velocities and fluid dispersal at the frontal interface is presented. The derivation of this weighting and its implementation into a standard, finite difference compositional model is described in detail, and a method for treating zero permeability grid cells with the mixed method weighting is discussed. It is also shown that the use of this new weighting technique for both immiscible and miscible processes produces sharper saturation and composition profiles and more rapid frontal advances than standard finite differences with upstream weighting. (A.V.)

It is considered the classical nonlinear shallow-water model (SWM) of an ideal fluid. It is well known that the model conserves several integral characteristics such as the mass, total energy and potential enstrophy. It is extremely desirable to conserve the same characteristics in a fully discrete SWM (discrete both in space and time), because they guarantee the stability of calculations and correct description of the energy cascades in the discrete model. However, the full discretization usually destroys some, if not all, of the conservation laws, and therefore, the construction of conservative fully discrete SWMs is a non-trivial and actual scientific problem. For the last forty years there have been suggested various semi-discrete SWMs (discrete in space, but still continuous in time), which conserve one or all of the three above-mentioned integral characteristics. In particular, the model by Ringler and Randall (2002) uses rather complicated geodesic grids on a sphere, while that by Salmon (2004) applies a sophisticated stencil (containing 25 nodes) in a doubly periodic domain on the f-plane. Nevertheless, explicit time discretization used in both works resulted in the loss of all the conservation laws except the mass conservation. In this work, new fully discrete SWMs are suggested, which exactly conserve the mass and total energy. The splitting of the SWM operator in geometric coordinates provides substantial benefits in the computational cost of the solution, as well as in the applicability to a doubly periodic domain on the plane, in a periodic channel on a rotating sphere, and on the whole sphere. Each split one-dimensional fully discrete system conserves the mass and total energy, too. In fact, a family of finite-difference schemes of different approximation order is suggested, either linear or nonlinear, depending on the choice of certain parameters. Note that on a sphere and in a doubly periodic domain, our approach allows constructing various linear conservative schemes of arbitrary approximation order in space. Results of numerical experiments are discussed.

In this paper we demonstrate that some well-known finite-difference schemes can be interpreted within the framework of the local discontinuous Galerkin (LDG) methods using the low-order piecewise solenoidal discrete spaces introduced in (SIAM J. Numer. Anal. 1990; 27(6): 1466-1485). In particular, it appears that it is possible to derive the well-known MAC scheme using a first-order Nedelec approximation on rectangular cells. It has been recently interpreted within the framework of the Raviart-Thomas approximation by Kanschat (Int. J. Numer. Meth. Fluids 2007; published online). The two approximations are algebraically equivalent to the MAC scheme, however, they have to be applied on grids that are staggered on a distance h/2 in each direction. This paper also demonstrates that both discre...

This short note reports the extension of the f-waves approximate Riemann solver (Ahmad and Lindeman, 2007; LeVeque, 2002; Bale et al., 2002) for three-dimensional meso- and micro-scale atmospheric flows. The Riemann solver employs flux-based wave decomposition for the calculation of Godunov fluxes and does not require the explicit definition of the Roe matrix to enforce conservation. The other important feature of the Riemann solver is its ability to incorporate source term due to gravity without introducing discretization errors. The resulting finite volume scheme is second-order accurate in space and time. The finite-difference schemes currently used in atmospheric flow models are neither conservative nor able to resolve regions of sharp gradients. The finite volume scheme described in t...

Finite difference methods for solving the wave equation more accurately capture the physics of waves propagating through the earth than asymptotic solution methods. Unfortunately, finite difference simulations for 3D elastic wave propagation are expensive. The authors model waves in a 3D isotropic elastic earth. The wave equation solution consists of three velocity components and six stresses. The partial derivatives are discretized using 2nd-order in time and 4th-order in space staggered finite difference operators. Staggered schemes allow one to obtain additional accuracy (via centered finite differences) without requiring additional storage. The serial code is most unique in its ability to model a number of different types of seismic sources. The parallel implementation uses the MPI library, thus allowing for portability between platforms. Spatial parallelism provides a highly efficient strategy for parallelizing finite difference simulations. In this implementation, one can decompose the global problem domain into one-, two-, and three-dimensional processor decompositions with 3D decompositions generally producing the best parallel speedup. Because I/O is handled largely outside of the time-step loop (the most expensive part of the simulation) the authors have opted for straight-forward broadcast and reduce operations to handle I/O. The majority of the communication in the code consists of passing subdomain face information to neighboring processors for use as ghost cells. When this communication is balanced against computation by allocating subdomains of reasonable size, they observe excellent scaled speedup. Allocating subdomains of size 25 x 25 x 25 on each node, they achieve efficiencies of 94% on 128 processors. Numerical examples for both a layered earth model and a homogeneous medium with a high-velocity blocky inclusion illustrate the accuracy of the parallel code.

We present a numerical method to solve the optimal transport problem with a quadratic cost when the source and target measures are periodic probability densities. This method relies on a numerical resolution of the corresponding Monge-Amp\\`ere equation. We use an existing Newton-like algorithm that we generalize to the case of a non uniform final density. The main idea consists of designing an iterative scheme where the fully nonlinear equation is approximated by a non-constant coefficient linear elliptic PDE that we discretize and solve at each iteration, in two different ways: a second order finite difference scheme and a fast Fourier transform (FFT) method. The FFT method, made possible thanks to a preconditioning step based on the coefficient-averaged equation, results in an overall O(P log P)-operations algorithm, where P is the number of discretization points. In particular, we use fourth order finite differences to approximate the action of the densities on the solution iterates, which result in more a...

Direct numerical simulation (DNS) is used to study turbulent incompressible flow past a bluff body with a compliant surface. We use a 3D time-dependent coordinate transformation to account for the motion of the bluff body surface. Spatially, the flow domain is discretized using a dealiased pseudospectral method in the axial and azimuthal directions, while the radial (wall-normal) direction is discretized using a finite difference scheme. The grid is stretched in the azimuthal direction, which is handled spectrally. This leads to a unique challenge when solving the Poisson equation in the fractional step method for the time march, which we address with both multigrid and preconditioned BiCGStab algorithms. We are presently extending this flow code with a model for the compliant bluff body surface based on the ``tensegrity fabric'' paradigm which combines compressive members (bars) and tensile members (tendons) in a stable, flexible network.

Wave propagation in a stratified fluid/porous medium is studied here using analytical and numerical methods. The semi-analytical method is based on an exact stiffness matrix method coupled with a matrix conditioning procedure, preventing the occurrence of poorly conditioned numerical systems. Special attention is paid to calculating the Fourier integrals. The numerical method is based on a high order finite-difference time-domain scheme. Mesh refinement is applied near the interfaces to discretize the slow compressional diffusive wave predicted by the Biot theory. Finally, an immersed interface method is used to discretize the boundary conditions. The numerical benchmarks are based on realistic soil parameters and on various degrees of hydraulic contact at the fluid/porous boundary. The ti...

This study demonstrates an immersed boundary (IB) method which integrates a depth-averaged two dimensional flow model is proposed to tackle a typical fluid-solid phase problem in fluid dynamics field. The finite-difference scheme with curvilinear coordinate system is employed to discretize the shallow-water flow equations. Lagrangian markers and Eulerian grid are applied to portray the geometric contour of interior boundary and discretize the flow domain, respectively. The Dirac delta function is accordingly conducted to link both Lagrangian and Eulerian coordinate systems. The numerical simulations of single pier are performed and compared to examine the effect of marker's mesh width, grid size, and the various Dirac delta functions. Experimental data from literatures are compared with nu...

Abstract We consider the modeling of (polarized) seismic wave propagation on a rectangular domain via the discretization and solution of the inhomogeneous Helmholtz equation in 3D, by exploiting a parallel multifrontal sparse direct solver equipped with Hierarchically Semi-Separable (HSS) structure to reduce the computational complexity and storage. In particular, we are concerned with solving this equation on a large domain, for a large number of different forcing terms in the context of seismic problems in general, and modeling in particular. We resort to a parsimonious mixed grid finite differences scheme for discretizing the Helmholtz operator and Perfect Matched Layer boundaries, resulting in a non-Hermitian matrix. We make use of a nested dissection based domain decomposition, and in...

This paper presents and evaluates the numerical solution of a coupled system of equations that arises in a model for the formation and evolution of three-dimensional longshore sand ridges. The model is based on the interaction between surficial or internal weakly nonlinear shallow-water waves, having weak spanwise spatial dependence, and the deformable bottom topography. The presentation of the details concerning the discretization of the model is primarily motivated by two facts: (1) The model involves equations for which little is known regarding its solutions, and (2) the predictor-corrector scheme presented here, which combines finite difference techniques and fixed-point methods, is simple, fast, and general enough to be used in the discretization of partial differential equations with local nonlinearities whose solutions are smooth.

Large-eddy simulation (LES) of turbulent boundary layer induced by external oscillatory flow without or with a unidirectional component over rippled bottom is presented. These flows represent prototype turbulent flows of wave-current interactions in the coastal zone. The numerical method is based on a time-splitting scheme for the temporal discretization and a finite-differences approximation on orthogonal grid for the spatial discretization. The immersed boundary (IMB) method is utilized to represent complex boundary shapes, i.e., the rippled bed, on the orthogonal grid. Results are presented for oscillatory flow over ripples of three steepness values and for oscillatory-unidirectional flow of two current magnitudes, including comparisons to available experimental data. In general, the ef...

Steam flooding is a tertiary oil recovery mechanism with proven economic potential. Reliable numerical simulations of candidate injection schemes can aid in the optimization of process parameters. State-of-the-art numerican models exhibit varying degrees of grid effects which affect their reliability with respect to the modeling of pattern floods. This paper presents a numerical steam flood model which does not exhibit a significant amount of grid orientation. This model utilizes a numerical discretization technique which is an admixture of finite difference and finite element methods. Fluid properties are determined as a function of primary variables in a manner that allows a straight forward implementation of Newton's method for solving the nonlinear system of discretized equations. Fractional flows are defined using mass rather than volumetric mobilities, and are used to upwind the convective flow terms.

In the field of computational fluid dynamics, many numerical algorithms have been developed to simulate inviscid, compressible flows problems. Among those most famous and relevant are based on flux vector splitting and Godunov-type schemes. Previously, this system was developed through computational studies by Mawlood [1]. However the new test cases for compressible flows, the shock tube problems namely the receding flow and shock waves were not investigated before by Mawlood [1]. Thus, the objective of this study is to develop a high-order compact (HOC) finite difference solver for onedimensional Euler equation. Before developing the solver, a detailed investigation was conducted to assess the performance of the basic third-order compact central discretization schemes. Spatial discretization of the Euler equation is based on flux-vector splitting. From this observation, discretization of the convective flux terms of the Euler equation is based on a hybrid flux-vector splitting, known as the advection upstream splitting method (AUSM) scheme which combines the accuracy of flux-difference splitting and the robustness of flux-vector splitting. The AUSM scheme is based on the third-order compact scheme to the approximate finite difference equation was completely analyzed consequently. In one-dimensional problem for the first order schemes, an explicit method is adopted by using time integration method. In addition to that, development and modification of source code for the one-dimensional flow is validated with four test cases namely, unsteady shock tube, quasi-one-dimensional supersonic-subsonic nozzle flow, receding flow and shock waves in shock tubes. From these results, it was also carried out to ensure that the definition of Riemann problem can be identified. Further analysis had also been done in comparing the characteristic of AUSM scheme against experimental results, obtained from previous works and also comparative analysis with computational results generated by van Leer, KFVS and AUSMPW schemes. Furthermore, there is a remarkable improvement with the extension of the AUSM scheme from first-order to third-order accuracy in terms of shocks, contact discontinuities and rarefaction waves.

Phase field models offer a systematic physical approach for investigating complex multiphase systems behaviors such as near-critical interfacial phenomena, phase separation under shear, and microstructure evolution during solidification. However, because interfaces are replaced by thin transition regions (diffuse interfaces), phase field simulations require resolution of very thin layers to capture the physics of the problems studied. This demands robust numerical methods that can efficiently achieve high resolution and accuracy, especially in three dimensions. We present here an accurate and efficient numerical method to solve the coupled Cahn-Hilliard/Navier-Stokes system, known as Model H, that constitutes a phase field model for density-matched binary fluids with variable mobility and viscosity. The numerical method is a time-split scheme that combines a novel semi-implicit discretization for the convective Cahn-Hilliard equation with an innovative application of high-resolution schemes employed for direct numerical simulations of turbulence. This new semi-implicit discretization is simple but effective since it removes the stability constraint due to the nonlinearity of the Cahn-Hilliard equation at the same cost as that of an explicit scheme. It is derived from a discretization used for diffusive problems that we further enhance to efficiently solve flow problems with variable mobility and viscosity. Moreover, we solve the Navier-Stokes equations with a robust time-discretization of the projection method that guarantees better stability properties than those for Crank-Nicolson-based projection methods. For channel geometries, the method uses a spectral discretization in the streamwise and spanwise directions and a combination of spectral and high order compact finite difference discretizations in the wall normal direction. The capabilities of the method are demonstrated with several examples including phase separation with, and without, shear in two and three dimensions. The method effectively resolves interfacial layers of as few as three mesh points. The numerical examples show agreement with analytical solutions and scaling laws, where available, and the 3D simulations, in the presence of shear, reveal rich and complex structures, including strings.

The faster growth curves in the speed of GPUs relative to CPUs in recent years and its rapidly gained popularity has spawned a new area of development in computational technology. There is much potential in utilizing GPUs for solving evolutionary partial differential equations and producing the attendant visualization. We are concerned with modeling tsunami waves, where computational time is of extreme essence, for broadcasting warnings. In order to test the efficacy of the GPU on the set of shallow-water equations, we employed the NVIDIA board 8600M GT on a MacBook Pro. We have compared the relative speeds between the CPU and the GPU on a single processor for two types of spatial discretization based on second-order finite-differences and radial basis functions. RBFs are a more novel method based on a gridless and a multi- scale, adaptive framework. Using the NVIDIA 8600M GT, we received a speed up factor of 8 in favor of GPU for the finite-difference method and a factor of 7 for the RBF scheme. We have also studied the atmospheric dynamics problem of swirling flows over a spherical surface and found a speed-up of 5.3 using the GPU. The time steps employed for the RBF method are larger than those used in finite-differences, because of the much fewer number of nodal points needed by RBF. Thus, in modeling the same physical time, RBF acting in concert with GPU would be the fastest way to go.

A semi-implicit finite-difference scheme is proposed for solving the nonlinear viscous compressible Navier-Stokes equations. Coordinate transformations are constructed that yield a uniform mesh in the computational plane even though the physical domain under consideration is time-varying and curvilinear. The finite-difference scheme was tested using model examples.

This paper presents a comparison of the results of numerical simulations obtained by two different numerical methods for one specific case of stably stratified incompressible flow. The focus in this paper is on the numerical results obtained using some of the compact finite-difference discretizations introduced in the paper [1]. The numerical scheme itself follows the principle of semi-discretisation, with high order compact discretisation in space, while the time integration is carried out by the Strong Stability Preserving Runge-Kutta scheme. Results are compared against the reference solution obtained by the AUSM finite volume method. The test case used to demonstrate the capabilities of selected numerical methods represents the flow of stably stratified fluid over low, smooth, hill-lik...

A systematic analysis shows how results from the finite difference code SEAWAT are sensitive to choice of grid dimension, time step, and numerical scheme for unstable flow problems. Guidelines to assist in selecting appropriate combinations of these factors are suggested. While the SEAWAT code has been tested for a wide range of problems, the sensitivity of results to spatial and temporal discretization levels and numerical schemes has not been studied in detail for unstable flow problems. Here, the Elder-Voss-Souza benchmark problem has been used to systematically explore the sensitivity of SEAWAT output to spatio-temporal resolution and numerical solver choice. A grid size of 0.38 and 0.60% of the total domain length and depth respectively is found to be fine enough to deliver results wi...

Convection in a porous medium may produce strong nonuniqueness of patterns. We study this phenomena for the case of a multicomponent fluid and develop a mimetic finite-difference scheme for the three-dimensional problem. Discretization of the Darcy equations in the primitive variables is based on staggered grids with five types of nodes and on a special approximation of nonlinear terms. This scheme is applied to the computer study of flows in a porous parallelepiped filled by a two-component fluid and with two adiabatic lateral planes. We found that the continuous family of steady stable states exists in the case of a rather thin enclosure. When the depth is increased, only isolated convective regimes may be stable. We demonstrate that the non-mimetic approximation of nonlinear terms leads...

Mandatory emission trading schemes are being established around the world. Participants of such market schemes are always exposed to risks. This leads to the creation of an accompanying market for emission-linked derivatives. To evaluate the fair prices of such financial products, one needs appropriate models for the evolution of the underlying assets, emission allowance certificates. In this paper, we discuss continuous time diffusion and jump-diffusion models, the latter enabling one to model information shocks that cause jumps in allowance prices. We show that the resulting martingale dynamics can be described in terms of non-linear partial differential and integro-differential equations and use a finite difference method to investigate numerical properties of their discretizations. The results are illustrated by a small numerical study.

We consider the numerical approximation of general semilinear parabolic stochastic partial differential equations (SPDEs) driven by additive space-time noise. In contrast to the standard time stepping methods which uses basic increments of the noise and the approximation of the exponential function by a rational fraction, we introduce a new scheme, designed for finite elements, finite volumes or finite differences space discretization, similar to the schemes in \\cite{Jentzen3,Jentzen4} for spectral methods and \\cite{GTambue} for finite element methods. We use the projection operator, the smoothing effect of the positive definite self-adjoint operator and linear functionals of the noise in Fourier space to obtain higher order approximations. We consider noise that is white in time and either in $H^1$ or $H^2$ in space and give convergence proofs in the mean square $L^{2}$ norm for a diffusion reaction equation and in mean square $ H^{1}$ norm in the presence of an advection term. For the exponential integrator...

The direct numerical simulation (DNS) of the Taylor-Couette flow in the fully turbulent regime is described. The numerical method extends the work by Quadrio and Luchini [M. Quadrio, P. Luchini, Eur. J. Mech. B/Fluids 21 (2002) 413-427], and is based on a parallel computer code which uses mixed spatial discretization (spectral schemes in the homogeneous directions, and fourth-order, compact explicit finite-difference schemes in the radial direction). A DNS is carried out to simulate for the first time the turbulent Taylor-Couette flow in the turbulent regime. Statistical quantities are computed to complement the existing experimental information, with a view to compare it to planar, pressure-driven turbulent flow at the same value of the Reynolds number. The main source for differences in ...

We discuss a numerical scheme to solve the continuum Kardar-Parisi-Zhang equation in generic spatial dimensions. It is based on a momentum-space discretization of the continuum equation and on a pseudospectral approximation of the nonlinear term. The method is tested in (1+1) and (2+1) dimensions, where it is shown to reproduce the current most reliable estimates of the critical exponents based on restricted solid-on-solid simulations. In particular, it allows the computations of various correlation and structure functions with high degree of numerical accuracy. Some deficiencies that are common to all previously used finite-difference schemes are pointed out and the usefulness of the present approach in this respect is discussed. PMID:11909192

A fast-slow factored scheme is presented for use with shallow-water primitive equation numerical weather prediction models. The technique was developed to reduce the rotational mode errors which arise when the fast and slow terms of the governing differential equations are treated simultaneously. The method factors out the fast and slow terms along the coordinate directions by means of a modified Crank-Nicolson scheme. A finite-difference spatial discretization is carried out in the zonal and meridional directions to reduce the factorization error to near-zero, and that time steps of 60-90 min can be used to obtain acceptably accurate results, even in the presence of fine spatial structures in the flow.

The macroscopic motion of a suspension of non-colloidal buoyant particles has been studied numerically. Phase separation takes place under the influence of centrifugal or gravitational force fields, resulting in a segmentation of the flow field into regions of pure fluid, mixture and sediment. A ``mixture model'', including collision terms, has been formulated in terms of volume averaged velocities. Efficient computer code to solve the equations in complex geometries has been developed by combining a variety of discretization schemes. Temporal splitting and projection methods are used to enhance the performance. The Stokes part of the momentum equation is spatially discretized using a Galerkin finite element method, and the advective part is discretized using a high-order upwind finite difference scheme. For the numerical treatment of discontinuities (kinematic shocks) associated with interfaces between regions of pure fluid, mixture and sediment, an upwind finite volume scheme has been developed. Implementations of the code support numerical simulations for arbitrary two-dimensional and axisymmetric geometries on unstructured grids. A few benchmark problems are presented to elucidate the ability of the code to handle a variety of phenomena occurring in basic mixture flows.

In this second report, the conservation properties of finite difference schemes in generalized curvilinear coordinate system constructed in the first report are examined by numerical tests of two-dimensional inviscid flow with periodic domains. It is confirmed that the modified finite difference scheme in a colocated grid layout conserves the mass and kinetic energy properly even in nonorth-ogonal nonuniform computational grid. The conservation properties of the finite difference scheme in a staggered grid layout are also proper only when computational grids are orthogonal. In addition, direct numerical simulations of a plane channel flow with fairly coarse mesh at R{sub e{tau}} = 180 are conducted using the finite difference schemes in generalized curvilinear coordinate system. It is found that the mesh dependency such as nonorthogonality and nonuniformity on results is reduced by the modified finite difference schemes in the colocated grid layout. (author)

Studies on applications of the finite element approach to transonic flow calculations are reported. Different discretization techniques of the differential equations and boundary conditions are compared. Finite element analogs of Murman's mixed type finite difference operators for small disturbance formulations were constructed and the time dependent approach (using finite differences in time and finite elements in space) was examined.

May 1, 1995 ... Title: Nonlinear wave propagation using three different finite difference ... used to examine the computation of one-dimensional nonlinear wave propagation. ... condition implementation, and discretization of governing equations. ... COMPUTATIONAL FLUID DYNAMICS; DAMPING; DISCONTINUITY; FINITE ...

A few families of counterexamples are provided to "A proof that the discrete singular convolution (DSC)/Lagrange-distributed approximation function (LDAF) method is inferior to high order finite differences", Journal of Computational Physics, 214, 538-549 (2006).

Finite difference schemes currently applied to the modeling of two-phase flows in flow networks exhibit difficulties in properly simulating certain spatial and temporal discontinuities. These discontinuities include points along the one-dimensional flow axis where density and other thermophysical properties become discontinuous or experience rapid state domain changes. A methodology for treating spatial and temporal discontinuities is presented. This methodology consists of three main features: (a) subnode time-averaged donoring of thermodynamic properties, (b) a variable pressure-at-discontinuity staggered mesh discretization, and (c) a variable point state equation linearization. The proposed scheme is similar in form to standard semi-implicit, staggered mesh discretizations, requires little extra overhead, and results in substantially improved accuracy and code execution times. Comparisons are made with standard time and spatial discretizations, as well as with two simpler alternate methods for recognizing and tracking discontinuities. The first of these attempts is to adjust the time-step size such that the fluid discontinuity arrives at a node boundary, or a change in fluid state occurs precisely at the end of a time advancement. The second attempts to redistribute mass and energy to correct for improperly donored values when a discontinuity crosses a node boundary during a time step. Neither of these alternatives proved adequate.

The connection between the solution of linear systems of equations by iterative methods and explicit time stepping techniques is used to accelerate to steady state the solution of ODE systems arising from discretized PDEs which may involve either physical or artificial transient terms. Specifically, a class of Runge-Kutta (RK) time integration schemes with extended stability domains has been used to develop recursion formulas which lead to accelerated iterative performance. The coefficients for the RK schemes are chosen based on the theory of Chebyshev iteration polynomials in conjunction with a local linear stability analysis. We refer to these schemes as Chebyshev Parameterized Runge Kutta (CPRK) methods. CPRK methods of one to four stages are derived as functions of the parameters which describe an ellipse {Epsilon} which the stability domain of the methods is known to contain. Of particular interest are two-stage, first-order CPRK and four-stage, first-order methods. It is found that the former method can be identified with any two-stage RK method through the correct choice of parameters. The latter method is found to have a wide range of stability domains, with a maximum extension of 32 along the real axis. Recursion performance results are presented below for a model linear convection-diffusion problem as well as non-linear fluid flow problems discretized by both finite-difference and finite-element methods.

A class of {open_quotes}fully-Lagrangian{close_quotes} methods for solving systems of conservation equations is defined. The key step in formulating these methods is the definition of a new set of field variables for which Lagrangian discretization is trivial. Recently popular lattice-Boltzmann simulation schemes for solving such systems are shown to be a useful sub-class of these fully-Lagrangian methods in which (a) the conservation laws are satisfied at each grid point, (b) the Lagrangian variables are expanded perturbatively, and (c) discretization error is used to represent physics. Such schemes are typically derived using methods of kinetic theory. Our numerical analysis approach shows that the conventional physical derivation, while certainly valid and fruitful, is not essential, that it often confuses physics and numerics and that it can be unnecessarily constraining. For example, we show that lattice-Boltzmann-like methods can be non-perturbative and can be made higher-order, implicit and/or with non-uniform grids. Furthermore, our approach provides new perspective on the relationship between lattice-Boltzmann methods and finite-difference techniques. Among other things, we show that the lattice-Boltzmann schemes are only conditionally consistent and in some cases are identical to the well-known Dufort-Frankel method. Through this connection, the lattice-Boltzmann method provides a rational basis for understanding Dufort-Frankel and gives a pathway for its generalization. At the same time, that Dufort-Frankel is no longer much used suggests that the lattice-Boltzmann approach might also share this fate.

It has long been recognized that a common ground exists between governing equations used for describing various flow and transport phenomena in porous media. Put another way they are all generally based on the same form of mass and/or energy conservation laws. This implies that there may exist a unified formulation and numerical scheme applicable to modeling all of these physical processes. This paper explores such a possibility and proposes a generalized framework, as well as a mathematical formulation for modeling all known transport phenomena in porous media. Based on this framework, a unified numerical approach is developed and tested using multidimensional, multiphase flow, isothermal and nonisothermal reservoir simulators. In this approach, a spatial domain of interest is discretized with an unstructured grid, then a time discretization is carried out with a backward, first-order, finite-difference method. The final discrete nonlinear equations are handled fully implicitly, using Newton iteration. In addition, the fracture medium is handled using a general dual-continuum concept with continuum or discrete modeling methods. A number of applications are discussed to demonstrate that with this unified approach, modeling a particular porous-medium flow and transport process simply becomes a matter of defining a set of state variables, along with their interrelations or mutual influence.

In this paper, a new smoothness indicator is proposed to improve the finite-difference lattice Boltzmann method (FDLBM). The necessary and sufficient conditions for convergence are derived. A detailed analysis reveals that the convergence order is higher than that of the previous finite-difference scheme. The coupled double distribution function (DDF) model is used to describe discontinuity flows and verify the improvement. Numerical simulations of compressible flows with shock wave show that the improved finite-difference lattice Boltzmann scheme is accurate and has less dissipation. The numerical results are found to be in good agreement with the analytical results and better than those of the previous scheme.

A computational method is developed involving the simultaneous integration of the Navier-Stokes and structural equations for the purpose of studying the stability of concentric annular passages conducting incompressible laminar flows. It is assumed that one side of the annulus, i.e. the centre-body, is fixed and the outer cylindrical duct is flexibly supported. The outer cylinder is displaced from its equilibrium position and is then released. In this situation, the fluid part of the problem is solved by an accurate method using a three-point backward implicit scheme, followed by a pseudo-time iteration using an artificial compressibility factor. The fluid equations are discretized in space based on a finite-difference formulation and primitive variables, for which stretched staggered grid...

A parallel fully implicit PETSc-based fluid modeling equations solver for simulating gas discharges is developed. Fluid modeling equations include: the neutral species continuity equation, the charged species continuity equation with drift-diffusion approximation for mass fluxes, the electron energy density equation, and Poisson's equation for electrostatic potential. Except for Poisson's equation, all model equations are discretized by the fully implicit backward Euler method as a time integrator, and finite differences with the Scharfetter???Gummel scheme for mass fluxes on the spatial domain. At each time step, the resulting large sparse algebraic nonlinear system is solved by the Newton???Krylov???Schwarz algorithm. A 2D-GEC RF discharge is used as a benchmark to validate our solver by...

A numerical method for the direct numerical simulation of incompressible wall turbulence in rectangular and cylindrical geometries is presented. The distinctive feature resides in its design being targeted towards an efficient distributed-memory parallel computing on commodity hardware. The adopted discretization is spectral in the two homogeneous directions; fourth-order accurate, compact finite-difference schemes over a variable-spacing mesh in the wall-normal direction are key to our parallel implementation. The parallel algorithm is designed in such a way as to minimize data exchange among the computing machines, and in particular to avoid taking a global transpose of the data during the pseudo-spectral evaluation of the non-linear terms. The computing machines can then be connected to each other through low-cost network devices. The code is optimized for memory requirements, which can moreover be subdivided among the computing nodes. The layout of a simple, dedicated and optimized computing system based ...

In this work we consider the numerical solution of a heat conduction problem for a material with non-constant properties. By approximating the time derivative of the solution through a finite difference, the transient equation is transformed into a sequence of inhomogeneous Helmholtz-type equations. The corresponding elliptic boundary value problems are then solved numerically by a meshfree method using fundamental solutions of the Helmholtz equation as shape functions. Convergence and stability of the method are addressed. Some of the advantages of this scheme are the absence of domain or boundary discretizations and/or integrations. Also, no auxiliary analytical or numerical methods are required for the derivation of the particular solution of the inhomogeneous elliptic problems. Numeric...

A two-dimensional lattice Boltzmann model with 19 discrete velocities for compressible Euler equations is proposed (D2V19-LBM). The fifth-order Weighted Essentially Non-Oscillatory (5th-WENO) finite difference scheme is employed to calculate the convection term of the lattice Boltzmann equation. The validity of the model is verified by comparing simulation results of the Sod shock tube with its corresponding analytical solutions. The velocity and density gradient effects on the Kelvin-Helmholtz instability (KHI) are investigated using the proposed model. Sharp density contours are obtained in our simulations. It is found that, the linear growth rate $\\gamma$ for the KHI decreases with increasing the width of velocity transition layer ${D_{v}}$ but increases with increasing the width of density transition layer ${D_{\\rho}}$. After the initial transient period and before the vortex has been well formed, the linear growth rates, $\\gamma_v$ and $\\gamma_{\\rho}$, vary with ${D_{v}}$ and ${D_{\\rho}}$ approximately i...

We report the results of a detailed study of the spectral properties of Laplace and Stokes operators, modified with a volume penalization term designed to approximate Dirichlet conditions in the limit when a penalization parameter, $\\eta$, tends to zero. The eigenvalues and eigenfunctions are determined either analytically or numerically as functions of $\\eta$, both in the continuous case and after applying Fourier or finite difference discretization schemes. For fixed $\\eta$, we find that only the part of the spectrum corresponding to eigenvalues $\\lambda \\lesssim \\eta^{-1}$ approaches Dirichlet boundary conditions, while the remainder of the spectrum is made of uncontrolled, spurious wall modes. The penalization error for the controlled eigenfunctions is estimated as a function of $\\eta$ and $\\lambda$. Surprisingly, in the Stokes case, we show that the eigenfunctions approximately satisfy, with a precision $O(\\eta)$, Navier slip boundary conditions with slip length equal to $\\sqrt{\\eta}$. Moreover, for a gi...

The immersed boundary method with discrete direct forcing approach was combined with a direct numerical simulation (DNS) code to study the time-dependent response of boundary-layer flow over a flat plate to the active vortex generators (AVG) that consist of a pair of deployable circular wing lip type blades driven by a pre-defined duty cycle. The DNS code solves the three-dimensional Navier-Stokes equations for compressible flow in general curvilinear coordinates using a fully implicit LU-SGS method. A fourth-order finite difference scheme is used to compute the spatial derivatives. The underlying curvilinear mesh near the vortex generator is designed properly with a large portion of the immersed boundary intersecting with the grid nodes to reduce the number of near-boundary nodes that req...

We put forward the use of total-variation-diminishing implicit-explicit Runge-Kutta methods for the time integration of the equations of motion associated with the semiconvection problem in the simulation of stellar convection. The fully compressible Navier-Stokes equation augmented by continuity and total energy equations and an equation of state describing the relation between the thermodynamic quantities is semi-discretized in space by essentially non-oscillatory schemes and dissipative finite difference methods, and subsequently integrated in time by Runge-Kutta methods which are constructed such as to reduce the total variation in the spatial profile in the course of time integration under certain restrictions on the time step-size. We analyze the stability, accuracy and dissipativity of the time integrators and demonstrate that the most successful methods yield a substantial gain in computational efficiency as compared to classical explicit Runge-Kutta methods.

Numerical optimization techniques combined with a three-dimensional thin-layer Navier-Stokes solver are presented to find an optimum shape of a stator blade in an axial compressor through calculations of single stage rotor-stator flow. Governing differential equations are discretized using an explicit finite difference method and solved by a multi-stage Runge-Kutta scheme. Baldwin-Lomax model is chosen to describe turbulence. A spatially-varying time-step and an implicit residual smoothing are used to accelerate convergence. A steady mixing approach is used to pass information between stator and rotor blades. For numerical optimization, searching direction is found by the steepest decent and conjugate direction methods, and the golden section method is used to determine optimum moving distance along the searching direction. The object of present optimization is to maximize efficiency. An optimum stacking line is found to design a custom-tailored 3-dimensional blade for maximum efficiency with the other parameters fixed.

We present finite difference schemes for solving the variable coefficient Poisson and heat equations on irregular domains with Dirichlet boundary conditions. The computational domain is discretized with non-graded Cartesian grids, i.e., grids for which the difference in size between two adjacent cells is not constrained. Refinement criteria is based on proximity to the irregular interface such that cells with the finest resolution is placed on the interface. We sample the solution at the cell vertices (nodes) and use quadtree (in 2D) or octree (in 3D) data structures as efficient means to represent the grids. The boundary of the irregular domain is represented by the zero level set of a signed distance function. For cells cut by the interface, the location of the intersection point is foun...

Although the numerical solution of one-dimensional phase-change, or Stefan, problems is well documented, a review of the most recent literature indicates that there are still unresolved issues regarding the start-up of a computation for a region that initially has zero thickness, as well as how to determine the location of the moving boundary thereafter. This paper considers the so-called boundary immobilization method for four benchmark melting problems, in tandem with three finite-difference discretization schemes. We demonstrate a combined analytical and numerical approach that eliminates completely the ad hoc treatment of the starting solution that is often used, and is numerically second-order accurate in both time and space, a point that has been consistently overlooked for this type...

A new algorithm for finding numerical solutions of optimal feedback control based on dynamic programming is developed. The algorithm is based on two observations: (1) the value function of the optimal control problem considered is the viscosity solution of the associated Hamilton-Jacobi-Bellman (HJB) equation and (2) the appearance of the gradient of the value function in the HJB equation is in the form of directional derivative. The algorithm proposes a discretization method for seeking optimal control-trajectory pairs based on a finite-difference scheme in time through solving the HJB equation and state equation. We apply the algorithm to a simple optimal control problem, which can be solved analytically. The consistence of the numerical solution obtained to its analytical counterpart in...

In this work the network simulation method (NSM) is proposed as a tool for solving the two-dimensional unsteady free convection and mass transfer flow of a viscous dissipative fluid along a semi-infinite vertical plate, taking into account the variation of the viscosity, thermal conductivity, and mass diffusivity with temperature. A spatial discretization, as in finite-difference schemes, is applied to the governing equations, with the time remains as a real continuous variable, and the different analogies (electrical-motion, electrical-thermal, and electrical-mass) developed in this work are applied to solve the problem in an adequate electric simulation program. The effect of variations in viscosity, thermal conductivity, and mass diffusivity with temperature is discussed. The effects of...

This paper presents a new numerical method, the Laplace Transform MultiQuadrics (LTMQ) method, developed for the solution of the diffusion-type parabolic Partial Differential Equation (PDE) of fluid flow through porous media. LTMQ combines a MultiQuadrics (MQ) approximation scheme with a Laplace transform formulation. The use of MQ in the spatial approximations allows the accurate description of problems in complex porous media with a very limited number of nodes. The Laplace transform formulation eliminates the need for time discretization, thus allowing an unlimited time step size without any loss of accuracy. LTMQ is tested against results from three test problems of groundwater flow obtained from a standard Finite Difference (FD) model, as well as from analytical solutions. An excellent agreement between the LTMQ and the analytical and FD solutions is observed, while significant reductions in computer execution times may be achieved.

In this work, a nonlinear model to predict actuation characteristics in lateral electrostatically-actuated DC-contact MEMS switches is proposed. In this case a parallel-plate approximation cannot be applied. The model is based on the equilibrium equation for an elastic beam. It is demonstrated that the contribution of fringing fields is essential. The model relies on finite-difference discretization of the structures, applying boundary conditions and is solved with a Gauss-Seidel relaxation iteration scheme. Its usefulness is demonstrated in a series MEMS switch with lateral interdigital electrostatic actuation. Copyright 2007 Wiley Periodicals, Inc. Microwave Opt Technol Lett 49: 1238-1241, 2007; Published online in Wiley InterScience (www.interscience.wiley.com). DOI 10.1002/mop.22450

We present a hierarchical transform that can be applied to Laplace-like differential equations such as Darcy’s equation for single-phase flow in a porous medium. A finite-difference discretization scheme is used to set the equation in the form of an eigenvalue problem. Within the formalism suggested, the pressure field is decomposed into an average value and fluctuations of different kinds and at different scales. The application of the transform to the equation allows us to calculate the unknown pressure with a varying level of detail. A procedure is suggested to localize important features in the pressure field based only on the fine-scale permeability, and hence we develop a form of adaptive coarse graining. The formalism and method are described and demonstrated using two synthetic toy problems.

Fluid-solid interaction is a difficult computational problem, primarily because solids and fluids are described in different perspectives --- solid laws are written in Lagrangian frame while fluids are represented in Eulerian. Our work attempts to resolve this dilemma using a new method for Eulerian solid mechanics. We study the interaction of a large-deformation elastic solid with a Newtonian fluid in a single computational framework. We use a level set to track the interface between the two phases. The standard projection method is used to impose incompressibility in both phases, and the equations are discretized with an explicit, staggered finite-difference scheme. In the current implementation, a smeared Heaviside function is used to blur material properties across the interface. Simulations of various test cases will be presented in this talk.

Finite element and finite difference methods are examined in order to bring out their relationship. It is shown that both methods use two types of discrete representations of continuous functions. They differ in that finite difference methods emphasize the discretization of independent variable, while finite element methods emphasize the discretization of dependent variable (referred to as functional approximations). An important point is that finite element methods use global piecewise functional approximations, while finite difference methods normally use local functional approximations. A general conclusion is that finite element methods are best designed to handle complex boundaries, while finite difference methods are superior for complex equations. It is also shown that finite volume difference methods possess many of the advantages attributed to finite element methods.

This short note reports the extension of the f-waves approximate Riemann solver (A hmad and L indeman, 2007; L eV eque, 2002; B ale et al., 2002) for three-dimensional meso- and micro-scale atmospheric flows. The Riemann solver employs flux-based wave decomposition for the calculation of Godunov fluxes and does not require the explicit definition of the Roe matrix to enforce conservation. The other important feature of the Riemann solver is its ability to incorporate source term due to gravity without introducing discretization errors. The resulting finite volume scheme is second-order accurate in space and time. The finite-difference schemes currently used in atmospheric flow models are neither conservative nor able to resolve regions of sharp gradients. The finite volume scheme described in this paper is fully conservative and has the ability to resolve regions of sharp gradients without introducing spurious oscillations in the solution. The scheme shows promise in accurately resolving flows on the meso- and micro-scales and should be considered for implementation in the dynamical cores of next generation meso- and micro-scale atmospheric flow models.

We describe a numerical scheme for solving Maxwell's equations in the frequency domain on a conformal, structured, non-orthogonal, multi-block mesh. By considering Maxwell's equations in a volume parameterized by dimensionless curvilinear coordinates, we obtain a set of tensor equations that are a continuum analogue of common circuit equations, and that separate the metrical and metric-free parts of Maxwell's equations and the material constitutive relations. We discretize these equations using a new formulation that treats the electric field and magnetic induction using simple basis-function representations to obtain a discrete form of Faraday's law of induction, but that uses finite integral representations for the displacement current and magnetic field to obtain a discrete form of Ampere's law, as in the finite integration technique [T. Weiland, A discretization method for the solution of Maxwell's equations for six-component fields, Electron. Commun. (AE U) 31 (1977) 116; T. Weiland, Time domain electromagnetic field computation with finite difference methods, Int. J. Numer. Model: Electron. Netw. Dev. Field 9 (1996) 295-319]. We thereby derive new projection operators for the discrete tensor material equations and obtain a compact numerical scheme for the discrete differential operators. This scheme is shown to exhibit significantly reduced numerical dispersion when compared to the standard linear finite element method. We take advantage of the mesh structure on a block-by-block basis to implement these numerical operators efficiently, and achieve computational speed with modest memory requirements when compared to explicit sparse matrix storage. Using the Jacobi-Davidson [G.L.G. Sleijpen, H.A. van der Vorst, A Jacobi-Davidson iteration method for linear eigenvalue problems, SIAM J. Matrix Anal. Appl. 17 (2) (1996) 401-425; S.J. Cooke, B. Levush, Eigenmode solution of 2-D and 3-D electromagnetic cavities containing absorbing materials using the Jacobi-Davidson algorithm, J. Comput. Phys. 157 (1) (2000) 350-370] and quasi-minimal residual [R.W. Freund, N.M. Nachtigal, QMR: a quasi-minimal residual method for non-Hermitian linear systems, Numer. Math. 60 (1991) 315-339] iterative matrix solution algorithms, we solve the resulting discrete matrix eigenvalue equations and demonstrate the convergence characteristics of the algorithm. We validate the model for three-dimensional electromagnetic problems, both cavity eigenvalue solutions and a waveguide scattering matrix calculation.

In engineering fields there have recently been a number of studies that have applied the finite-difference time-domain (FDTD) numerical method to elastodynamic problems. These studies include those of ultrasonic transducers and electro-mechanical devices. In FDTD formulations, the first-order set of partial differential equations given by the constitutive equations are discretized using a leap-frog finite-difference scheme. When a high-contrast discontinuity, especially a free surface, is present, some difficulties arise due to the spatially staggered nature of the grid of the FDTD approach, with neither all the velocity variables nor all the stress variables appearing on the same grid lines. The present study considers the following modifications to the FDTD approach: a standard staggered grid for anisotropic elastic wave fields is rotated so that the diagonal directions of the standard grid lie parallel to the axes of the analysis region. This configuration, called the diagonally staggered grid (DSG), improves the accuracy of the implementation of free boundaries without requiring virtual grids in a vacuum area. The effectiveness of DSG was verified by applying this method to model problems of isotropic and anisotropic solid materials.

Galerkin finite element methods. IResearch partly supported by the National. Science Foundation ...... finite difference. ENO schemes and finite element discontinuous. Galerkin methods, see [17], [5] and [6]. .... Notes in Physics, v371, Springer- ...

charge cup. The high thermal conductivity of alumina at cryogenic temperatures was overlooked in ... The objective of this test was to determine if the oxidation level of zirconium affected the .... We employed an explicit finite difference scheme ...

and Integral. Equations. Lafon, F. and S. Osher: High order filtering methods for approximating .... of implicit finite difference schemes in computational fluid dynamics. ICASE. Report. No. ...... of the Korteweg- de Vries equation. ICASE. Report ...

An efficIent numerical mesh generation scheme capable of creatmg orthog- onal or nearly orthogonal grids about moderately complex three dlmensional con- ... Body conforming curvilinear grids are often used m finite difference flo\\\\.

In this work, we present a numerical method to consistently approximate solutions of a spatially discrete, double sine-Gordon chain which considers the presence of external damping. In addition to the finite-difference scheme employed to approximate the solution of the difference-differential equations of the model under investigation, our method provides positivity-preserving schemes to approximate the local and the total energy of the system, in such a way that the discrete rate of change of the total energy with respect to time provides a consistent approximation of the corresponding continuous rate of change. Simulations are performed, first of all, to assess the validity of the computational technique against known qualitative solutions of coupled sine-Gordon and coupled double sine-Gordon chains. Secondly, the method is used in the investigation of the phenomenon of nonlinear transmission of energy in double sine-Gordon systems; the qualitative effects of the damping coefficient on the occurrence of the nonlinear process of supratransmission are briefly determined in this work, too.

The compressible Navier Stokes equations can be extended to model multi-species, chemically reacting gas flows. The result is a large system of convection-diffusion equations with stiff source terms. In this paper we develop the framework needed to apply modern high accuracy numerical methods from computational gas dynamics to this extended system. We also present representative computational results using one such method. The framework developed here is useful for many modern numerical schemes. We first present an enthalpy based form of the equations that is well suited both for physical modeling and for numerical implementation. We show how to treat the stiff reactions via time splitting, and in particular how to increase accuracy by avoiding the common practice of approximating the temperature. We derive simple, exact formulas for the characteristics of the convective part of the equations, which are essential for application of all characteristic-based schemes. We also show that the common practice of using approximate analytical expressions for the characteristics can potentially produce spurious oscillations in computations. We implement these developments with a particular high accuracy characteristic-based method, the finite difference ENO space discretization with the 3rd order TVD Runge-Kutta time discretization, combined with the second order accurate Strang time splitting of the reaction terms. We illustrate the capabilities of this approach with calculations of a 1-D reacting shock tube and a 2-D combustor. 16 refs., 11 figs.

Optical tomography is a medical imaging technique based on light propagation in the near infrared (NIR) part of the spectrum. We present a new way of predicting the short-pulsed NIR light propagation using a time-dependent two-dimensional-global radiative transfer equation in an absorbing and strongly anisotropically scattering medium. A cell-vertex finite-volume method is proposed for the discretization of the spatial domain. The closure relation based on the exponential scheme and linear interpolations was applied for the first time in the context of time-dependent radiative heat transfer problems. Details are given about the application of the original method on unstructured triangular meshes. The angular space (4?Sr) is uniformly subdivided into discrete directions and a finite-differences discretization of the time domain is used. Numerical simulations for media with physical properties analogous to healthy and metastatic human liver subjected to a collimated short-pulsed NIR light are presented and discussed. As expected, discrepancies between the two kinds of tissues were found. In particular, the level of light flux was found to be weaker (inside the medium and at boundaries) in the healthy medium than in the metastatic one. PMID:22894479

In this Letter, a three-dimensional (3D) lattice-Boltzmann model is presented following the non-free-parameter lattice-Boltzmann method of Qu et al. [K. Qu, C. Shu, Y.T. Chew, Phys. Rev. E 75 (2007) 036706]. A simple function, which satisfies the zeroth- through third-order moments of the Maxwellian distribution function, is introduced to replace the Maxwellian distribution function as the continuous equilibrium distribution function in 3D space. The function is then discretized to discrete-velocity directions via a 25-point Lagrangian interpolation polynomial. To simulate compressible flows with shock waves, an implicit-explicit finite-difference scheme based on the total variation diminishing flux limitation is adopted to solve the discrete Boltzmann-BGK equation in order to capture the shock waves in compressible flows with a finite number of grid points. The model is validated by its application to some typical inviscid compressible flows ranging from 1D to 3D, and the numerical results are found to be in excellent agreement with the analytical solutions and/or other numerical results.

A finite?difference method for integro?differential equations arising from Lévy driven asset processes in finance is discussed. The equations are discretized in space by the collocation method and in time by an explicit backward differentiation formula. The discretization is shown to be second?order...

Vortex methods have a history as old as finite differences. They have since faced difficulties stemming from the numerical complexity of the Biot-Savart law, the inconvenience of adding viscous effects in a Lagrangian formulation, and the loss of accuracy due to Lagrangian distortion of the computational elements. The first two issues have been successfully addressed, respectively, by the application of the fast multipole method, and by a variety of viscous schemes which will be briefly reviewed in this article. The standard method to deal with the third problem is the use of remeshing schemes consisting of tensor product interpolation with high-order kernels. In this work, a numerical study of the errors due to remeshing has been performed, as well as of the errors implied in the discretization itself using vortex blobs. In addition, an alternative method of controlling Lagrangian distortion is proposed, based on ideas of radial basis function (RBF) interpolation (briefly reviewed here). This alternative is formulated grid-free, and is shown to be more accurate than standard remeshing. In addition to high-accuracy, RBF interpolation allows core size control, either for correcting the core spreading viscous scheme or for providing a variable resolution in the physical domain. This formulation will allow in theory the application of error estimates to produce a truly adaptive spatial refinement technique. Proof-of-concept is provided by calculations of the relaxation of a perturbed monopole to a tripole attractor.

Abstract in english The present work is inserted into an effort to develop a Chimera flow simulation code capable of handling general launch vehicle configurations. The paper is primarily concerned with presenting results of laminar and turbulent viscous simulations of flows over the first Brazilian satellite launch vehicle, the VLS, during its first-stage flight. The finite difference method is applied to the governing equations written in conservation-law form for general body conforming c (more) urvilinear coordinates. The spatial discretization is accomplished with a central difference scheme in which artificial dissipation terms, based on a scalar, non-isotropic model, are added to the numerical scheme to maintain stability. The time march process is accomplished with a 5-stage, 2nd-order accurate, Runge-Kutta scheme. The results here included are indicative of the current status of the Chimera flow simulation capability under development by the authors. The results also highlight interesting features of the flow over the complete VLS and point out the importance of the inclusion of viscous effects for flow simulation over such complex vehicles.

In this paper, several discrete schemes for Gaussian curvature are surveyed. The convergence property of a modified discrete scheme for the Gaussian curvature is considered. Furthermore, a new discrete scheme for Gaussian curvature is resented. We prove that the new scheme converges at the regular vertex with valence not less than 5. By constructing a counterexample, we also show that it is impossible for building a discrete scheme for Gaussian curvature which converges over the regular vertex with valence 4. Finally, asymptotic errors of several discrete scheme for Gaussian curvature are compared.

A 3-D hydrodynamic dispersion model for tracer transport is developed and implemented into the TOUGH2 EOS3 (T2R3D) module. The model formulation incorporates a full dispersion tensor, based on a 3-D velocity field with a 3-D, irregular grid in a heterogeneous geological system. Two different weighting schemes are proposed for spatial average of 3-D velocity fields and concentration gradients to evaluate the mass flux by dispersion and diffusion of a tracer or a radionuclide. This new module of the TOUGH2 code is designed to simulate processes of tracer/radionuclide transport using an irregular, 3-D integral finite difference grid in non-isothermal, three-dimensional, multiphase, porous/fractured subsurface systems. The numerical method for this transport module is based on the integral finite difference scheme, as in the TOUGH2 code. The major assumptions of the tracer transport module are: (a) a tracer or a radionuclide is present and transported only within the liquid phase, (b) transport mechanisms include molecular diffusion and hydrodynamic dispersion in the liquid phase in addition to advection, and (c) first order decay and linear adsorption on rock grains are taken into account. The tracer or radionuclide is introduced as an additional mass component into the standard TOUGH2 formulation, time is discretized fully implicitly, and non-linearities of the conservation equations are handled using the Newton/Raphson iteration. We have verified this transport module by comparison with results of a 2-D transport problem for which an analytical solution is available. In addition, a field application is described to demonstrate the use of the proposed model.

High accuracy solution of PDEs requires proper error analysis. Previous analysis for a non-dispersive system [T.K. Sengupta, A. Dipankar, P. Sagaut, Error dynamics: beyond von Neumann analysis, J. Comput. Phys. 226 (2007) 1211-1218] identified sources of error correctly. Here, the aim is to extend the spectral analysis for the model linearized rotating shallow water equations (LRSWE), as an example of dispersive system. We perform the analysis when high accuracy compact schemes are used to solve the LRSWE relevant to geophysical fluid dynamics, using different grid arrangements proposed in Mesinger and Arakawa [F. Mesinger, A. Arakawa, Numerical Methods Used in Atmospheric Models, GARP Publ. Ser. No. 17, vol. 1, WMO, Geneva, 1976, pp. 43-64] and Randall [D.A. Randall, Geostrophic adjustment and the finite-difference shallow-water equations, Mon. Wea. Rev. 122 (1994) 1371-1377]. Compact schemes are used for fluid dynamical problem, as these afford near-spectral accuracy in solving non-periodic problems. However, higher accuracy methods also suffer from errors, those are often filtered by low order methods. For example, dispersion error is present in all numerical methods and extreme form of it leads to q-waves, which appear at higher wavenumbers for compact schemes as compared to lower order method. We also evaluate a compact scheme specifically designed for use with staggered grids. Here, two and four time-level temporal discretization methods have been compared for solving LRSWE by considering classical fourth-order, four-stage Runge-Kutta (RK4), two time-level forward-backward (FB) and four time-level generalized FB temporal integration schemes.

As part of the continuing effort at NASA LeRC to improve both the durability and reliability of hot section Earth-to-orbit engine components, significant enhancements must be made in existing finite element and finite difference methods, and advanced techniques, such as the boundary element method (BEM), must be explored. The BEM was chosen as the basic analysis tool because the critical variables (temperature, flux, displacement, and traction) can be very precisely determined with a boundary-based discretization scheme. Additionally, model preparation is considerably simplified compared to the more familiar domain-based methods. Furthermore, the hyperbolic character of high speed flow is captured through the use of an analytical fundamental solution, eliminating the dependence of the solution on the discretization pattern. The price that must be paid in order to realize these advantages is that any BEM formulation requires a considerable amount of analytical work, which is typically absent in the other numerical methods. All of the research accomplishments of a multi-year program aimed toward the development of a boundary element formulation for the study of hot fluid-structure interaction in Earth-to-orbit engine hot section components are detailed. Most of the effort was directed toward the examination of fluid flow, since BEM's for fluids are at a much less developed state. However, significant strides were made, not only in the analysis of thermoviscous fluids, but also in the solution of the fluid-structure interaction problem.

In the past, arguments have been advanced suggesting that certain finite-difference solutions of the 3+1 form of Einstein's equations suffer from a fundamental inconsistency. Specifically, it has been claimed that freely evolved solutions, where the constraint equations are not explicitly imposed after the initial time, will generally satisfy discrete versions of the constraints to a lower order in the basic mesh spacing {ital h} than the truncation order of the discretized evolution equations. This issue is reexamined here, and using the key observation, originally due to Richardson, that a numerical differentiation need not produce an {ital O}({ital h}{sup {ital p}{minus}1}) quantity from an {ital O}({ital h}{sup {ital p}}) one, it is argued that there should be {ital no} such inconsistency for convergent difference schemes. Numerical results from a study of spherically symmetric solutions of a massless scalar field minimally coupled to the gravitational field are presented in support of this claim. These results show that the expected convergence of various residual quantities can be achieved in practice.

We present the development and application of compatible finite element discretizations of electromagnetics problems derived from the time dependent, full wave Maxwell equations. We review the H(curl)-conforming finite element method, using the concepts and notations of differential forms as a theoretical framework. We chose this approach because it can handle complex geometries, it is free of spurious modes, it is numerically stable without the need for filtering or artificial diffusion, it correctly models the discontinuity of fields across material boundaries, and it can be very high order. Higher-order H(curl) and H(div) conforming basis functions are not unique and we have designed an extensible C++ framework that supports a variety of specific instantiations of these such as standard interpolatory bases, spectral bases, hierarchical bases, and semi-orthogonal bases. Virtually any electromagnetics problem that can be cast in the language of differential forms can be solved using our framework. For time dependent problems a method-of-lines scheme is used where the Galerkin method reduces the PDE to a semi-discrete system of ODE's, which are then integrated in time using finite difference methods. For time integration of wave equations we employ the unconditionally stable implicit Newmark-Beta method, as well as the high order energy conserving explicit Maxwell Symplectic method; for diffusion equations, we employ a generalized Crank-Nicholson method. We conclude with computational examples from resonant cavity problems, time-dependent wave propagation problems, and transient eddy current problems, all obtained using the authors massively parallel computational electromagnetics code EMSolve.

A numerically stable method to solve the discretized Boltzmann-Enskog equation describing the behavior of nonideal fluids under inhomogeneous conditions is presented. The algorithm employed uses a Lagrangian finite-difference scheme for the treatment of the convective term and a forcing term to account for the molecular repulsion together with a Bhatnagar-Gross-Krook relaxation term. In order to eliminate the spurious currents induced by the numerical discretization procedure, we use a trapezoidal rule for the time integration together with a version of the two-distribution method of He [J. Comput. Phys.JCTPAH0021-999110.1006/jcph.1999.6257 152, 642 (1999)]. Numerical tests show that, in the case of a one-component fluid in the presence of a spherical potential well, the proposed method reduces the numerical error by several orders of magnitude. We conduct another test by considering the flow of a two-component fluid in a channel with a bottleneck and provide information about the density and velocity field in this structured geometry.

A third-order ordinary differential equation (ODE) for thin film flow with both Neumann and Dirichlet boundary conditions is transformed into a second-order nonlinear ODE with Dirichlet boundary conditions. Numerical solutions of the nonlinear second-order ODE are investigated using finite difference schemes. A finite difference formulation to an Emden-Fowler representation of the second-order nonlinear ODE is shown to converge faster than a finite difference formulation of the standard form of the second-order nonlinear ODE. Both finite difference schemes satisfy the von Neumann stability criteria. When mapping the numerical solution of the second-order ODE back to the variables of the original third-order ODE we recover the position of the contact line. A nonlinear relationship between the position of the contact line and physical parameters is obtained.

Abstract in english Mixing layers are present in very different types of physical situations such as atmospheric flows, aerodynamics and combustion. It is, therefore, a well researched subject, but there are aspects that require further studies. Here the instability of two-and three-dimensional perturbations in the compressible mixing layer was investigated by numerical simulations. In the numerical code, the derivatives were discretized using high-order compact finite-difference schemes. A (more) stretching in the normal direction was implemented with both the objective of reducing the sound waves generated by the shear region and improving the resolution near the center. The compact schemes were modified to work with non-uniform grids. Numerical tests started with an analysis of the growth rate in the linear regime to verify the code implementation. Tests were also performed in the non-linear regime and it was possible to reproduce the vortex roll-up and pairing, both in two-and three-dimensional situations. Amplification rate analysis was also performed for the secondary instability of this flow. It was found that, for essentially incompressible flow, maximum growth rates occurred for a spanwise wavelength of approximately 2/3 of the streamwise spacing of the vortices. The result demonstrated the applicability of the theory developed by Pierrehumbet and Widnall. Compressibility effects were then considered and the maximum growth rates obtained for relatively high Mach numbers (typically under 0.8) were also presented.

This study enables the use of very high-order finite-difference schemes for the solution of conservation laws on stretched, curvilinear, and deforming meshes. To illustrate these procedures, we focus on up to 6th-order Pade-type spatial discretizations coupled with up to 10th-order low-pass filters. These are combined with explicit and implicit time integration methods to examine wave propagation and wall-bounded flows described by the Navier-Stokes equations. It is shown that without the incorporation of the filter, application of the high-order compact scheme to nonsmooth meshes results in spurious oscillations which inhibit their applicability. Inclusion of the discriminating low-pass high-order filter restores the advantages of high-order approach even in the presence of large grid discontinuities. When three-dimensional curvilinear meshes are employed, the use of standard metric evaluation procedures significantly degrades accuracy since freestream preservation is violated. To overcome this problem, a simple technique is adopted which ensures metric cancellation and thus ensures freestream preservation even on highly distorted curvilinear meshes. For dynamically deforming grids, an effective numerical treatment is described to evaluate expressions containing the time-varying transformation metrics. With these techniques, metric cancellation is guaranteed regardless of the manner in which grid speeds are defined. The efficacy of the new procedures is demonstrated by solving several model problems as well as by application to flow past a rapidly pitching airfoil and past a flexible panel.

Numerical method is presented for simulation of the deformation, drainage and rupture of axisymmetric film (gap) between colliding drops in the presence of film soluble surfactants under the influence of van der Waals forces at small capillary and Reynolds numbers and small surfactant concentrations. The mathematical model is based on the lubrication equations in the gap between drops and the creeping flow approximation of Navier-Stokes equations in the drops, coupled with velocity and stress boundary conditions at the interfaces. A non-uniform surfactant concentration on the interfaces, related with that in the film, leads to a gradient of the interfacial tension which in turn leads to additional tangential stress on the interfaces (Marangoni effects). Both film and interface surfactant concentrations, related via adsorption isotherm, are governed by a convection-diffusion equation. The numerical method consists of: Boundary integral method for the flow in the drops; Finite difference method for the flow in the gap, the position of the interfaces and the surfactant concentration on the interfaces, as well as in the film. Second order approximation of the spatial terms on adaptive non-uniform mesh is constructed in combination with Euler explicit scheme for the time discretization. For the convection-diffusion equation in the film first order implicit and Crank-Nicolson time integration schemes are used as well. Tests and comparisons are performed to show the accuracy and stability of the presented numerical method.

In this paper we develop a finite-difference scheme to approximate radially symmetric solutions and (1 + 1)-dimensional solutions of the initial-value problem with smooth initial conditions 62w6t2 -12w-b 66t1 2w+g6w 6t+m2w+G' w=0, subjectto: wx¯,0 =fx¯ ,x¯?D, 6w6t x¯,0+y x¯, x¯?D, in an open sphere D around the origin, where the internal and external damping coefficients beta and gamma, respectively, are constant. The functions ? and psi are radially symmetric in D, they are small at infinity, and r?( r) and rpsi(r) are also assumed to be small at infinity. We prove that our scheme is consistent order O (Deltat2) + O (Deltar2) for G ' identically equal to zero, and provide a necessary condition for it to be stable order n. A cornerstone of our investigation will be the study of potential applications of our model to discrete versions involving nonlinear systems of coupled oscillators. More concretely, we make use of the process of nonlinear supratransmission of energy in these chain systems and our numerical techniques in order to transmit binary information. Our simulations show that, under suitable parametric conditions, the transmission of binary signals can be achieved successfully.

Galerkin and Petrov-Galerkin finite element methods are used to obtain new finite difference schemes for the solution of linear two-dimensional convection-diffusion problems. Numerical estimates are made of the rates of convergence of these schemes, uniformly with respect to the perturbation parameter, and these uniform rates are shown to compare favourably with those of established methods.

The aim of this work is to develop general optimization methods for finite difference schemes used to approximate linear differential equations. The specific case of the transport equation is exposed. In particular, the minimization of the numerical error is taken into account. The theoretical study of a related linear algebraic problem gives general results which can lead to the determination of the optimal scheme.

The recently developed exponential discontinuous spatial differencing scheme for the discrete-ordinate equations has been extended to x-y-z geometry with hexahedral cells. This scheme produces strictly positive angular fluxes given positive discrete-ordinate sources. The exponential discontinuous scheme has been developed and implemented into the three-dimensional, discrete-ordinate code. THREEDANT. Numerical results are given which show that the exponential discontinuous scheme is very accurate for deep-penetration transport problems with optically thick spatial meshes.

The numerical solution of multidimensional wave-propagation problems is considerably more complex than solutions for one-dimensional problems, but improving a method`s accuracy produces more significant increases in efficiency. This is particularly true for first-order accurate methods. Higher-order methods, which are already comparatively efficient, are especially difficult to construct on irregular meshes. The passive-scalar advection equation provides an ideal vehicle for investigating the relationships between accuracy, efficiency, and physical acceptability for a variety of finite-difference schemes. Development and analysis of methods for this equation provide insight into solving more complicated problems involving material and/or neutron transport. We found that three approaches to two-dimensional upwind differencing led to schemes equivalent to the dimensionally split Lax-Wendroff method. Assuming stability, we were able to prove that both the split-Lax-Wendroff and a new two-dimensional, predictor-corrector scheme yielded second-order convergence on nonuniform tensor-product grids, despite first-order discretization errors for such grids. We were able to prove that in some cases, Roe`s upwinding for triangulations (which is locally inconsistent) converges with first-order accuracy. The sensitivity of fourth-order accurate, conservative, mimetic methods to the roughness of nonuniform grid-spacing was explored analytically and numerically in one-dimension and for logically-rectangular two-dimensional grids. We also considered the application of a novel set of unknowns towards conservative second-order methods on irregular polyhedral grids. Queried, we noted that convergent, stable, split schemes for ill-posed problems could not exist. This prompted us to reconsider conditions implying that splitting is stable for the differential operators and for related difference methods.

Simulating any realistic seismic scenario requires incorporating physical basis into the model. Considering both the dynamics of the rupture process and the anelastic attenuation of seismic waves is essential to this purpose and, therefore, we choose to extend the hp-adaptive Discontinuous Galerkin finite-element method to integrate these physical aspects. The 3D elastodynamic equations in an unstructured tetrahedral mesh are solved with a second-order time marching approach in a high-performance computing environment. The first extension incorporates the viscoelastic rheology so that the intrinsic attenuation of the medium is considered in terms of frequency dependent quality factors (Q). On the other hand, the extension related to dynamic rupture is integrated through explicit boundary conditions over the crack surface. For this visco-elastodynamic formulation, we introduce an original discrete scheme that preserves the optimal code performance of the elastodynamic equations. A set of relaxation mechanisms describes the behavior of a generalized Maxwell body. We approximate almost constant Q in a wide frequency range by selecting both suitable relaxation frequencies and anelastic coefficients characterizing these mechanisms. In order to do so, we solve an optimization problem which is critical to minimize the amount of relaxation mechanisms. Two strategies are explored: 1) a least squares method and 2) a genetic algorithm (GA). We found that the improvement provided by the heuristic GA method is negligible. Both optimization strategies yield Q values within the 5% of the target constant Q mechanism. Anelastic functions (i.e. memory variables) are introduced to efficiently evaluate the time convolution terms involved in the constitutive equations and thus to minimize the computational cost. The incorporation of anelastic functions implies new terms with ordinary differential equations in the mathematical formulation. We solve these equations using the same order of interpolation as for the elastic equations (i.e. the so-called P0, P1 or P2 interpolations functions). We compare solutions from several numerical strategies (e.g. Finite Difference and Discontinuous Galerkin methods). For the second extension, the dynamic rupture formulation requires explicit boundary conditions on discontinuous surface edges bounding the fracture. These conditions have been implemented for the different interpolation orders we consider and are based on the conservation of a discrete energy. The fault shear stress follows a linear slip-weakening law, although any other friction law could be implemented. We validate our mathematical and computational model by comparing synthetic seismograms with those yielded by the semi-analytical Discrete Wavenumber method for the attenuation effect and with a well-verified Finite Difference method (SGSN) for the dynamic rupture model. This work will allow us shortly to perform realistic simulations of possible physics-based seismic scenarios in the Valley of Mexico to study the associated hazard.

Integro-differential equations of Volterra type arise, naturally, in many applications such as for instance heat conduction in materials with memory, diffusion in polymers and diffusion in porous media. The aim of this paper is to study a finite difference discretization of the mentioned integro-differential equations. Second convergence order with respect to the Formula Not Shown norm is established which means that the discretization proposed is supraconvergent in finite difference methods language. As the finite difference method can be seen as a piecewise linear finite element method combined with special quadrature formulas, our result establishes the supercloseness of the gradient in the finite element language. Numerical results illustrating the discussed theoretical results are inc...

Nonlocal conditions arise in mathematical models of various physical, chemical or biological processes. Therefore, interest in developing computational techniques for the numerical solution of partial differential equations (PDEs) with various types of nonlocal conditions has been growing fast. We construct and analyse a weighted splitting finite-difference scheme for a two-dimensional parabolic equation with nonlocal integral conditions. The main attention is paid to the stability of the method. We apply the stability analysis technique which is based on the investigation of the spectral structure of the transition matrix of a finite-difference scheme. We demonstrate that depending on the parameters of the finite-difference scheme and nonlocal conditions the proposed method can be stable ...