Numerical solutions

In the numerical solution of the one-dimensional system, we used a three-point central differencing scheme such that the first and second spatial derivatives of C at point i (for a grid with non-uniform spacing) were 

In the i2-direction there are two finite bomidaries that must be considered, namely the R = 0 axis, and the singularity where the outer edge of the electrode meets the insulating smface, at i = 1. The spatial grid must therefore expand outwards from these two points, giving a high density of points at both edges of the electrode and a lower density in the middle of it [1]  

We apply the same three-point central differencing scheme as was used in the one-dimensional case to this two-dimensional space, which gives the form of the derivatives for the concentration at point i, j) as 

In Chapter 3, we chose to use an implicit method of solution, as opposed to an explicit one for reasons of stability and simulation efHciency (despite the greater complexity of the implicit method). The implicit discretisation 

In the ADI method [4-6], the explicit and implicit methods are combined. The timesteps, AT, are divided into two half-timesteps, ATf2. For the first of these half-timesteps, the derivatives along the. -coordinate 

A great many stagnation-flow solutions have been shown in the previous sections. The solution algorithms are similar to those discussed in Chapter 16. However, there are some differences and some complications that arise to deal with the special needs of stagnation flows. One issue has to do with computing the pressure-gradient eigenvalue in finite-gap problems. Another has to do with velocity reversal in opposed flows (Section 6.10). 

Consider first how the eigenvalue can be determined, which is mainly concerned with the continuity and momentum equations. The continuity equation 

This difference formula propagates axial-velocity u information from the lower boundary (e.g., stagnation surface) up toward the inlet-manifold boundary. At the stagnation surface, a boundary value of the axial velocity is known, u = 0. A dilemma occurs at the upper boundary, however. At the upper boundary, Eq. 6.106 can be evaluated to determine a value for the inlet velocity u j. However, in the finite-gap problem, the inlet velocity is specified as a boundary condition. In general, the velocity evaluated from the discrete continuity equation is not equal to the known boundary condition, which is a temporary contradiction. The dilemma must be resolved through the eigenvalue, which is coupled to the continuity equation through the V velocity and the radial-momentum equation. 

It is easy to see that changing the value of Ar produces a different V profile through the solution of the momentum equation. Then, in turn, different V profiles produce different u profiles through the solution of the continuity equation. It is possible to find a value of Ar such that the continuity equation is satisfied at the top boundary and the known axial-velocity boundary conditions is satisfied, that is, uj = f/iniet However, it is inefficient to carry out such an iteration explicitly. Rather, it it is more efficient to implement an implicit boundary condition. 

To maintain a banded Jacobian structure, the eigenvalue is represented through a trivial differential equation, 

and (2.31) with Eq. (2.9) can be solved simultaneously without simplification. Since the analytical solution of the preceding simultaneous differential equations are not possible, we need to solve them numerically by using a computer. Among many software packages that solve simultaneous differential equations, Advanced Continuous Simulation Language (ACSL, 1975) is very powerful and easy to use. 

The heart of ACSL is the integration operator, INTEG, that is, 

Original set of differential equations are converted to a set of first-order equations, and solved directly by integrating. For example, Eq. (2.14) can be solved by integrating as 

For more details of this simulation language, please refer to the ACSL User Guide (ACSL, 1975). 

You can also use Mathematica (Wolfram Research, Inc., Champaign, IL) or MathCad (MathSoft, Inc., Cambridge, MA), to solve the above problem, though they are not as powerful as ACSL. 

The mathematical model developed in the preceding section consists of six coupled, three-dimensional, nonlinear partial differential equations along with nonlinear algebraic boundary conditions, which must be solved to obtain the temperature profiles in the gas, catalyst, and thermal well the concentration profiles and the velocity profile. Numerical solution of these equations is required. 

A comparison of the benefits and drawbacks of common numerical solution techniques for complex, nonlinear partial differential equation models is given in Table II. Note that it is common and in some cases necessary to use a combination of the techniques in the different dimensions of the model. 

Numerical Solution Techniques for Partial Differential Equations Arising in Packed Bed Reactor Modeling 

Finite difference Has simple construction Is easily extended to multiple dimensions Solution is stable for sharp gradients or in response to a concentration or temperature front Leads directly to state-space representation Is often computationally prohibitive since accurate solutions may require a large number of grid points 

Methods of Generally needs few grid Can have difficulties with 

We employ the Newton-Raphson method to iterate toward a set of values for the unknown variables (nw, mi, mp)r for which the residual functions become vanishingly small. 

To do so, we calculate the Jacobian matrix, which is composed of the partial derivatives of the residual functions with respect to the unknown variables. Differentiating the mass action equations for aqueous species Aj (Eqn. 4.2), we note that, 

