SIMULATION OF SPECIES CONCENTRATION DISTRIBUTION IN REACTIVE FLOWS WITH UNSTEADY BOUNDARY CONDITIONS

The determination of species concentration profiles in reactive flows with variable inlets is a problem of practical interest to many fields such as in flow reactor transient operation and in cyclic degradable pollutants disposals in watercourses. In these cases, the inflow condition often consists of a time-dependent function, which may imply unsteady outflows, not always well represented by the usual boundary conditions (BC) used so far. A new approach, using an outlet condition in the form of a material derivative, termed Material Derivative Boundary Condition (MDBC), is introduced and a numerical model to solve convection-diffusion-reaction equations in twodimensional (2-D) incompressible flows is developed. Upon reviewing the literature, it is noted that the Finite Element Method (FEM) is rarely used in the simulation of reactive flows, in spite of its ability of consistently coping with variable BCs. The above facts are reasons to explore its use along with a semi-discrete formulation with the Galerkin Method in our simulations. Results are obtained for various conditions, in order to show features of the code, and are compared to existing solutions. Use of the MDBC is shown to provide a better approximation of the exit concentrations and use of FEM in reactive flows is further enhanced.


INTRODUCTION Preliminaries
The determination of species concentration profiles in incompressible reactive flows presents practical interest to many engineering applications, such as tubular continuous chemical reactors design and operation, concentration evolution prediction of degradable and non-buoyant contaminants in rivers, downstream industrial wastewater or domestic sewage discharge, etc.
While reactants in chemical reactors are subject to transformation due to chemical or biochemical reactions, pollutants in rivers may also disappear by physical processes, such as volatilization or reactive decay, all of which can be accounted for in the transport equation by addition of a reaction term r (van der Perk, 2013):

Brazilian Journal of Chemical Engineering
ISSN 0104-6632 Printed in Brazil www.scielo.br/bjcewhere we define for 2-D flows: After a certain initial time interval, when the mixing processes are completed, species concentration along the flow can be modeled by the use of equation 1.In ideal tube reactors, often treated as plug flow devices, molecular diffusion and radial/lateral velocities terms may be dropped (Levenspiel, 1999), leading to a one-dimensional (1-D) pure advective-reactive model.In other cases, these terms must be taken into account, requiring 2-D models to describe the flow.It is also reasonable to assume 1-D convective and diffusive flows for small rivers and channels when the length is ten or more times larger than its width (Kachiashvili et al., 2007).In larger watercourses, in turn, where the river depth is significantly small compared to its width, depth-averaged concentrations assuming vertically well-mixed species could be employed (Lee and Seo, 2007), making it possible to apply a 2-D model derived from equation 1.
Thus, it is all about solving equation 1 in the applicable dimensions, subject to proper initial and boundary conditions.Usually, three types of BC apply: where c , q and g may be homogeneous, constant valued or functions of time and the greek letters Γ denote the corresponding surface where the BC applies.Equation 2 is usually referred to as the Dirichlet or Essential Boundary Condition (EBC), equation 3, as the Neumann or Natural Boundary Condition (NBC) and equation 4, as the Robin or Cauchy Boundary Condition.

Scope
In this paper, we are particularly interested in 2-D simulations of reacting species transport, where the inlet boundary concentration is a pulse, a series of pulses or a continuous periodic function.These inlet conditions apply to cases of flow chemical reactors operating under variable inlet feed and variable species concentration spills in rivers and channels.
This class of problems has motivated studies pursuing analytical solutions of convection-diffusion-reaction equation subjected to time-dependent BCs, like the ones from van Genuchten and Alves (1982), Logan and Zlotnik (1995), Logan (1996), Aral and Liao (1996), Golz and Dorroh (2001), Chen and Liu (2011) and Pérez Guerrero et al. (2013).However, these studies either are restricted to 1-D cases, or adopt conditions that may not represent time dependence close to the domain exit.
We emphasize that, in the case of time-dependent inlet conditions, special attention must be given to the outlet BC.Since the exit concentration or the species flux is an unknown, assuming prescribed values at the outlet is not consistent.
Up to the present, as in the works cited above, this indeterminacy is treated either by considering that the outlet concentration gradients are zero, which may be physically unrealistic (Ziskind et al., 2011), or by using Robin type BCs, best suited to represent inlet conditions.

