Simulation of species concentration distribution in reactive flows with unsteady boundary conditions

Abstract

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 two-dimensional (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.

Keywords:
Concentration Profile Simulation; 2-D Reactive Flows; Finite Element Method; Material Derivative; Unsteady Boundary Conditions.

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):

$$\frac{\partial C}{\partial t}=-{\overline{u}}_i\frac{\partial C}{\partial x_i}+\frac{\partial }{\partial x_i}\left[D_{ij}\frac{\partial C}{\partial x_j}\right]\pm r$$ (1)

where we define for 2-D flows:

$$D_{ij}=\left[\begin{array}{cc}{D}_x& 0\\ 0& D_y\end{array}\right]$$

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:

$$C=\overline{c}\begin{array}{ccc}& & \end{array}on{\Gamma }_e$$ (2)

$$\frac{\partial C}{\partial x_i}=\overline{q}\begin{array}{ccc}& & \end{array}on{\Gamma }_n$$ (3)

$$\overline{u}C+\overline{D}\frac{\partial C}{\partial x_i}=\overline{g}\begin{array}{ccc}& & \end{array}on{\Gamma }_r$$ (4)

where $\overline{c}$ , $\overline{q}$ and $\overline{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. Kachiashvili et al., (2007) provided a consistent model for river reactive flow problems in one, two and three dimensions and used dimension-splitting 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, 2010). 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, 2008) 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:

$$\frac{\partial C}{\partial t}=-{\overline{u}}_i\frac{\partial C}{\partial x_i}+D_{ij}\frac{{\partial }^{2}C}{\partial x_i\partial x_j}-kC$$ (5)

with initial condition given by:

$$C\left(x_i,0\right)=0$$ (6)

BCs used at the inlet or upstream are prescribed in one of the two forms below:

$$\begin{array}{l}{C}_{inj}\left(\mathrm{0,}y,t\right)=\mathrm{0,}t\ne n\tau \hfill \\ C_{inj}\left(\mathrm{0,}y,t\right)=C^{*}{}_{inj},t=n\tau \hfill \end{array}\}$$ (7)

in order to represent short injections at arbitrary times nτ, or to represent a periodic injection, we assume :

$$C_{inj}\left(\mathrm{0,}y,t\right)=C_I\left(1+cosm\pi t\right)$$ (8)

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:

$${\displaystyle \underset{\Omega }{\int }\left(\frac{\partial C}{\partial t}+{\overline{u}}_i\frac{\partial C}{\partial x_i}-D_{ij}\frac{{\partial }^{2}C}{\partial x_i\partial x_j}+kC\right)}wd\Omega =0$$ (9)

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:

$${\displaystyle \underset{\Omega }{\int }\left(w\frac{\partial C}{\partial t}+w{\overline{u}}_x.\frac{\partial C}{\partial x}+w{\overline{u}}_y.\frac{\partial C}{\partial y}+D_x\frac{\partial w}{\partial x}.\frac{\partial C}{\partial x}+D_y\frac{\partial w}{\partial y}.\frac{\partial C}{\partial y}+kwC\right)}d\Omega ={\displaystyle \underset{\Gamma }{\int }w\left(n_x.D_x\frac{\partial C}{\partial x}+n_y.\text{}{D}_y\frac{\partial C}{\partial y}\right)}d\Gamma$$ (10)

where $n_x=\overrightarrow{n}\cdot {\overrightarrow{e}}_x,\begin{array}{cc}& \end{array}{n}_y=\overrightarrow{n}\cdot {\overrightarrow{e}}_y\begin{array}{cc}& \end{array}and\begin{array}{cc}& \end{array}\Gamma ={\Gamma }_{in}\cup {\Gamma }_1\cup {\Gamma }_2\cup {\Gamma }_{out}$ .

Γ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:

$$\overrightarrow{n}={\overrightarrow{e}}_x$$ (11)

and, therefore, the r.h.s. of equation 10 becomes:

$${\displaystyle \underset{{\Gamma }_{out}}{\int }w\left(D_x\frac{\partial C}{\partial x}\right)}d{\Gamma }_{out}$$ (12)

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 :

