## 741 Complex Geometry

Many flows relevant to chemical reactor engineering practice involve complex geometries. Although the principles of discretization and solution methods for algebraic systems described earlier may be used, some modifications are required to handle such complex geometries. The properties of the solution algorithm depend on choice of the grid and the arrangement of variables on the grid. Some of these issues are discussed in this section.

For a regular (for example, rectangular or circular) geometry, the grid lines usually follow the co-ordinate directions. In complicated geometries, the choice of grid is not trivial. The grid is subject to constraints imposed by the discretization method and solution algorithm. If the algorithm is designed for structured quadrilaterals, an unstructured grid consisting of triangles cannot be used. When the geometry is complex, some compromises have to be made to fulfill the constraints. Body-fitted non-orthogonal grids are most often used to calculate flows in complex geometries (most commercial codes use such grids). In such grids, grid lines follow the boundaries of the solution domain, which makes implementation of boundary conditions easier. The transformed equations for non-orthogonal grids, however, contain additional terms leading to difficulty in programming and increased computational costs. Despite this disadvantage, these grids are used in most applications. Grid generation for complex geometries will not be discussed here. Some relevant comments are included in Chapter 1. More details of grid generation can be found in Thompson et al. (1985) and Arcilla et al. (1991). Some general comments relevant to numerical solutions are included here.

Though complex geometry demands that the grid be non-orthogonal, it is useful to make it as orthogonal as possible. In finite volume methods, orthogonality of grid lines at corners (vertices) of computational cells (CV) is not important. The angle between the cell face normal vector and the line connecting the CV centers on either side is important. Cell topology is also important. If the midpoint rule integral approximation and linear interpolation is used, then the accuracy will be higher if CVs are quadrilaterals than if CVs are triangles (in 2D). Accuracy is also improved if one set of grid lines closely follows the streamlines of flow (especially for convective flows). This can be achieved with a structured grid (of quadrilaterals in 2D) but not with triangular grids. Complex geometries also demand non-uniform grids. A finer grid should be used in the regions where strong variations are expected to occur. When generating non-orthogonal and non-uniform grids, three measures of grid quality should be kept in mind:

• ratio of adjacent cell sizes: preferably less than two;

• aspect ratio of a computational cell (ratio of adjacent edges of a cell): preferably less than five;

• skewness of a computational cell (angle between adjacent edges of a cell): preferably greater than 45°.

When the geometry is complex, grid generation may require a significant fraction of the time necessary to complete the development and application of the computational model. Since the accuracy of the flow solution depends as much on the grid quality as on the approximations used for discretization of the equations, time spent generating a quality grid is a worthwhile investment.The basic principles of the finite volume method discussed in the previous chapter are independent of the type of grid used. There are, however, some new features, which need to be introduced to handle arbitrary non-orthogonal grids. Some of these are discussed here.

As discussed in the previous chapter, with the finite volume method, we need to approximate the surface and volume integrals to calculate fluxes and sources. Consider the two-dimensional control volume shown in Fig. 7.21. Following Eq. (6.6), the mass FIGURE 7.21 A typical control volume (2D) for non-orthogonal grid.

flux through face e can be written as:

where n is a unit normal vector to the east face. Unlike with a regular grid, the surface vector has components in more than one Cartesian direction and all the velocity components contribute to the mass flux. Mass flux can be calculated by summing products of each Cartesian velocity component and the corresponding surface vector component (projection of the cell face on a Cartesian co-ordinate plane):

where superscripts or subscripts indicate Cartesian co-ordinate direction. The con-vective flux of any transported quantity can be calculated as the product of mass flux through the face and the value of transported property at the center of the cell face. Various interpolations discussed in the previous chapter can be used to estimate the value of transported property at the center of the cell face. The integrated diffusive flux in a general coordinate system can be written:

The gradient of <p at the cell face center can be expressed either in terms of the derivatives with respect to global Cartesian co-ordinates or local orthogonal co-ordinates. If the local orthogonal system attached to the cell face center is used, then only the derivative in the normal (n) direction contributes to the diffusive flux:

On a Cartesian grid, n = x at the e face and the usual schemes can be used to estimate gradient at e. For example, a central differencing scheme will give:

d n Je Lpe where LPE is the distance between nodes E and P. When the grid is non-orthogonal, the above expression must be corrected by a term containing the difference between the gradients in the f and n direction. Mujaferija (1994) has proposed the correction as:

where if is the unit vector in the f direction (Fig. 7.21). The second term of the right-hand side is usually evaluated explicitly from previously known values.

Apart from the convective and diffusive fluxes, it is also necessary to evaluate source terms. As mentioned in the previous chapter, the volume integral can be calculated as a product of the CV center value of the integrand and the CV volume. This approximation is independent of the CV shape. For non-orthogonal grids, the calculation of the cell volume becomes more complicated. For 2D quadrilaterals, the volume of the cell can be calculated by taking the vector product of the two diagonals. The expression for cell volume of the cell shown in Fig. 7.21 becomes:

^cell = ^ [(Xne - ^sw^nw - Jse) - few - *seX^ne - Jsw)] (7.24)

In three-dimensional solution domains, the cell faces are not necessarily planar. Suitable approximations are necessary. Further details of discretization on non-orthogonal grids will not be presented here. The purpose of this section was just to highlight some of the relevant issues. Reader may refer to Ferziger and Peric (1995) and references cited therein for more details. Even if a reactor engineer is not interested in developing an in-house CFD code, familiarity with these issues is essential for developing computational models of complex flows in complex geometry. More often than not, poor quality grids lead to divergence, and suitable corrective measures need to be taken to obtain convergence. Some examples of the simulation of complex flows in industrial equipment are discussed in Chapters 9 to 14.

### 7.4.2. Enhancing Convergence Performance

Overall convergence performance depends on several factors. Factors such as, large number of computational cells, overly conservative under-relaxation factors, strong non-linearity and coupling between different equations, hinder convergence. Grid quality also affects the convergence rate significantly. It is often necessary to use various 'tricks' to enhance the convergence performance. Different classes of flow problems (with varying degrees of complex flow physics and complex geometry) may require different strategies. Some of the general ways of enhancing convergence are discussed here.

Supplying an initial guess for important flow variables often enhances overall convergence behavior. The process of providing an initial guess also allows the user to examine the main characteristic space and time scales. Examination of these scales is necessary to select appropriate numerical parameters such as time step and under-relaxation factors. The process is also useful to verify the adequacy of the generated grid. Another commonly used technique to enhance convergence performance of complex flow problems is to break down the overall problem into a sequence of problems with increasing complexity at each stage of the sequence. For example, when solving a non-isothermal problem, it may be a good idea to first obtain a reasonable solution (need not be completely converged) to the isothermal problem, which can act as a good initial guess for the non-isothermal problem. Solution of the simplified isothermal problem can be started by setting the initial guess to the temperature field, which will remain unchanged during the solution process. After obtaining reasonable convergence for the isothermal problem, solution of the enthalpy equation can be started. Sometimes it is advantageous (computationally) to switch off the solution of momentum and continuity equations while solving the enthalpy equations. Once the enthalpy equation starts converging, all the equations can be included in the solution process. Several examples of such a step-by-step procedure to improve the convergence of complex flow problems are discussed in Chapter 9 and later chapters.

As mentioned in the previous chapter, it is advantageous to use under-relaxation factors to improve the convergence behavior of a set of non-linear coupled equations. Under-relaxation needs to be applied not only to direct flow variables such as velocity, pressure, temperature etc. but also to indirect variables like fluid properties (density, viscosity, heat capacity), when these variables are composition or temperature dependent. It is indeed impossible to provide generalized guidelines for setting under-relaxation parameters since several factors such as the algorithm, parameters of the linear equation solver, the extent of non-linearity and coupling affect suitable values of under-relaxation factors. Usually, the solution is started with rather conservative values of under-relaxation factors. The values may be increased as the solution progresses with the help of continuous monitoring of residual history. The influence of under-relaxation factors on convergence rate is demonstrated in Fig. 6.18. Unless repeated simulations of similar cases have to be carried out, it may not be computationally economical to optimize under-relaxation factors.

Apart from the formulation of discretized governing equations and algorithms to treat various couplings, the solution method of the resulting linear algebraic equations also has a significant impact on convergence rate. In many CFD simulations of complex flows, iterative line-by-line solvers are used. Two parameters of line-by-line solvers, namely, the direction of sweeping lines in the solution domain and the number of sweeps for each equation, govern the overall convergence behavior. In general, lines which are normal to the primary flow direction, are solved in the direction of the primary flow. When there is no single dominant flow direction, it is useful to use alternating sweep directions. An increase in the number of sweeps for any equation leads to more computations per iteration; however, this may improve the local and therefore, global convergence behavior. Of course, there will be an optimum number of internal iterations or of sweeps throughout the solution domain, to minimize overall convergence time. Generally, a higher number of sweeps needs to be specified for pressure and for any equation which is difficult to converge (e.g. species equations in the presence of non-linear chemical reactions). When the flow problem involves large body forces, it is often necessary to increase the number of sweeps of the pressure equation.

In addition to these parameters of line-by-line solvers, several other techniques have been proposed to accelerate the convergence rate. Some of these methods are discussed in Section 6.2.2 and only brief comments will be added here. One-dimensional block corrections based on an additive correction philosophy (Kelkar and Patankar, 1989) reduce long wavelength errors in the direction in which they are applied. However, these methods are not suitable if very steep gradients exist in the solution. Multigrid methods are effective in reducing long wavelength errors in such situations. As mentioned earlier, choice of such parameters as number of grid levels, number of iterations on each level, the order in which various levels are visited, interpolations between various levels and so on will affect convergence behavior. The values of these parameters need to be tuned to accelerate the overall convergence performance. It is not possible to discuss all these details here, and the reader is referred to Hackbusch (1985) and Ferziger and Peric (1995).