For the K and Freundlich models, as mentioned, there is no basis entry Ap and hence we do not write a residual function of the form Equation 9.42, nor do we carry Jacobian entries for Equations 9.46 or 9.49-9.52. 

To evaluate the Jacobian matrix, we need to compute values for dmq/dnw, dmq/dnii, and dmq/dmp. For the Kt and Freundlich models, 

No value for dmq/dmp is needed to evaluate the and Freundlich models. For the Langmuir model and the ion exchange model under the Gaines-Thomas and Gapon conventions, 

The procedure for solving the governing equations parallels the technique described in Chapter 5, with the added complication of accounting for electrostatic effects. We begin as before by identifying the nonlinear portion of the problem to form the reduced basis 

For each basis entry we cast a residual function, which is the difference between the right and left sides of the mass balance equations (Eqns. 8.11a, b, and e) 

At each step in the iteration, we evaluate the residual functions and Jacobian matrix. We then calculate a correction vector as the solution to the matrix equation 

Quasi-single models typically set the axial velocity at the inlet to the gas punch velocity thus, 

Two-phase models typically assume that only the gas enters the bath at the nozzle such that, 

Assuming a 10-pct turbulence intensity generated by the incoming gas bubbles, Schwarz and Turner [40] and Ilegbusi et al. [8] used the following boundary condition for the k and e at the nozzle  

Integration of the governing partial differential equations presented in Sect. 9.3 over a control volume leads to an algebraic equation of the form  

In multi-fluid models, there is a separate solution field for each phase. Transported quantities interact via interphase terms. The pressure in both phases is assumed to be the same within a computational cell. Field equations for each phase are weighted with the volume fraction of that phase. The model is solved using the inter-phase slip algorithm (IPSA) embodied in the PHOENICS computational code with modifications for the interfacial parameters and other source terms. This approach and its derivatives have been adopted in other commercial computational codes. 

Process engineering and design using Visual Basic 

The theoretically calculated values and values obtained through numerical analysis are presented in Table 1.11. 

---

No a priori information about the unknown profile is used in this algorithm, and the initial profile to start the iterative process is chosen as (z) = 1. Moreover, the solution of the forward problem at each iteration can be obtained with the use of the scattering matrices concept [8] instead of a numerical solution of the Riccati equation (4). This allows to perform reconstruction in a few seconds of a microcomputer time. The whole algorithm can be summarized as follows ... 

The deconvolution is the numerical solution of this convolution integral. The theory of the inverse problem that we exposed in the previous paragraph shows an idealistic character because it doesn t integrate the frequency restrictions introduced by the electro-acoustic set-up and the mechanical system. To attenuate the effect of filtering, we must deconvolve the emitted signal and received signal. 

Lane improved on these tables with accurate polynomial fits to numerical solutions of Eq. 11-17 [16]. Two equations result the first is applicable when rja 2... 

We will describe integral equation approximations for the two-particle correlation fiinctions. There is no single approximation that is equally good for all interatomic potentials in the 3D world, but the solutions for a few important models can be obtained analytically. These include the Percus-Yevick (PY) approximation [27, 28] for hard spheres and the mean spherical (MS) approximation for charged hard spheres, for hard spheres with point dipoles and for atoms interacting with a Yukawa potential. Numerical solutions for other approximations, such as the hypemetted chain (EfNC) approximation for charged systems, are readily obtained by fast Fourier transfonn methods... 

The solutions to this approximation are obtained numerically. Fast Fourier transfonn methods and a refomuilation of the FINC (and other integral equation approximations) in tenns of the screened Coulomb potential by Allnatt [M are especially useful in the numerical solution. Figure A2.3.12 compares the osmotic coefficient of a 1-1 RPM electrolyte at 25°C with each of the available Monte Carlo calculations of Card and Valleau [ ]. 

Using the Omstein-Zemicke equation, numerical solutions for the restricted primitive model can be. 

Shampine S 1994 Numerical Solutions of Ordinary Differential Equations (New York Chapman and Hall)... 

Numerical solution of this set of close-coupled equations is feasible only for a limited number of close target states. For each N, several sets of independent solutions F.. of the resulting close-coupled equations are detennined subject to F.. = 0 at r = 0 and to the reactance A-matrix asymptotic boundary conditions,... 

Wlien describing the interactions between two charged flat plates in an electrolyte solution, equation (C2.6.6) cannot be solved analytically, so in the general case a numerical solution will have to be used. Several equations are available, however, to describe the behaviour in a number of limiting cases (see [41] for a detailed discussion). Here we present two limiting cases for the interactions between two charged spheres, surrounded by their counterions and added electrolyte, which will be referred to in further sections. This pair interaction is always repulsive in the theory discussed here. 

Theorem 1 ([8]). Let H be analytic. There exists some r > 0, so that for all T < Tt the numerical solution Xk = ) Xo and the exact solution x of the perturbed system H (the sum being truncated after N = 0 1/t) terms) with x(0) = Xq remain exponentially close in the sense that... 