$${\frac{\partial C}{\partial x}|}_{{\Gamma }_{out}}=0\Rightarrow {\displaystyle \underset{{\Gamma }_{out}}{\int }w\left(D_x\frac{\partial C}{\partial x}\right)}d{\Gamma }_{out}=0$$ (13)

We must have in mind that, for a developed profile, equation 13 also implies, taking into account equation 5 , that:

$${\frac{\partial C}{\partial t}|}_{{\Gamma }_{out}}=-kC$$ (14)

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 time-dependent 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 , one has:

$${\left(\frac{\partial C}{\partial t}+{\overline{u}}_i.\frac{\partial C}{\partial x_i}\right)|}_{{\Gamma }_{out}}=-{kC|}_{{\Gamma }_{out}}$$ (15)

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, $\overrightarrow{u}=\overrightarrow{n}U$, where $\overline{U}=\sqrt{{\overline{u}}_x^2+{\overline{u}}_y^2}$ , then, equation 15 can be expressed as:

$${\left(\frac{\partial C}{\partial t}+\overline{U}{n}_i\frac{\partial C}{\partial x_i}\right)|}_{{\Gamma }_{out}}=-{kC|}_{{\Gamma }_{out}}$$ (16)

where: $n_i=\overrightarrow{n}\cdot {\overrightarrow{e}}_i$

Combining equations 10 , 12 and 16 , it is then possible to write:

$${\displaystyle \underset{\Omega }{\int }\left(w\frac{\partial C}{\partial t}+w{\overline{u}}_x.\frac{\partial C}{\partial x}+w{\overline{u}}_y.\frac{\partial C}{\partial y}+D_x\frac{\partial w}{\partial x}.\frac{\partial C}{\partial x}+D_y\frac{\partial w}{\partial y}.\frac{\partial C}{\partial y}+kwC\right)}d\Omega =-{\displaystyle \underset{{\Gamma }_{out}}{\int }w\frac{D_x}{\overline{U}}\left(kC+\frac{\partial C}{\partial t}\right)}d{\Gamma }_{out}$$ (17)

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):

$$\frac{C\left(x,t\right)}{C_{inj}}=\frac{1}{2}\exp \left(\frac{-kx}{{\overline{u}}_x}\right)erfc\left[\frac{x-{\overline{u}}_{x}t\left(1+H_x\right)}{\sqrt{4D_{x}t}}\right]$$ (18)

and:

$$C\left(x,t\right)=\frac{M_{inj}}{\sqrt{4\pi D_{x}t}}\exp \left[-kt-\frac{{\left(x-{\overline{u}}_{x}t\right)}^2}{4D_{x}t}\right]$$ (19)

where $H_x=\raisebox{1ex}{$2k{D}_x$}\!\left/ \!\raisebox{-1ex}{${\overline{u}}_x^2$}\right.$ 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):

$$C\left(x,y,t\right)=\frac{M_{inj}}{4\pi t\sqrt{D_x{D}_y}}\exp -\left[kt+\frac{{\left(x-{\overline{u}}_{x}t\right)}^2}{4D_{x}t}+\frac{y^2}{4D_{y}t}\right]$$ (20)

When the inlet BC is given by equation 8 , an one-dimensional 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:

$$C\left(x,t\right)=e^{\widehat{\alpha }x+\widehat{\beta }t}$$ (21)

where $\widehat{\alpha }$ and $\widehat{\beta }$ are complex valued, thus:

$$\widehat{\alpha }={\alpha }_R+i{\alpha }_I\begin{array}{ccc}& & \end{array}and\begin{array}{ccc}& & \end{array}\widehat{\beta }={\beta }_R+i{\beta }_I$$ (22)

Then, substituting equations 21 - 22 in the 1-D form of equation 5 , one obtains:

$$\widehat{\beta }=-k-{\overline{u}}_x\widehat{\alpha }+D_x{\widehat{\alpha }}^2$$ (23)

Once the periodic BC forces the inlet concentration at a fixed value, β _R = 0 and the solution may be expressed as:

$$C\left(x,t\right)=R\left[e^{\left({\alpha }_R+i{\alpha }_I\right)x+i{\beta }_{I}t}\right]$$ (24)

where R means the real part of equation 24 and:

$${\alpha }_I=\pm \sqrt{{\alpha }_R^2-\frac{{\overline{u}}_x}{D_x}{\alpha }_R-\frac{k}{D_x}}\begin{array}{cc}& \end{array}and\begin{array}{cc}& \end{array}{\beta }_I=\pm \sqrt{{\alpha }_R^2-\frac{{\overline{u}}_x}{D_x}{\alpha }_R-\frac{k}{D_x}}\left(2D_x{\alpha }_R-{\overline{u}}_x\right)$$ (25)

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:

$$C\left(x,t\right)=C_o+R\left[e^{\left({\alpha }_R+i{\alpha }_I\right)x+i{\beta }_{I}t}\right]$$ (26)

For this constant forcing, obviously β _R = β _I =0 and, therefore, with the use of equation 25 :

$${\widehat{\alpha }}_o=-\frac{{\overline{u}}_x}{2D{}_x}\pm \sqrt{\frac{{\overline{u}}_x^2}{4D_x^2}-\frac{k}{D_x}}$$ (27)

which implies that:

$$C_o=R\left[e^{\left({\widehat{\alpha }}_o\right)x}\right]$$ (28)

Thus, given k, ${\overline{u}}_x$ 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:

$$C_{appr}\left(x_i,t\right)={\displaystyle \sum _{j=1}^{NN}{C}_j\left(t\right)S_j}\left(x_i\right)$$ (29)

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:

$${\displaystyle \sum _{j=1}^{NN}\{\left({\displaystyle \underset{\Omega }{\int }{S}_i{S}_j}d\Omega +\frac{D_x}{\overline{U}}{\displaystyle \underset{{\Gamma }_{out}}{\int }{S}_i{S}_{j}d{\Gamma }_{out}}\right)\frac{d{C}_j}{dt}+{\displaystyle \underset{\Omega }{\int }[S_i\left({\overline{u}}_x\frac{\partial S_j}{\partial x}+{\overline{u}}_y\frac{\partial S_j}{\partial y}\right)+(D_x\frac{\partial S_i}{\partial x}\frac{\partial S_j}{\partial x}+}}{D}_y\frac{\partial S_i}{\partial y}\frac{\partial S_j}{\partial y})]d\Omega C_j+k\left({\displaystyle \underset{\Omega }{\int }{S}_i{S}_{j}d\Omega }+\frac{D_x}{\overline{U}}{\displaystyle \underset{{\Gamma }_{out}}{\int }{S}_i{S}_{j}d{\Gamma }_{out}}\right)C_j\}=0$$ (30)

where the boundary integral (r.h.s of equation 17 ) was approximated through:

$${\displaystyle \sum _{j=1}^{NN}\left[\frac{D_x}{\overline{U}}\left({\displaystyle \underset{{\Gamma }_{out}}{\int }{S}_i{S}_{j}d{\Gamma }_{out}}\right)\left(-k{C}_j+\frac{d{C}_j}{dt}\right)\right]}$$ (31)

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:

$$\left[{\text{M}}_{\text{1}}\right]\left\{\stackrel{•}{C}\right\}+\left[\text{K}\right]\left\{C\right\}+k\left[{\text{M}}_{\text{1}}\right]\left\{C\right\}=0$$ (32)

where M₁ 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:

$$\left\{C_{t+1}\right\}={\left(\left[{\text{M}}_{\text{1}}\right]+\frac{\Delta t}{2}\left(\text{K}-k{\text{[M}}_{\text{1}}\text{]}\right)\right)}^{-1}\left(\left[{\text{M}}_{\text{1}}\right]-\frac{\Delta t}{2}\left(\text{K}+k{\text{[M}}_{\text{1}}\text{]}\right)\right)C_t$$ (33)

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:

$$\left\{C_{t+1}\right\}={\left(\left[\text{M}\right]+\frac{\Delta t}{2}\left[{\text{K}}_{\text{1}}\right]\right)}^{-1}\left[\left(\left[\text{M}\right]-\frac{\Delta t}{2}\left[{\text{K}}_{\text{1}}\right]\right)C_t+\frac{\Delta t}{2}\left({\left\{B\right\}}_t+{\left\{B\right\}}_{t+1}\right)\right]$$ (34)

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 ${\displaystyle \sum _{j=1}^{NN}{\displaystyle \underset{\Omega }{\int }{S}_i{S}_j}d\Omega }$ and [K₁] is a modified stiffness matrix, now including the decay term, last term on the left of equation 17 , or:

$${\displaystyle \sum _{j=1}^{NN}\{{\displaystyle \underset{\Omega }{\int }[S_i\left({\overline{u}}_x\frac{\partial S_j}{\partial x}+{\overline{u}}_y\frac{\partial S_j}{\partial y}\right)+(D_x\frac{\partial S_i}{\partial x}\frac{\partial S_j}{\partial x}+}}{D}_y\frac{\partial S_i}{\partial y}\frac{\partial S_j}{\partial y}+k{S}_i{S}_j)]d\Omega C_j\}$$ (35)

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 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):

$$\Delta t_i\le \frac{{\left(\Delta x_i\right)}^2}{2D_{x_i}+k.{\left(\Delta x_i\right)}^2}\begin{array}{ccc}& and& \end{array}\Delta x_i\le \frac{2D_{x_i}}{{\overline{u}}_i}$$ (36)

RESULTS AND DISCUSSION

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 (${\overline{u}}_x$ = 1.0; D _x = 4.0; Da = 2.0).

Figure 1
1-D Plot of Analytical Solution ( Equation 18 ).

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):

$$C\left(x,t\right)=C_{\text{o}}\left(x,t\right)\exp \left(-kt\right)$$ (37)

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 :

$$C_{\text{o}}\left(\mathrm{0,}t\right)=\frac{C_{inj}}{2}\left\{erfc\left[\frac{-{\overline{u}}_{x}t\left(1+H_x\right)}{\sqrt{4D_{x}t}}\right]\right\}$$ (38)

for equation 19 :

$$C_{\text{o}}\left(\mathrm{0,}t\right)=\frac{M_{inj}}{\sqrt{4\pi D_{x}t}}\exp \left[-kt-\frac{{\left(-{\overline{u}}_{x}t\right)}^2}{4D_{x}t}\right]$$ (39)

and:

for equation 20 :

$$C_{\text{o}}\left(\mathrm{0,}y,t\right)=\frac{M_{inj}}{4\pi t\sqrt{D_x{D}_y}}\exp \left[-kt-\frac{{\left(-{\overline{u}}_{x}t\right)}^2}{4D_{x}t}-\frac{y^2}{4D_{y}t}\right]$$ (40)

Having that in 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: ${\overline{u}}_x$ = 10, Dx = 4.0 and k = 1.0; for Pe = 200 are: ${\overline{u}}_x$ = 10, Dx = 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, Dy was set equal to 0.2 and its 1-D plot (graph C of Figure 2) represents the centerline concentration profile (y = 0.0).

Figure 2
1-D Analytical and Numerical Solution of Equations 18 - 20 Cases.

Figure 3
2-D Analytical and Numerical Solution of Equation 20 Case. (1250 Element Mesh; GQ 9 points; Pe = 50)

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 ${\overline{u}}_x$ = 5.0, Dx =1.0 and k = 0.1, implying Da = 0.4.

Figure 4
1-D Analytical and Numerical Solutions of Equation 5 with periodic inlet BC. (1250 Element Mesh; GQ 9 points).

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 (${\overline{u}}_x$ = 1.0; Dx = 4.0; k = 0.1), Pe = 25 (${\overline{u}}_x$ = 5.0; Dx = 4.0; k = 0.1) and Pe = 100 (${\overline{u}}_x$ = 5.0; Dx = 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.

Figure 5
Centerline Concentration Profile Differences - Numerical Solutions.

(Inlet BC: Equation 8 ; Outlet EBC = 0.5; Outlet NBC: Equation 13 ; Outlet MDBC: Equation 16 ).

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:

$$\text{RMSD}=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^{nd}{\left(C{}_i-C_i^a\right)}^2}}{nd}}$$ (41)