For time-dependent flows, in addition to the factors discussed above, the values of time step and number of iterations per time step, govern the overall convergence behavior. In general, the selected time step should be at least an order of magnitude smaller than the smallest relevant time scale of the modeled flow process. If the explicit method is used to march in time, iterations within a time step are not required. However, in such a case, the value of time step is severely constrained. For an implicit method, it may be necessary to use more than one iteration per time step to obtain a converged solution at each time step. In general, the value of time step is selected in such a way that the number of iterations per time step does not increase beyond 20. For many practical transient flows, a very fast 'start-up' phase exists, which may decay rapidly. In such cases, the value of time step may be gradually increased as the calculations proceed. When the problem is being solved as a transient problem and the steady or pseudo-steady solution is of interest, it may not be necessary to obtain complete convergence for each time step at the beginning. Once the solution is developed, the internal iterations can be increased to achieve the desired accuracy of the solution. This simple technique is useful for solving complex multiphase flows, which are often solved as time-dependent flows. In addition, several problem-dependent techniques can be used to accelerate the overall convergence behavior of multiphase flows. Examples of these are discussed in Chapters 9 to 13.

7.4.3. Error Analysis of Complex Simulations

Errors in CFD simulations arise mainly from two sources:

• inherent errors in representing reality by the set of model equations (physically deficient representation); and

• errors arising from inexact solution methods (numerically deficient representation).

It is essential to identify and separate these two types of errors to avoid confusion. If numerical errors are not isolated, they may lead to undesirable spurious model calibration exercises. It is, therefore, necessary to devise systematic methods to quantify numerical errors. The basic idea behind error analysis is to obtain a quantitative measure of numerical errors, to devise corrective measures to ensure that numerical errors are within tolerable limits and the results obtained are almost independent of numerical parameters. Having established adequate control of numerical errors, the simulated results may be compared with experimental data to evaluate errors in physical modeling. The latter process is called model validation. Several examples of model validation are discussed in Chapters 10 to 14. In this section, some comments on error analysis are made.

Low-order numerical methods contribute to the robustness and computational efficiency of the CFD code. However, this same robustness and speed make it imperative that error estimates be available so that plausible looking results are not confused with accurate results. In many complex reactor-engineering applications, the reactor engineer may be interested in capturing the trends rather than absolute values. Even then, it is essential to verify that captured key flow features are independent of numerical parameters. Formal error estimates for grid-based numerical methods may be based on Taylor series expansions (Roache, 1976). These analytic approaches are valuable for development and evaluation of numerical methods but are rarely used to assess complex flow simulations. The usual method of assessing the numerical accuracy of complex flows is through grid refinement. In this approach, the computations are repeated on progressively finer grids and resolution is presumed to be adequate when results do not change significantly with grid density. Estimates of the local converged solution may be obtained by extrapolating solutions obtained on two or more grids to an infinitely dense grid (zero grid spacing). An example of such a method is discussed in Chapter 6. It must, however, be remembered that routine application of grid refinement to complex flow simulations (with complex flow physics and complex geometry) is problematic for several reasons. Most flow simulations relevant to industrial reactor engineering routinely use a few hundred thousand computational cells (see examples in Chapter 9 to 14). Moreover, for complex geometric configurations, numerical errors resulting from non-uniform distribution of grids, and from departures from grid orthogonality further complicate the quantification of numerical errors. In such cases, systematic assessment of grid sensitivity requires variation in both grid density as well as grid distribution and topology. More often than not, difficulties in grid generation and the magnitude of computational resources required for such complex flow simulations preclude the routine use of grid sensitivity tests.

When systematic grid refinement is not possible to assess numerical errors, global balances of numerically non-conserved quantities (quantities that are not conserved at the control volume level in the construction of the numerical scheme) can be used. Haworth et al. (1993) described this approach in detail. The method is based on volume-integrated partial differential equations for primary or derived physical quantities of interest. Balances can be applied to the full computational domain or to any sub-domain down to the single-cell level. Comparison of relative magnitudes of terms in the balances provides insight into the physics of the flow being computed. For quantities that are not conserved on a computational cell level, the imbalance provides a direct measure of numerical inaccuracy using a single grid simulation. Haworth et al. (1993) recommended imbalance in mean kinetic energy, which they show to be a good measure of low-order spatial discretization error. For more details, the original paper may be consulted. Such an approach, combined with grid refinement studies (wherever possible), can provide a useful measure of numerical errors.

For multiphase flows, the analysis of numerical errors becomes more difficult but even more important. In many practical multiphase flows, it may not be possible to obtain grid-independent solutions. More often than not, the reactor engineer has to rely on CFD simulations, which may show grid dependence, when making important practical design decisions. In such cases, it is essential to ensure that grid dependence is not affecting key conclusions on which the engineering decision is being based (even though the overall flow field may show grid dependence). Some examples of developing useful engineering decisions based on flow simulations obtained with modest grid requirements are discussed in Chapters 9 to 14.

0 0