Literature Review
A number of papers address the advection-dispersion equation, with or without the reaction term, providing both analytical and numerical solutions for cases of pollutant discharge.O' Loughlin and Bowmer (1975), for instance, applied analytical solutions to equation 1 in 1-D channel flows with decaying species, later extended by Chapman (1979) to non-uniform steady rivers, both considering only pulse or continuous inlet concentrations and homogeneous NBC for the concentrations at the outlet.Comparison with the results obtained in the experimental works of Vilhena and Leal (1981) for non-reacting pollutants in point source injection showed good agreement with them.Czernuszenko (1987), also working with dispersion of conservative species, proposed a numerical solution for the 2-D advection-diffusion equation, using a conditionally stable finite differences (FD) scheme.But, since the study was restricted to mixing far from the pollution source, leaving convection in the background, the equation was bounded by NBCs, not encompassing unsteady BCs.Piasecki and Katopodes (1997), interested in the sensitivity of contaminant concentration profiles to timely changes in their load, a similar aspect of our own concern, treated the problem by the use of a FEM scheme, but the unsteady load was a zeroth order production term of the transport equation and the problem was subjected to Dirichlet and Neumann type BCs.Kaschiashvili et al. (2007) provided a consistent model for river reactive flow problems in one, two and three dimensions and used dimensionsplitting FD numerical schemes, with unsteady upstream BC and a NBC downstream.But, due to the equilibrium condition at the outlet, consisting of a constant spatial concentration gradient, this BC no longer applies and is modified, sometimes, with the introduction of an additional parameter in order to better reproduce experimental data.The fact supports our remark that time-dependent inlet conditions may imply difficulties for prescribing values for the outlet conditions.Lee and Seo (2007) used a 2-D finite element model, based on the Streamline-Upwind Petrov-Galerkin Method (SUPG) together with a Crank-Nicholson FD scheme for the time derivative, as in this paper, but restricted to rivers where the process is diffusion dominated and the downstream BC was a prescribed diffusion flux.Two years later, the same authors employed this same method for accidental mass release in rivers (Lee and Seo, 2009).Similar to Piasecki and Katopodes (1997), the accidental mass release was represented by a zeroth order production term of the transport equation which was subjected to Dirichlet inlet BC and Neumann outlet BC, once more not considering unsteady BCs.
The literature survey detailed above, related to watercourse pollutant spills, shows that FEM has not been widely used to obtain solutions of reactive flows, in spite of its ability of consistently coping with differential BCs (Logan, 2007).This might be explained by the existence of the advective term in the transport equation that makes the system of equations nonsymmetric and prone to numerical oscillations (Yu and Singh, 1995).Several authors addressed the problem by focusing the development of consistent and stable FEM schemes for these flows (Yu and Singh, 1995;Galeão et al., 2004;John and Schmeyer, 2011) but rarely holding their attention on unsteady BCs.We also quote the studies of Konzen et al (2007) by which a convective-diffusive-reactive problem formulated through vorticity and stream-function is numerically solved, employing Galerkin FEM (GFEM) together with a Runge-Kutta scheme for the time stepping.But, owing to the formulation adopted, the BCs were assumed to be homogeneous Neumann type and the flow, taking place in a closed cavity, is not subjected to inflows and outflows rates, as in rivers and continuous chemical reactors.
Modeling work on fluid dynamics by FEM in chemical reactors is also not commonly found in the literature.Ranade's (2002) book on computational fluid modeling of reactors employs the finite volume method in the examples and applications presented.Sometimes, commercial packages using the FEM in their built-in routines are employed for the study of chemical reactors models performance (Galante, 2012;Mushtaq, 2014).However, in addition to being proprietary, these routines often focus simulations of chemical reaction media, rather than flow dynamics.Yet, it is possible to verify, in the works by Skrzypacz and Tobiska (2005) and Skrzypacz (2010), a FEM scheme to solve a simple 1-D reactive flow in packed bed reactors.Even though these two studies assume steady flow, BCs are of the Dirichlet type and the reaction term is not explicitly solved, the convenience of using FEM in chemical reactors flow modeling is pointed out.

