FORMULATION OF A REACTIVE-DIFFUSIVE-CONVECTIVE MODEL SOLVED BY THE FINITE ELEMENT METHOD COUPLED BY SYMMETRIZED STRANG SPLIT

Abstract. This work develops a mass transport model coupled with a reactive kinetics model through the Finite element method, to simulate the behavior of pollutants substances in an urban atmosphere. In this work a Symmetrized Strang Split method is used to couple the mass transport model with the reactive kinetics model. This split method solves separately the diffusion, advection and reaction effects in a cascade process. The formulation of the mass transport process is showed and its results are compared with a fully coupled finite element method. This comparison is important to show the solution time and computational cost and validity of the split model.


INTRODUCTION
One of the most common physical phenomena that occur is the mass transport.This phenomenon is responsible for the propagation of particles that could be pollutant or inert in a liquid or gaseous fluid.This process can occur in atmospheric, river, sea and groundwater contamination, in water purification processes, in rivers' erosion, in transport of substances into blood system, in bacteria's growth, in the fabrication processes of pharmaceutics [1], in the obtaining of metals and when the achieving of a product with a desired physicals and chemicals characteristics is important [2], among others.The numerical simulation of mass transport is widely studied due to the need of its better understanding and the properly optimization of each one of the processes mentioned above.
The reaction is one of the processes include in the physical phenomena mentioned above.Even this process isn't very complex to represent mathematically, it causes that the equations system that represent the mass transport phenomena becomes coupled and in some cases nonlinear [3] [4].One proper alternative to solve a coupled equations system is to use a split method [5], it implies that the parts of equation which symbolize the coupled processes could be decoupled, so it could be solved separately from the others processes, substantially reducing the computational cost and the simulation time in relation to the coupled solution.
The atmospheric contamination is one of the physical phenomenon of mass transport that have a great importance for the research community, due to the harmful impact over the public health, the flora and fauna carried for the atmospheric pollutants [6], [7].The use of a numerical model simplifies the studies over pollutant substances behavior in urban atmos-phere [8], [9], [10], because, in relation with the field studies, the numerical model obtains the results faster, with lower cost and makes possible the development of preventive models.This kind of phenomenon includes not only the advection and diffusion processes but also reaction and deposition in some cases.
This paper develops a mathematically representation of transport phenomena for the reactive species NO, NO2 and O3, which are part of the pollutant present in an urban center and that are related to the emission gases of automobiles.The paper first presents the formulation of the physical phenomena and the considerations associated with this, and then explains the fundaments of the finite element method, the SUPG technique and the Symmetrized Strang Split methodology.At last presents the factorial experiment design used to analyze and compare the results of the split method with those obtained with a coupled method.

FORMULATION
The implementation of the mass transport phenomenon through the finite element method applying the SUPG technique is presented below.The system of differential equations that represents this phenomenon is expressed in (1).
With initial and boundary condition defined as: ( Where is the concentration value for i chemical species, ̅ is the velocity field associated with the advective process, is the molecular diffusion coefficient for each substances, is the turbulent diffusion coefficient and is the reactive term that is defined for each substance depending of the reaction mechanism, define the value of initial condition, is the function that define the value of scalar field and is the function that define the value of flow in the boundary.Using the method of finite elements is obtained a system of weak discretized spatial equations as: Where ̂ is the discrete approximation of concentration defined by: is the shape function, is the elemental domain, is the elemental boundary and is the weight function.The weight function is modified using the SUPG technique for improve the stabilization of the system (5).(5) In this last equation h is the characteristic element size and α is a perturbation parameter calculated with (6) [11].
Where Pe is the dimensionless Peclet number defined by | ̅| .For simplicity the system of weak discretized spatial equations is presented in the matrix form (7). Where.
∫ 12) is the capacitive matrix, is the advective Matrix, is the diffusion trix, is the reactive matrix and is the load vector.Besides is necessary the implementation of a temporal discretization scheme (13), in this case is used the implicit scheme because is unconditionally stable [3], [12].
Where .Next step is solving the equations system considering it as a system of linear equations.