The typical dependence of the stable stationary solutions to (4) on the control parameter of the model Xe is presented in Fig. 1. These results have been obtained as numerical solutions of (4) with equal to... 

D. Beglov and B. Roux. Numerical solutions of the hypernetted chain equation for a solute of arbitrary geometry in three dimensions. J. Chem. Phys., 103 360-364, 1995. 

We further discuss how quantities typically measured in the experiment (such as a rate constant) can be computed with the new formalism. The computations are based on stochastic path integral formulation [6]. Two different sources for stochasticity are considered. The first (A) is randomness that is part of the mathematical modeling and is built into the differential equations of motion (e.g. the Langevin equation, or Brownian dynamics). The second (B) is the uncertainty in the approximate numerical solution of the exact equations of motion. 

Consider a numerical solution of the Newton s differential equation with a finite time step - At. In principle, since the Newton s equations of motion are deterministic the conditional probability should be a delta function... 

We have in mind trajectory calculations in which the time step At is large and therefore the computed trajectory is unlikely to be the exact solution. Let Xnum. t) be the numerical solution as opposed to the true solution Xexact t)- A plausible estimate of the errors in X um t) can be obtained by plugging it back into the differential equation. 

The last approximation is for finite At. When the equations of motions are solved exactly, the model provides the correct answer (cr = 0). When the time step is sufficiently large we argue below that equation (10) is still reasonable. The essential assumption is for the intermediate range of time steps for which the errors may maintain correlation. We do not consider instabilities of the numerical solution which are easy to detect, and in which the errors are clearly correlated even for large separation in time. Calculation of the correlation of the errors (as defined in equation (9)) can further test the assumption of no correlation of Q t)Q t )). 

L.R. Petzold, L.O. Jay, and J. Yen. Numerical solution of highly oscillatory ordinary differential equations. Acta Numerica, pages 437-484, 1997. 

Computational issues that are pertinent in MD simulations are time complexity of the force calculations and the accuracy of the particle trajectories including other necessary quantitative measures. These two issues overwhelm computational scientists in several ways. MD simulations are done for long time periods and since numerical integration techniques involve discretization errors and stability restrictions which when not put in check, may corrupt the numerical solutions in such a way that they do not have any meaning and therefore, no useful inferences can be drawn from them. Different strategies such as globally stable numerical integrators and multiple time steps implementations have been used in this respect (see [27, 31]). 

The weighted residual method provides a flexible mathematical framework for the construction of a variety of numerical solution schemes for the differential equations arising in engineering problems. In particular, as is shown in the followmg section, its application in conjunction with the finite element discretizations yields powerful solution algorithms for field problems. To outline this technique we consider a steady-state boundary value problem represented by the following mathematical model... 

The comparison between the finite element and analytical solutions for a relatively small value of a - 1 is shown in Figure 2.25. As can be seen the standard Galerkin method has yielded an accurate and stable solution for the differential Equation (2.80). The accuracy of this solution is expected to improve even further with mesh refinement. As Figmre 2.26 shows using a = 10 a stable result can still be obtained, however using the present mesh of 10 elements, for larger values of this coefficient the numerical solution produced by the standard... 

The finite element results obtained for various values of (3 are compared with the analytical solution in Figure 2.27. As can be seen using a value of /3 = 0.5 a stable numerical solution is obtained. However, this solution is over-damped and inaccurate. Therefore the main problem is to find a value of upwinding parameter that eliminates oscillations without generating over-damped results. To illustrate this concept let us consider the following convection-diffusion equation... 

Johnson, C., 1987. Numerical Solution of Partial Deferential Equations by the Finite Element Method, Cambridge University Press, Cambridge. 

Lapidus, L. and Pinder, G. F., 1982. Numerical Solution of Partial Differential Equations in Science and Engineering, Wiley, New York. 

Taylor, C. and Hood, P., 1973. A numerical solution of the Navier-Stokes equations using the finite element technique. Comput. Fluids 1, 73-100. 

Further details of the BB, sometimes referred to as Ladyzhenskaya-Babuska-Brezi (LBB) condition and its importance in the numerical solution of incompressible flow equations can be found in textbooks dealing with the theoretical aspects of the finite element method (e.g. see Reddy, 1986), In practice, the instability (or checker-boarding) of pressure in the U-V-P method can be avoided using a variety of strategies. 

Morton, K. W., 1996. Numerical Solution of Convection Diffusion Problems, Chapman Hall, London. 

---