Aims and Objective
Thus, additional motivation exists for the study of concentration fields using FEM, to simulate problems modeled by equation 1 and subjected to unsteady BCs.
Our proposal, and what depicts the main contribution of this work, is to use an outlet BC in the form of a material derivative, directly representing the concentration gradient or the species flux time dependence, an usual feature for such models.
To the authors' knowledge, no analytical solution considering a material derivative as the outlet BC was yet constructed.So, a computer code prototype is developed in MATLAB, through a semi-discrete formulation with GFEM and implicit FD scheme for the simulations.The inlet, or upstream, unsteady BC behavior is assumed either as time periodic, or as pulse functions, providing a variable condition.At the outlet, or downstream, to better represent the equilibrium condition among diffusion, advection and reaction in unsteady conditions, the outlet flux is evaluated by the species concentration material derivative.

MATHEMATICAL FORMULATION
Considering the objectives of the present study, of addressing isothermal reactive flows, an average hydrodynamic field is assumed, so turbulence models are not introduced in the evolution equations.We emphasize that averaging the concentration field along one of the three directions, in order to construct 2-D models, requires that reactants or pollutants be mixed at a much faster rate than the reaction rate, as in the microfluid idealization (Levenspiel, 1999).
The reaction term in equation 1 may vary considerably, depending on the process.For simplicity, it was decided to analyze only a first order reaction model and the diffusion tensor was considered constant.The transport equation then becomes: with initial condition given by: BCs used at the inlet or upstream are prescribed in one of the two forms below:

�
or to represent a periodic injection, we assume: where C I is the mean amplitude of the species concentration at the inlet.In equations 7-8, the y coordinate dependence is applicable to 2-D flows and may represent the injection in part or along all its length.
As already mentioned, analytical solutions for this kind of problem exist and will be used in order to validate numerical results.These solutions assume either prescribed or Neumann's outlet BCs mostly at semi-infinite domains.Moreover, even the solutions for finite domains that accept one or other of those BC are subjected to criticism (Ziskind, 2011).
Equation 5 is solved by a FEM scheme, with a Galerkin formulation.So, a weighted residual statement of that equation reads: By applying the divergence theorem to the third term of the above equation, and substituting the result in equation 9, the following weak form is obtained: . Γ 1 and Γ 2 represent lateral surfaces and the related fluxes are zero.Γ in , in turn, represents the inlet boundary, subject to specified, but time-dependent, BCs, as given by equations 7-8.In this case, the weight functions are zero for Γ in , implying that the surface integral is only evaluated along Γ out .
For the outlet surface, we can assume that: x e n   = and, therefore, the r.h.s. of equation 10 becomes: By looking again at equation 12, it can be verified that the weak formulation boundary term represents the species flux by Fick's Law.Yu and Singh (1995) sustain that this formulation should only be applied to situations where there are exclusively diffusion fluxes at the outlet boundary.But in the problems under consideration, advection effectively occurs at the outlet, and must be taken into account in the BC expression.
In fact, there are cases where gradients normal to the outlet surface are zero, bringing the formulation back into consistency, even in the presence of convection because it eliminates the surface integral.Again considering equation 12: We must have in mind that, for a developed profile, equation 13 also implies, taking into account equation 5, that: We emphasize that this condition does not hold when the gradients at the outlet are not zero.It is well known that flow problems involving the transport of chemical species with homogeneous NBC fail to satisfy the conservation law for species concentrations within the domain (Golz and Dorroh, 2001).In particular, prescribed constant outlet fluxes also do not lead to correct description of timedependent problems.
So, for the sake of generality another outlet BC must be assumed.We point out that, in the flows under consideration, the species dispersion is mainly due to vertical and transverse velocity gradients, while molecular and turbulent diffusions are generally negligible (Launay et al., 2015).So, adding the advection term to equation 14, Equation 15 is in fact a nonhomogeneous material derivative that automatically evaluates the spatial gradients at the outlet boundary.We propose to term it Material Derivative Boundary Condition, or MDBC, as previously mentioned.
Assuming that, at the boundary Γ out , , then, equation 15 can be expressed as: Combining equations 10, 12 and 16, it is then possible to write: Then, equation 17 is the one to be numerically implemented by GFEM, in order to obtain the species concentration profiles.
The numerical procedure may be tested by comparing the results with existing analytical solutions.In the simplest case of 1-D flow, analytical solutions for continuous and pulse mass injection, are, respectively (O'Loughlin and Bowmer, 1975;Chapman, 1979): and M inj is the total mass injected per unit area.And for a 2-D case with pulse injection where there is a transversal diffusion D y and zero lateral component of velocity (Vilhena and Sefidvash, 1985):


When the inlet BC is given by equation 8, an onedimensional analytical solution may be obtained.By following the work of Logan and Zlotnik (1996), it is possible to establish that equation 5 clearly admits a solution of the form: where α ˆand β ˆ are complex valued, thus: Then, substituting equations 21-22 in the 1-D form of equation 5, one obtains: Once the periodic BC forces the inlet concentration at a fixed value, β R = 0 and the solution may be expressed as: where R means the real part of equation 24 and: Also, considering that the concentration at x = 0 cannot take negative values, it is necessary to add a constant forcing, such that this restriction is satisfied, and equation 24 becomes: For this constant forcing, obviously β R = β I =0 and, therefore, with the use of equation 25: u and D x , as well as an abitrary α R , the analytical solution may be constructed, employing equations 25-28.

NUMERICAL PROCEDURE
By using the Galerkin formulation, the concentration profile is approximated by: Substituting this approximation into the weak form given by equation 17, where, according to the GFEM, the weight functions are the same as the shape functions (Zienkiewicz and Taylor, 2000), one has: where the boundary integral (r.h.s of equation 17) was approximated through: Equation 30 encompasses a stiffness matrix and a modified mass matrix, which is related to the concentration time derivative and the reaction term.It can be put under matrix form as: where M 1 and K are, respectively, the modified mass and stiffness matrices.
In order to solve equation 32, we employ a numerical scheme, using the Crank-Nicholson Method (Lewis et al., 2005), which reads: It must be observed that it is also possible to look for another solution without modifying the original mass matrix, as suggested above.In this case, the use of the Crank-Nicholson scheme on the GFEM approximation of equation 10, implies: where {B} t and {B} t+1 are the boundary terms arising from the line integral approximation on the r.h.s of equation 10 at times t and t+1, [M] is 1 and [K 1 ] is a modified stiffness matrix, now including the decay term, last term on the left of equation 17, or: In this case, the boundary vectors ({B} t and {B} t+1 ) must be evaluated using equation 31.Being dependent on the concentration and its time derivative in past and present time steps, these vectors must be continuously updated, making the numerical scheme for solving equation 33 simpler than the one required for solving equation 34.Thus, we opted for the first scheme.
The code was implemented in MATLAB, taking advantage of its matrix calculation resources.The integrals in equation 30 were evaluated by the Gauss Quadrature (GQ).The solution domain was discretized in regular triangular or quadrangular element meshes by routines within the program, depending on the case run.The program is also capable of performing GQ calculations in diversified number of interval points.Linear shape functions were used throughout this work, so precision of the scheme was controlled by properly refining the mesh.
It is well known that simple GFEM presents numerical oscillations and instabilities in problems where advection is important.So, more elaborated FEM schemes would be required to solve problems with small diffusion coefficients.However, considering that the role of the unsteady BC along (31) with the outlet BC represented by a material derivative were the main aspects to be investigated, this method was employed with restrictions.Aware that some of the major factors causing these issues are improper choice of a time step size and also of element size and shape (Yu and Singh, 1995), we adopted, as a basis for the time step and element size control, respectively, (Chapra and Canale, 2010):

Preliminary Tests
A more detailed look at the analytical solution presented by equation 18 reveals that, actually, the assumed constant upstream BC is not time independent, as it may appear to be at a first glance.Assuming unitary injection concentration (C inj = 1.0), the analytical solution results in the plots of Figure 1, obtained for Pe = 5.0 (u x = 1.0;D x = 4.0; Da = 2.0).As it can be verified, within the stream limits, the BC shows an unsteady profiles characterized by the inlet concentration correction due to particular advective and diffusion effects.Obviously, the analytical solution follows the general form of the concentration profile for this kind of problem (Vilhena and Sefidvash, 1985): where C o is the corrected species concentration at initial time.
It can be also easily seen, by inspection of equations 19 and 20, that the solution for pulse injections also follows equation 37, in order to correct the inlet concentration values.
So, in order to check the code results, the inlet BCs to be applied at x=0 must carry on the initial shape of the defined concentration, as suggested by Yu and Li (1998).This implies that: for equation 18: for equation 20:   The numerical solution of equation 5, for the periodic inlet BC (equation 8), may be compared with the 1-D analytical solution constructed from equations 25-28 through a plot extracted from the centerline concentration profile.Figure 4 shows the outcome for Pe = 100, where u x = 5.0, D x =1.0 and k = 0.1, implying Da = 0.4.In order to obtain the plots of Figures 2 to 4, we ran the code and then compared the results with the analytical solution corresponding to the time run.Numerical calculation was performed, respecting the stability restrictions posed by equations 36.The plots show good agreement between analytical and numerical solutions even for high Péclet Numbers.
It is possible to observe a better agreement between analytical and numerical solutions for the continuous injection case (equation 15), as shown in plot A of Figure 2. The plots B and C of Figure 2 (equations 19 and 20) and the plot of Figure 3 (equation 20), show that the numerical curves are slightly delayed compared to the exact solutions.This delay results from the fact that the discrete time integration cannot completely follow the instant moment of mass release (Lee and Seo, 2010).

Comparing Analytical and Numerical Solutions
Figure 5 compares simulated concentration profiles for sorted conditions, such as Pe = 5 (u x = 1.0;D x = 4.0; k = 0.1), Pe = 25 (u x = 5.0; D x = 4.0; k = 0.1) and Pe = 100 (u x = 5.0; D x = 1.0; k = 0.1).In order to obtain the plots, we solved equation 5 subjected to a time periodic inlet BC (equation 8), changing the outlet BC type.First, we employed an EBC arbitrarily set to a given constant value, then, we employed a homogeneous NBC and last, our proposed MDBC.Since the meshes used were the same in all simulations, we compared the centerline node values obtained, plotting the concentration differences (Dif C).
As we can see, profiles obtained when the adopted outlet condition is either EBC or the homogeneous NBC, compared to those obtained by the adoption of the MDBC, concentrate larger differences around the exit.
In order to check the validity of the above proposition, we numerically evaluated concentration 1-D profiles for various flow and reaction parameters.The results were compared to the analytical solution and analyzed by the Root-Mean-Square Deviation (RMSD), or: where a i C is the analytical solution at node i for a given total number of nodes nd at the exit region.

Outcome Analysis
We observe that the numerical solutions with outlet EBC provide the poorest approximations in all Péclet and Damköhler Numbers considered and that MDBC solutions result in better approximations than NBC in almost all cases.This is possibly due to the fact that MDBC better captures specific features of the flow because it encompasses, in its formulation, physical effects of the problem which are not present in the usual types of BCs.
Results in Table 1 also point at examples where the advantages of using MDBC instead of homogeneous NBC are not clear.Such situations arise from particular flow conditions that imply very small concentration gradients at the outlet, as a consequence of Péclet and Damköhler Number combinations.These cases approach patterns that can be treated conveniently by the homogeneous NBC (equation 3) and so, when we compare the outcomes obtained both with the use of NBC and MDBC, we verify analogous deviations from the analytical solution.16).
However, these are special cases of the problem and the use of the MDBC for more general formulations is established.