OPERATOR SPLIT METHODOLOGY
When the reaction mechanisms are coupled in the finite element discretization, the resulting system becomes stiff and occasionally nonlinear.This causes that the solution process requires smaller time steps, more iterations at each time step and more computational cost [13].
Operator splitting methods are a useful tool to solve this kind of problems without the disadvantages of the couple methodology.This technique was originally developed for simplify the multidimensional problem into several one-dimensional problems and treat them individually using specialized numerical algorithms making the solution easier.The implementation of the operator splitting methodology divides a system of advection-diffusionreaction equations into a number of sub problems, which are solved independently on every time step.With the implementation of this technique is possible to use different numerical discretization for every single part of the partial differential equations.
These kinds of methods are also called time splitting methods or fractional methods, because the split divide the equations in parts that are solved in different time steps.There are several kinds of split methods like the Yanenko and the Symmetrized Strang [14].
The Yanenko method, an operator splitting method, resolves the problem dividing it in two sub problems that are solved in cascade.Below present the split configuration: (14) Where and are the subproblems that represent the terms to the differential equations.The Yanenko method has the procedure presented in the figure 1 to solve the two operators.In this method is possible to use different discretization and apply explicit or implicit schemes for each sub problem.The Yanenko method first order accuracy, only one step for each sub problem, causes operator-splitting errors resulting from both the operator-splitting approach itself and from the time and spatial discretization, which is clearly a weakness.
The second method is the Symmetrized Strang Split, this method has the same characteristics of the last one, however Strang splitting has advantages over Yanenko method because this method is second order accurate.The split configuration is presented in (14), and its procedure is presented below in figure 2. and are the subproblems that represent the terms to the differential equations.In some cases the number of terms in the equations are more than two, therefore it can be grouped in different ways.For example if there are three terms then is possible to use any of the next setups.
The two methods are useful to implement and solve models of transport phenomena, but considering that the symmetrized strang split has a best accuracy, it was selected.

NUMERICAL EXPERIMENTATION (TEST PROBLEM)
In the implementation of the operator splitting method the numerical representation is divided in two parts, the first part represent the advection diffusion processes and is solved with (16), the second part, the reaction process, is represented with (17) and is solved with any ordinary differential equation solver.

First Part
A multilevel factorial design with the test problem presented in figure 3 is implemented to compare the coupled solution with the symmetrized strang split method.In the test problem is simulated the transport of the reactive species NO, NO2 and O3 which are part of the pollutant species presented in an urban center and that are related with the gases emission of automobiles.The domain for this problem has a rectangular shape with dimensions of 20 x 20 meters and the conditions for all the frontiers are presented in (18).
(18) Also is included a velocity field of a rotational flow as showing in figure 3, this field is based in a Gaussian bell presented in figure 4 and expressed numerically in (19).Additional to the field velocity is included a homogenous diffusion coefficient in the entire domain to introduce the effect of diffusion into the test problem, the values of this coefficient are presented in table 1.
The initial condition for the test problem presented in figure 5, has approximated values to the concentration presented in an urban center.The second reaction mechanism level is a complete reaction for the same species, and includes two more substances.The mechanism is presented in table 3.  The diffusion coefficient is also included as a factor, one produced by the molecular high concentration in a part of the domain, and other produced by the turbulence.This factor is included cause the mass diffusion could affect the reaction process and perturb the adequately simulation.For the same reason of the diffusion, the advection is included into the factors.For this process the factor has two possible maximum velocities , the low one in cases of small flow and the other one in cases of high flow.
The last factor is the most important because with this is included the option to solve the model with each one of the split and coupled method and later can compare the solution between them.
All the factors for this factorial design has two levels, then if the number of cases is 2 n where n is the number of factors, the number of cases that must be solved are 64.The response variables have been selected keeping in mind the comparison of the 2 numerical models (coupled and split), they are presented in table 4.

FACTORIAL DESIGN ANALYSIS
With all the cases solved it is possible to create graphs that present the effects of the factors over the response variables that were selected for the simulation of the mass transport phenomena.This graphs are called effects graphs and are presented below for observe the influence of the factors over the behavior of the substances NO, NO2 and O3 in the five time steps selected below.The graphs presented in the table 5 show the effects of the factors over the substances behavior in four different times.Because the raw effects can be poor indicator of outliers due to their nonconstant variance, these effects are presented in the absolute standardized form.
As seen in the table 5 the most significant effects over all the substances are the advection and diffusion process as well as its interaction, also is possible to see that the mesh affect the substance behavior and this is because the resolution of the mesh affect the precision of the model and hence the substances simulation.For the final step the effects graphs are presented in the table 7.
Table 7. Effects of the factors and their interactions over the response variables of concentration of the substances NO, NO2 and O3 in the final time step.The final concentration is affected by the same factors that affect the behavior of the substances through the simulation time.Additional to this effects the substance behavior are also affected in the mole balance of oxygen and nitrogen elements, with this is possible to know which factors affects the mole conservation of the chemical elements and most important if the reaction mechanism and the numerical method has an significant effect over this.Below is presented the effects graphs for these response variables.Table 8.Effects of the factors and their interactions over the mole balance of the oxygen and nitrogen chemical elements.

Oxygen Nitrogen SD = Standard deviation
Finally are presented the effects graphs for the CPU Time, this graphs are also important because one of the possible advantage of the split method over the coupled method is the time of simulation and with the solution of the factorial design is possible to confirm or reject this hypothesis.Table 9.Effects of the factors and their interactions over the CPU Time.
CPU Time SD = Standard deviation

ANALYSIS AND DISCUTION
When the processes of advection and diffusion are included into in the simulation the substances behavior is affected because the spatial dispersion homogenized the concentration of all the substances generating that the substances react more than without this physical process, an example of this last is presented in figure 6 and 7.Where is observe that a 600 seconds of simulations the advections and diffusion processes affect the generation rate of NO, O3 like the elimination rate of NO2.Further the interaction of these two processes presents a similar effect over the behavior reactions as seen in the table 10.
Table 10.Effects of advection-diffusion interaction over the NO, NO2 and O3 concentration in the 600 seconds of simulation.