where $C_i^a$ is the analytical solution at node i for a given total number of nodes nd at the exit region.

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

Pe = 100			RMSD
Δx	Δt	Da	An. ˗ EBC	An. - NBC	An. - MDBC
0.2	0.02	0.1	0.8193	0.0315	0.0226
		1.0	0.2640	0.1839	0.1744
		2.0	0.0745	0.0044	0.0022
Pe = 50			RMSD
0.2	0.005	0.1	0.8628	0.0098	0.0041
	0.05	1.0	0.4982	0.0106	0.0032
		2.0	0.1312	0.0809	0.0798
Pe = 25			RMSD
0.2	0.02	0.1	0.8846	0.0777	0.0537
0.4	0.01	1.0	0.4010	0.0676	0.0679
0.2	0.02	2.0	0.1476	0.0387	0.0259
Pe = 5			RMSD
0.2	0.05	0.1	0.6155	0.0275	0.0191
	0.2	1.0	0.5672	0.0071	0.0072
	0.1	2.0	0.0880	0.0178	0.0034

• (Inlet BC: Equation 8 ; Outlet EBC = 0.0; Outlet NBC: Equation 13 ; Outlet MDBC: Equation 16 ).

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. 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 (${\overline{u}}_x$ = 5.0, ${\overline{u}}_y$ = 0.5, Dx =1.0, Dy =0.1 and k = 0.1), implying Pe = 100 and Da = 0.4.

