Numerical Conservation Properties

In this section the numerical conservation properties of finite volume schemes for inviscid incompressible flow are examined. Emphasis is placed on the theory of kinetic energy conservation. Numerical issues associated with the use of kinetic energy non-conservative schemes are discussed [158, 49, 47]. 

The early meteorological finite difference studies of long-term numerical time integrations of the equations of fluid motion, which involve non-linear convection terms, revealed the presence of non-linear instabilities due to aliasing errors [143, 144, 7,145, 210]. To avoid the occurrence of these non-linear instabilities, Arakawa [7] was the first to recognize the importance of the use of numerical schemes which conserve kinetic energy. 

If the local kinetic energy equation (1.133), derived in Sect. 1.2.4, is integrated over a grid volume, and transforming some of the resulting terms from volume to surface integrals by the Gauss theorem, the result is  

The first three terms on the RHS are integrals over the surface of the macroscopic volume. It is noted that the kinetic energy in the macroscopic volume is not changed by the action of convection and pressure. The fourth term on the RHS is an inte- 

It is noted that kinetic energy conservation is associated with numerical stability but not to the convergence or accuracy properties of a method. However, the physical dissipation used in the simulation is determined by the sum of the dynamic model dissipation and the apparent dissipation added. 

---

For long-term simulations, it generally proves advantageous to consider numerical integrators which pass the structural properties of the model onto the calculated solutions. Hence, a careful analysis of the conservation properties of QCMD model is required. A particularly relevant constant of motion of the QCMD model is the total energy of the system... 

This discretization method obeys a conservation property, and therefore is called conservative. With the exception of the first element and the last element, every element face is a part of two elements. The areas of the coincident faces and the forces on them are computed in exactly the same way (except possibly for sign). Note that the sign conventions for the directions of the positive stresses is important in this regard. The force on the left face of some element is equal and opposite to the force on the right face of its leftward neighbor. Therefore, when the net forces are summed across all the elements, there is exact cancellation except for the first and last elements. For this reason no spurious forces can enter the system through the numerical discretization itself. The net force on the system of elements must be the net force caused by the boundary conditions on the left face of the first element and the right face of the last element. 

Equation (12.277) is not necessary conservative due to the finite (i.e., in practice rather coarse) size grid resolution, and some sort of numerical trick must be used to enforce the conservative properties. It is mainly at this point in the formulation of the numerical algorithm that the class method of Hounslow et al [74], the discrete method of Ramkrishna [151] and the multi-group approach used by Carrica et al [24], among others, differs to some extent as discussed earlier. 

Although generalized coordinates are often useful for understanding molecular systems, as we shall learn in later chapters, it is desirable to avoid the use of formulations with a configmation-dependent mass matrix since it complicates the design of numerical methods with good conservation properties. This is one reason that molecular simulation methods are typically described in a Cartesian coordinate framework. 

Let us emphasize that the issues arising in the design and analysis of numerical methods for molecular dynamics are slightly different than those confronted in other application areas. For one thing the systems involved are highly structured having conservation properties such as first integrals and Hamiltonian structure. We address the issues related to the inherent structure of the molecular N-body problem in both this and the next chapter wherein we shall learn that symplectic discretizations are typically the most appropriate methods. 

The Verlet method is a numerical method that respects certain conservation principles associated to the continuous time ordinary differential equations, i.e. it is a geometric integrator. Maintaining these conservation properties is essential in molecular simulation as they play a key role in maintaining the physical environment. As a prelude to a more general discussion of this topic, we demonstrate here that it is possible to derive the Verlet method from the variational principle. This is not the case for every convergent numerical method. The Verlet method is thus a special type of numerical method that provides a discrete model for classical mechanics. 

