Numerical methods for solving the reactive diffusion equation in complex geometries

Introduction

The reactive diffusion equation is a system of (in general nonlinear) partial differential equations and describes the time evolution of mixable several chemical species which can also react to new chemical species. For  distinct chemical species it has the form:

While  denotes the index of the chemical species,  is their concentration,  is their diffusion flux vector and  is their gain and loss due to chemical reactions. Since the diffusion flux and the source term  is a nonlinear function in terms of the concentrations of all other chemical species, equation (1) cannot be solved analytically. This equation must be discretized  by a suitable discretization method such that a difference equation can be obtained. The solution of this discretized equation is less accurate than the exact solution, but the more discrete elements (e.g. lattices; denoted by )  are used for discretization the more accurate is the numerical solution. Moreover, the discrepancy between exact and numerical solution depends on the geometry of the discrete elements that are chosen to discretize the region in which the reactive diffusion equations are solved (denoted by ). Especially in complex geometries (e.g. porous media) it becomes difficult to fit the geometry of the discrete elements to the geometry of the solution region . In this cases there must be used more complicated models (Tompson 1992) than usual discretization methods. Typical methods for solving the reactive diffusion equation is the finite element method (FEM) (John 2008). Another approach for the solution of reactive diffusion equations are used in (Liu 2005).

Description of the numerical method

In general, the numerical solution of (1) requires that  is decomposed into  discrete elements as:

The decomposition (2) is equivalent to the ansatz

for some coefficients  and the ansatz functions  that can also be assumed species dependent. When giving some boundary and initial conditions, substituting (3) into (1) and integrate over the whole region  and the time interval , the discrete reactive diffusion equation can be solved by iteration such that the set of coefficients  can be obtained. In this paper it is illustrated, how complex geometries can be discretized by decomposing  into elements of a mathematical graph. Mathematical graphs consist on vertices , edges , but also areas  and volumes  that are enclosed by the edges (or areas) of the graph. For simplicity, the time discretization is performed by setting  for some time scale ; to make the ansatz functions time independent, it is obvious to make the coefficients  dependent on time. Another simplification is the assumption of species independent ansatz functions  (species-dependent information is given by the coefficients ). To characterize the ansatz functions as graph elements, the decomposition (3) has the form:

These ansatz functions  have the value 1 when  lies on the corresponding discrete element and 0 otherwise.

When (4) is substituted into (1), every ansatz function has to be derived by the coordinates. Coordinate derivatives can be performed with the aid of the Gaussian law , where  is the unit normal vector corresponding to the surface element  and  is an infinitesimal 3-dimensional region around the discrete element . For example,  corresponds to an infinitesimal 3-dimensional ball around the vertex . The volume integral can be approximated to  (an analogous approximation holds for the surface integral). To obtain the discrete reactive diffusion equation, products of the ansatz functions must be integrated over the entire region . Per definition, the product of ansatz functions satisfies the simple relation: . Finally, one obtains the system of difference equations:

Equation (5) is satisfied if and only if the coefficients  vanish.

Conclusions

When a numerical simulation is applied to geometries with subregions that can be approximated as vertices (e.g. very small cavities), edges (e.g. tube-like structures), etc., the numerical method described above canbe applied. Discretizing these geometries by graphs can be useful, since the geometry has a structure that is difficult to treat only with 3-dimensional volume elements like in FEM. If a small subregion is subdivided into many 3-dimensional discrete elements, then more computation time than in the case of only one single discrete element is necessary. 

References

[1]Tompson, A. et al. Particle-grid methods for reacting flows in porous media with application to Fisher's equation, Applied Mathematical Modelling, vol. 16, i. 7 (1992), pp. 374-383.

[2]John,V. et al. Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion, Comput. Methods Appl. Mech. Engrg. 198 (2008) 475–494.

[3] Liu,J. et al. An Operator Splitting Method for Nonlinear Reactive Transport Equations and Its Implementation Based on DLL and COM, Current Trends in High Performance Computing and Its Applications (2005), pp. 93-102