2-D Simulation Results
Having in mind the satisfactory results obtained in the tests, we further used the code to investigate the behavior of 2-D systems.Velocities and diffusion constants were chosen as close as possible to real configurations.
For instance, Figure 6 shows the results of 2-D and 1-D simulations under conditions such that the inlet BC is the periodic concentration oscillation given by equation 8, lateral components of velocity and diffusivity are ten times smaller than the longitudinal components (u x = 5.0, u y = 0.5, D x =1.0, D y =0.1 and k = 0.1), implying Pe = 100 and Da = 0.4.16) In this case, corresponding to a high Pe, convective transport plays a major role overcoming diffusion transport and reaction decay.Parts A and B of Figure 6 show the oscillatory behaviour of the concentration profile along the domain at different time values for the concentration along all the domain.We also note the variable outlet concentration values that would not properly be captured by EBCs and possibly NBCs.
Figure 7 shows the outlet concentrations for Pe = 10 and Da = 0.4 (u x =0.5; u y =0.05; D x =1.0; D y =0.1; k=0.01), subject to the same BCs, implying a more important role for diffusive transport.In addition, smaller flow rates allow the chemical reaction to further evolve as the convective transport takes place.The oscillatory behavior of the inlet concentration is damped before reaching the domain outlet and the solution approaches the typical shape of pure diffusive transport problems subjected to oscillatory BC, known as periodic steady-state (Bird et al, 2002).Figure 9 depicts the evolution of the species cloud deformed due to the existence of lateral components of velocity and diffusion.But the mass injection occurs uniformly at the inlet cross section area, a condition most found in chemical reactors or in small channels.
So, in order to demonstrate the code ability to simulate conditions more likely to happen in large watercourses, we modify the inlet BC as follows.Considering that in 2-D analysis the inlet may also be dependent on y (equation 8) we are able to obtain results: a) for centered point source injection: Figure 10 16)