Substance NO2
Substance NO and O3 Other factor that presents a significant effect over the substances behavior is the mesh, because depending of the mesh resolution used for solve the model the results obtained has more or less accuracy.The mesh resolution and the accuracy are directly proportional as seen in table 11, so in the cases when coarse meshes are used, the accuracy of the simulation is low.11.Effects of mesh factor over the NO, NO2 and O3 concentration in the 600 seconds of simulation.

Substance NO2
Substance NO and O3 The substances behaviors are affected largely by the interaction between the mesh, the advection and the diffusion factors.In the cases with rough mesh the advection and the diffusion processes are predicted poorly and in cases with large velocity flows and a low diffusion processes are developed high numerical oscillations that ending in a concentration sink.
The concentration sink leads to the elimination substance affecting also the mole conservation of the chemical elements.The interactions effect are presented in the table 12 and 13, and in the figure 8 is seen the mole conservation of the case 17 that present an example of a sink of concentration in the center of the domain.
Table 12.Effects of mesh-diffusion interaction over the NO, NO2 and O3 concentration in the 600 seconds of simulation.

Substance NO2 Substance NO and O3
Fine Coarse The numerical method also affect the substance behavior as seen in table 14 and 15, this mean that the two numerical methods selected for solve the problem has differences between their solution.Below are presented two tables that show mole concentration differences between the two numerical methods and its percentage obtained dividing the effect by the mean concentration of the substance for each time step.The oxygen mole balance is mostly affected by the reaction mechanism, because the simplified mechanism has an oxygen generation to satisfy the oxygen demand of the mechanism and the complete mechanism does not need any generation and hence the moles.This can be seen in the figure 9 where the simplified method has an oxygen increment of approximately 90% meanwhile the complete method does not increase the oxygen.The other factors that affect the mole balance of oxygen altered it in maximum 6.742%; this value is small for the case of simplified mechanism and in the case of complete method this percentage are insignificant because the 6.742 % of a small value like 1E-8 is irrelevant.
In the case of the nitrogen mole conservation, although some factors appear as significant, its percentage are less than 4.374 % and is presented for the presence of the high numerical oscillation talked above.
The CPU time is perturbed for some factors among these the mesh, the time step factors and their interactions as seen in figure 10.The reaction mechanism and its interactions also affect the CPU time, because the complete mechanism has more substances than simplified mechanism, and need more memory space which increase the computational cost.
The most useful effect on the CPU time is caused by the numerical method (figure 11).With this is possible to corroborate that split the nonlinear part of the system equation and solve it separately using other solution processes like in symmetrized strang split method reduce the CPU time in relation to the coupled method.

CONCLUSION
One of the major problematic of this kind of simulation is the conservation of mass specially when is included the reaction processes, however at seen previously even if the substance behavior is slightly disturbed, the mole conservation is preserved especially using the complex reaction mechanism, therefore any of the split and coupled methods are adequate to model this kind of problems.
The split method in one of the most adequate method for simulate and predict problems based in an advection-diffusion-reaction processes, because solves the nonlinear terms separately using additional solver methodologies resulting on low computational cost, small CPU time and a very small difference with the coupled method.
The behavior of the substances modeled by the mass transport equation using a finite element discretization and a SUPG methodology is affected for small oscillations especially in areas of low concentration, because the concentration gradients are more significant than in areas of high concentration.
For solve similar problems is adequately to use first a fine enough mesh and a time step among the range presented in the article to avoid generation of numerical problems like sink or oscillations effects.

Table 5 .
Effects of the factors and their interactions over the response variables of concentration of the substances NO, NO2 and O3 in the internal time step.deviation Where the factors are represented by letters as showed in table 6.

Figure 6 .
Figure 6.Mean effect of advection and diffusion over the concentration of NO and O3in the 600 seconds of simulation.

Figure 7 .
Figure 7. Mean effect of advection and diffusion over the concentration of NO2 in the 600 seconds of simulation.

Table 13 .Figure 8 .
Figure 8. Mole conservation of oxygen and nitrogen for the case 17.

Figure 9 .
Figure 9. Mean effect of reaction mechanism over the mole balance of oxygen.

Figure 10 .
Figure 10.Mean effect of mesh and time step over the CPU Time.

Figure 11 .
Figure 11.Mean effect of numerical method over the CPU Time.

Table 1 .
[15]gn factors.The first two factors are related directly with the methods used for solve the problems and are important because a change of this factor could affect the results obtained with the simulation.The first reaction mechanism level is presented in table 2, this is a simplified mechanism for NO, NO2 and O3 presented in[15]byBaik et al.

Table 14 .
Value and percentage effect of numerical method over the NO2 concentration.

Table 15 .
Value and percentage effect of numerical method over the NO and O3 concentration.