The equations of DPD are a little more complicated than standard MD. It is therefore important to design accurate and efficient methods tailored to their special form. Most numerical methods for DPD simulation exhibit substantial statistical bias. Some methods have been developed to reduce the per timestep efficiency but at a price in terms of the conservation properties. The only way that the errors are tamed in some ad hoc schemes is by reducing the timestep size excessively (well below the stability threshold), but this may destroy the practical value of those methods. 

In the context of fluid flow, the Navier-Stokes Eqs. 14 and 17 are integrated numerically, together with the advectiOTi equation D/Dt = 0 for the level set function (i.e., Eq. 24 or 25 with a = 0 and b = 0), and the re-initializafimi Eq. 26 is applied periodically or at every time step to keep the level set function well behaved. While this is more straightforward than the procedures involved in the volume-of-fluid method, it has been found that the VOF method has better mass conservation properties than the simple level set method. 

One property of the exact trajectory for a conservative system is that the total energy is a constant of the motion. [12] Finite difference integrators provide approximate solutions to the equations of motion and for trajectories generated numerically the total energy is not strictly conserved. The exact trajectory will move on a constant energy surface in the 61V dimensional phase space of the system defined by. 

Numerous studies for the discharge coefficient have been pubHshed to account for the effect of Hquid properties (12), operating conditions (13), atomizer geometry (14), vortex flow pattern (15), and conservation of axial momentum (16). From one analysis (17), the foUowiag empirical equation appears to correlate weU with the actual data obtained for swid atomizers over a wide range of parameters, where the discharge coefficient is defined as — QKA (2g/ P/) typical values of range between 0.3 and 0.5. 

The major mechanism of a vapor cloud explosion, the feedback in the interaction of combustion, flow, and turbulence, can be readily found in this mathematical model. The combustion rate, which is primarily determined by the turbulence properties, is a source term in the conservation equation for the fuel-mass fraction. The attendant energy release results in a distribution of internal energy which is described by the equation for conservation of energy. This internal energy distribution is translated into a pressure field which drives the flow field through momentum equations. The flow field acts as source term in the turbulence model, which results in a turbulent-flow structure. Finally, the turbulence properties, together with the composition, determine the rate of combustion. This completes the circle, the feedback in the process of turbulent, premixed combustion in gas explosions. The set of equations has been solved with various numerical methods e.g., SIMPLE (Patankar 1980) SOLA-ICE (Cloutman et al. 1976). 

Chapter 8 presented the last of the computational approaches that I find widely useful in the numerical simulation of environmental properties. The routines of Chapter 8 can be applied to systems of several interacting species in a one-dimensional chain of identical reservoirs, whereas the routines of Chapter 7 are a somewhat more efficient approach to that chain of identical reservoirs that can be used when there is only one species to be considered. Chapter 7 also presented subroutines applicable to a generally useful but simple climate model, an energy balance climate model with seasonal change in temperature. Chapter 6 described the peculiar features of equations for changes in isotope ratios that arise because isotope ratios are ratios and not conserved quantities. Calculations of isotope ratios can be based directly on calculations of concentration, with essentially the same sources and sinks, provided that extra terms are included in the equations for rates of change of isotope ratios. These extra terms were derived in Chapter 6. 

Numerical experiments have shown that in many one dimensional systems with total momentum conservation, the heat conduction does not obey the Fourier law and the heat conductivity depends on the system size. For example, in the so-called FPU model, k IP, with (3 = 2/5, and if the transverse motion is introduced, / = 1/3. Moreover, in the billiard gas channels (with conserved total momentum), the value of P differs from model to model(Li and Wang, 2003). The question is whether one can relate / to the dynamical and statistical properties of the system. 

The excitation of oscillations with a quasi-natural system frequency and numerous discrete stationary amplitudes, depending only on the initial conditions (i.e. discretization of the processes of absorption by the system of energy, coming from the high-frequency source). A new in principle property is the possibility for excitation of oscillations with the system s natural frequency under the influence of an external high-frequency force on unperturbed linear and conservative linear and non-linear oscillating systems. 

---