CONCLUSION
In transient reactive flow problems subjected to unsteady BC the main issue is to achieve physical coherence in constructing the model to be solved.Some analytical solutions of this class of problems are found in the literature which, though being parabolic, usually assume that the outlet BC is in the form of a constant concentration or of a given concentration gradient.
As indicated by Piasecki and Katopodes (1997), simulations presented in this work confirmed that oscillatory inlet conditions result in time-dependent concentrations at the outlet, that cannot be accounted for by EBCs and NBCs.Also, NBCs may not represent the total equilibrium flux at the outlet (Yu and Singh, 1995), leading to physically incomplete models that could perform imprecise profile estimation.
A new procedure was then proposed, by which a material derivative is considered as the outlet BC.Our results show that these BCs provide a better picture of the process, updating the outlet equilibrium concentration.
A MATLAB code was developed with a numerical scheme subjected to prescribed stability restrictions (equations 36), using a semi-implicit GFEM scheme.Good agreement was obtained between simulations and existing analytical solutions, as can be seen on Figures 2 to 4 and Table 1.Also shown in Table 1 and Figure 5 are comparisons of numerical solutions using EBC, homogeneous NBC and the proposed MDBC, evidencing the positive aspects of applying the material derivative as the outlet BC. 2-D simulations were then performed in rectangular channels, assuming fully developed velocity profiles.
The code features a certain flexibility for automatically generating regular triangular and quadrangular meshes that could be selected for the applicable case.There was also the option of changing the number of GQ points to evaluate the model integrals, known to slightly affect the computational time.
Several simulations were run on an i5 CPU notebook, limited to a maximum of 2000 element meshes, all requiring a few minutes to run, showing that even more refined meshes could be used while keeping CPU times within acceptable limits.Our tests indicate that the numerical scheme is sufficiently tested to be implemented in codes written in lower level languages.
The use of FEM in reactive flow simulations was reinforced and, finally, a further improvement could be made in the code by future work, in the sense of adopting more elaborated FEM formulations, involving a SUPG or other more advanced stabilization technique, so as to combine the advantages of more stable schemes with the proposed adoption of the MDBC.

NOMENCLATURE
in order to represent short injections at arbitrary times nτ,
mind, one can apply equations 38-40 to the MATLAB code and compare the results with the analytical solutions for constant and pulse injection cases.In the following Figures 2-3, conditions for Pe = 5.0 are the same as for Figure 1; for Pe = 50 are: u x = 10, D x = 4.0 and k = 1.0; for Pe = 200 are: u x = 10, D x = 1.0 and k = 1.0, resulting in the same Damköhler Number (2.0) for all cases.For the tests with equation 20, which admits a lateral component of diffusion, D y was set equal to 0.2 and its 1-D plot (graph C of Figure 2) represents the centerline concentration profile (y = 0.0).
; b) for right margin (left bottom) point source injection: Figure 11; c) for left margin (upper left) point source injection: Figure 12.Figures 10 and 11 are obtained from the same parameters as those for Figure 9 and in Figure 12 the lateral component of the velocity is set from the upper margin downwards, assuming the negative of Figure 9 value for this same component.

Table 1 .
RMSD between 1-D Analytical and Numerical Solutions.