Figure 6
Concentration Profile for Decaying Species (900 Element Mesh - Inlet BC: Equation 8 ; Outlet BC: Equation 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 (${\overline{u}}_x$ =0.5; ${\overline{u}}_y$ =0.05; Dx =1.0; Dy =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 7
Concentration Profile for Decaying Species (2000 Element Mesh - Inlet BC: Equation 8; Outlet BC: Equation 16).

When equations 7 are applied as the inlet BC, resulting in pulse injection of time-dependent concentrations, the code shows the concentration profiles approaching the oscillatory profile as the interval time between each injection becomes shorter (part A of Figure 8), or the pulse injection profile (part B of Figure 8), in a Gaussian shape, as it becomes larger.

Figure 8
Concentration Profile for Decaying Species (1000 Element Mesh - Inlet BC: Equation 7 ; Outlet BC: Equation 16).

The code is able to simulate 2-D configurations. including flow predictions when the velocity profiles are steady but dependent on the spatial coordinates, such that ${\overline{u}}_x={\overline{u}}_x\left(y\right)$ and ${\overline{u}}_y={\overline{u}}_y\left(x\right)$. For example, if a steady parabolic profile is considered for the longitudinal velocity ( equation 42 ), for the same other parameters as those of Figure 4, Figure 9 is obtained:

$${\overline{u}}_x=5.0y-y^2$$ (42)

Figure 9
Concentration Profile for Decaying Species - Parabolic Longitudinal Velocity (700 Element Mesh - Inlet BC: Equation 8; Outlet BC: Equation 16)

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:

for centered point source injection: Figure 10;

for right margin (left bottom) point source injection: Figure 11;

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.

Figure 10
Concentration Profile for Decaying Species - Left centerline injection (1250 Element Mesh - Inlet BC: Equation 8; Outlet BC: Equation 16)

Figure 11
Concentration Profile for Decaying Species - Bottom left injection (1250 Element Mesh - Inlet BC: Equation 8; Outlet BC: Equation 16)

Figure 12
Concentration Profile for Decaying Species - Upper left injection (1250 Element Mesh - Inlet BC: Equation 8; Outlet BC: Equation 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

• C  - section-averaged species concentration
• C _{appr  -} approximate concentration given by the FEM formulation
• C _{inj  -} injected averaged concentration
• Dx, D _{y  -} averaged diffusion coefficient in the direction of the respective coordinate axis
• Da  - Damköhler Number
• k  - reaction constant or pollutant decay constant
• m, n  - arbitrary integers 1, 2, 3 …
• NN  - number of nodes in the finite element mesh
• Pe  - Péclet Number
• r  - reaction term
• S_j (x _i )  - shape function
• t  - time
•  - averaged flow velocity along coordinate x_i
• w  - arbitrary weight function
•  - coordinate in an arbitrary direction i
• Γ  - control surface
• Γ_s  - arbitrary boundary on surface s
• τ  - arbitrary time between injections
• Ω  - control volume

ACKNOWLEGEMENTS

The authors wish to thank the Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ) for its support of this research under the program number E-26/111.079/2013.

REFERENCES

• Aral, M. M. and Liao, B. Analytical Solution for Two-Dimensional Transport Equation with Time-Dependent Dispersion Coefficients., J. Hydrol. Engrg., 1, n^o 1, 20-32 (1996).
• Bird, R. B., Stewart, W. E. and Lightfoot, E. N., Transport Phenomena, 2^nd ed. John Wiley & Sons, New York (2002).
• Chapman, B. M., Dispersion of Soluble Pollutants in Non-uniform Rivers - I. Theory, J. Hydrol., 40, 139-152 (1979).
• Chapra, S.C and Canale, R. P., Numerical Methods for Engineers, 6^th ed. Mc Graw-Hill, New York (2010).
• Chen, J. S. and Liu, C. W., Generalized Analytical Solutions for Advection-Dispersion Equation in Finite Spatial Domain with Arbitrary Time-Dependent Inlet Boundary Condition. Hydrol. Earth Syst. Sci., 15, 2471-2479 (2011).
• Czernuszenko, W., Dispersion of Pollutants in Rivers, Hydr. Sci. J., 32, n^o 1, 59-67 (1987).
• Galante, R. L., Modeling and Simulation of a Continuous Tubular Reactor for the Biodiesel Production (in Portuguese), D.S.Thesis, University of Santa Catarina, Santa Catarina (2012).
• Galeão, A. C., Almeida, R. C., Malta, S. M. C. and Loula, A. F. D., Finite Element Analysis of Convection Dominated Reaction-Diffusion Problems, Appl. Num. Math, 48, 205-222 (2004).
• Golz, W. J. and Dorroh, J. R. The Convection-Diffusion Equation for a Finite Domain with Time Varying Boundaries, Appl. Math. Lett., 14, 983-988 (2001).
• John, V. and Schmeyer, E., Finite Element Methods for Time-Dependent Convection-Diffusion-Reaction Equations with Small Diffusion, Comp. Meth. Appl. Mech. Engrg, 198, 475-494 (2008).
• Kachiashvili, K., Gordeziani, D, Lazarov, R. and Melikdzhanian, D., Modelling and Simulation of Pollutants Transport in Rivers, Appl. Math. Mod., 31, 1371-1396 (2007).
• Konzen, P. H., De Bortoli, A. L. and Thompson, M., Finite Element Method Applied to the Solution of a Convective-Diffusive-Reactive Flow, In: XXX Congresso Nacional de Matemática Aplicada e Computacional - XXX CNMAC (Annals), Florianópolis (2007). Available at: http://www.sbmac.org.br/eventos/cnmac/xxx_cnmac/PDF/294.pdf
  » http://www.sbmac.org.br/eventos/cnmac/xxx_cnmac/PDF/294.pdf
• Launay, M., Le Coz, J., Camenen, B., Walter, C., Angot, H., Dramais, G., Faure, J.-B. and Coquery, M., Calibrating Pollutant Dispersion in 1-D Hydraulic Models of River Networks, J. Hydro-env. Res., 9, 120-132 (2015).
• Lee, M. E. and Seo, I. W., Analysis of Pollutant Transport in the Han River with Tidal Current Using a 2D Finite Element Model, J. Hydro-env. Res. , 1, n^o1, 30-42 (2007).
• Lee, M. E. and Seo, I. W., 2D Finite Element Pollutant Transport Model for Accidental Mass Release in Rivers, KSCE J. Civ. Eng., 14, n^o1, 77-86 (2010).
• Levenspiel, O., Chemical Reaction Engineering, 3^rd ed. John Wiley & Sons, New York (1999).
• Lewis, R. W., Nithiarasu, P. and Seetharamu, K. N., Fundamentals of the Finite Element Method for Heat and Fluid Flow. John Wiley & Sons, Chichester (2005).
• Logan, D. L. A First Course in the Finite Element Method - 4^th Edition. Thomson, Ontario (2007).
• Logan, J. D. Solute Transport in Porous Media with Scale-Dependent Dispersion and Periodic Boundary Conditions, J. Hydr., 184, n^o 3, 261-276 (1996).
• Logan, J. D. and Zlotnik, V. The Convection-Diffusion Equation with Periodic Boundary Conditions, Appl. Math. Lett. , 8, n^o 3, 55-61 (1995).
• Mushtaq, F., Analysis and Validation of Chemical Reactors Performance Models Developed in a Commercial Software Platform, M.S. Thesis, KTH School of Industrial Engineering and Management, Stockholm (2014).
• O’Loughlin, E. M. and Bowmer, K. H., Dilution and Decay of Aquatic Herbicides in Flowing Channels, J. Hydrol. , 26, 217-235 (1975).
• Pérez Guerrero, J. S., Pontedero, E. M., van Genuchten, M. Th. and Skaggs, T. H., Analytical Solutions of the One-Dimensional Advection-Dispersion Solute Transport Equation Subject to Time-Dependent Boundary Conditions. Chem. Eng. J., 221, 487-491 (2013).
• Piasecki, M. and Katopodes, N. D., Control of Contaminant Releases in Rivers I: Adjoint Sensitivity Analysis, J. Hydrau. Eng., 123, n^o 6, 486-492 (1997).
• Ranade, V. V., Computational Flow Modeling for Chemical Reactor Engineering. Academic Press, San Diego (2002).
• Skrzypacz, P. and Tobiska, L., Finite Element and Matched Asymptotic Expansion Methods for Chemical Reactor Flow Problems, Proc. Appl. Math. Mech., 5, 843-844 (2005).
• Skrzypacz, P., Finite Element Analysis for Flow in Chemical Reactors , Dr. Rer. Nat. Dissertation, Otto von Guericke Universität, Magdeburg (2010).
• van der Perk, M., Soil and Water Contamination, 2^nd ed. CRC/Balkema, Leiden (2013).
• van Genuchten, M. Th. and Alves, W. J. Analytical Solutions of the One-Dimensional Convective-Dispersive Solute Transport Equation, US Dept. of Agric. Tech. Bull, no. 1661, 151 (1982).
• Vilhena, M.T. and Leal, C. A., Dispersion of Non-Degradable Pollutants in Rivers, Int. J. Appl. Radiat. Isot., 32, 443-446 (1981).
• Vilhena, M.T. and Sefidvash, F., Two-dimensional Treatment of Dispersion of Pollutants in Rivers, Int. J. Appl. Radiat. Isot. , 36, n^o 7, 569-572 (1985).
• Yu, F. X. and Singh, V. P., Improved Finite Element Method for Solute Transport, J. Hydraul. Eng., 121, n^o 2, 145-158 (1995).
• Yu, T.S. and Li, C. W., Instantaneous Discharge of Buoyant Fluid in Cross-flow, J. Hydraul. Eng. , 124, n^o 1, 1161-1176 (1998).
• Zienkiewicz, O. C. and Taylor, R. L., The Finite Element Method, Vol 1: The Basis, 5^th ed. Butterworth-Heinemann, Oxford (2000).
• Ziskind, G., Shmueli, H. and Gitis, V., An Analytical Solution of the Convection-Dispersion-Reaction Equation for a Finite Region with a Pulse Boundary Condition, Chem. Eng. J. , 167, 403-408 (2011).

Publication Dates

• Publication in this collection
  Oct 2017

History

• Received
  21 Jan 2016
• Reviewed
  24 June 2016
• Accepted
  12 July 2016