Mathematical modeling of heterogeneous biological ... - UPCommons

4 downloads 0 Views 5MB Size Report
C6H12O6 + 6O2 ⇒ 6CO2 + 6H2O ... C6H12O6 ⇒ 2CO2 + 2CH3CH2OH ... rial to energy-rich carbon polymers called polyhydroxyalkanoates without growing.
Mathematical modeling of heterogeneous biological wastewater processes Jordi Pascual i Ferrer

ii

Chapter 1 Motivation These days, wastewater treatment is a usual practice. Drainage and stormwater collection started in ancient times, so did the collection of the wastewater. Numerous examples of collection can be found through the history. Since the 20th century b.C., when the Hindus started using them, to the present time, many civilizations have used them (the ancient Egypt, the Greeks, the Romans, etc) [14] . But it was not until the late 1800s and early 1900s when its treatment, as it is known today, began to be made. On this treatment, played an important role the germ theory developed by Koch and Pasteur in the latter half of the nineteenth century [2]. There are different ways on doing a wastewater treatment, but all of them have the aim of returning water after its use to the environment without the pollutants which can damage it and our health. Those are traditionally classified as primary, secondary and tertiary treatment. The first one refers to physical processes which take place in order to separate the solids suspended in water. The secondary treatment, which can be made by biological or chemical processes, removes most of the organic matter. And the last treatment is used to remove other constituents that are not significantly reduced on the previous ones. The numerical modeling of the subsurface flow constructed wetlands is the motivation of this master thesis. They are part of the named “natural systems” for wastewater treatment. They consist on lagoons or channels, which are planted with humid zone plants. In those wetlands water is treated simultaneously by physical, chemical and biological processes [13]. Due to the combination of these processes, that are common on the “natural systems”, water quality is similar to or even better than that coming out from a tertiary treatment [2].

2

1.1

Mathematical modeling of heterogeneous biological wastewater processes

How did all begin

At Secci´o d’Enginyeria Sanit`aria i Ambiental (Sanitary and Environmental Engineering Section) which belongs to the Departament d’Enginyeria Hidr`aulica, Mar´ıtima i Ambiental (DEHMA, Hydraulic, Marine and Environmental Engineering) at UPC, one of their research topics is related to the natural wastewater treatment systems to treat urban wastewater. Their present work is based on the study of the design parameters from subsurface flow constructed wetlands so that secondary and tertiary treatment can be there correctly achieved. Due to this research line, they work on a pilot plant which has eight different wetlands (about 500 m2 altogether) where secondary treatment is made to the wastewater which comes from an urbanization (about 150 inhabitants) located on Les Franqueses del Vall`es (Barcelona). There, studies are based on the influence of the wetlands shape, its grain size distribution as also trying to find out the main routes of the organic matter degradation. Concerning to those main routes, on the superior part of those wetlands, aerobic processes take place, while on the inferior part, anaerobic processes are the ones in action. One of the problems when trying to understand the wetlands operation is knowing whether aerobic or anaerobic conditions are given at each point of the wetlands. It is not a simple question, though oxygen not only enters the wetlands solved in water and through the free surface, but also through the plants existing on those wetlands. But different aspects of the processes given on those wetlands are still a black box [13]. Although a biological model of wastewater processes is already developed, and some mathematical modeling of it has been achieved, those math models are not valid for its use on these subsurface flow wetlands. In order to clarify some of these aspects, the Section contacted with the Laboratori de C`alcul Num`eric (LaC`aN, Numerical Computation Laboratory) at the Departament de Matem`atica Aplicada III (MAIII, Department on Applied Mathematics III) also at UPC. From their collaboration is expected that a new mathematical model is implemented so that unsolved questions related to subsurface flow wetlands can be responded.

Motivation

1.2

3

Outline of the master thesis

This master thesis is an initial approach to those mathematical models, trying to find out which ways shall be taken to reach a satisfactory answer to all those questions. Firstly, a brief review of the biochemical processes and its modeling will be made, highlighting models presented by the International Water Association. Also the physical procedure used normally on modeling those processes, as also some other aspects refereing to its physics, will be outlined. Next the mathematical model shall be presented, starting with an overview of the continuum mechanics concepts that may lead to a mathematical formulation of the physical processes. Special attention on the boundary conditions will be made. At the end of chapter 3 the complete mathematical expression that should be solved will be presented. Once the problem statement is clear, a numerical solver shall be developed for its resolution. As it is an unsteady problem, both spatial and time discretization shall be made. The method employed on these procedures will be explained, as also some stabilization techniques that may be used in order to obtain a better solution. On chapter 5 the model will be validated, it means observing and analyzing if the whole process has been done correctly. Also there, with first results on hand, a study of the adequate stabilization parameter will be presented. Although the model may be still not so accurate to obtain realistic results, some of the obtained by the process will be shown. Last but not least, conclusions shall be presented, as also the future work which should be made to obtain a better model will be outlined.

4

Mathematical modeling of heterogeneous biological wastewater processes

Chapter 2 State of the art The present chapter tries to make an approach to the formulation of processes which take place on the subsurface flow wastewater treatment wetlands. This refers not only to the biochemical processes but also to some physical aspects. Firstly, main biochemical processes will be presented as a brief introduction on these topics, which may lead to get in touch with those concepts that shall be modeled. Then, one of its different models, known as Activated Sludge Model No. 1, will be explained. It may be remarked that this model is one of the basic tools of this master thesis. Also a brief vision on the traditional procedure used to model the hydraulics which govern those wetlands shall be presented, with the common simplifications that are normally made. Finally, as wastewater flows through a porous medium, some aspects of adsorption and its mathematical formulation shall be explained to know its effects on the resolution.

2.1

Biological wastewater treatment

Biological wastewater treatment is a set of processes that use the enzymatic and metabolic action of microorganisms present in wastewater to clean it. In contraposition to the only-chemical processes -where it is necessary to use large quantities of reagent, good systems for dosage and qualified staff- the biological processes have enough inertia to assume the concentration modifications of water or other situations appeared during the process by themselves.

Mathematical modeling of heterogeneous biological wastewater processes

6

Next, some of the main reactions given in a biological treatment are going to be presented. Those are carbon oxidation, nitrification, denitrification and phosphorus removal.

2.1.1

Carbon oxidation

Under carbon oxidation, it is understood the reactions that reduce the organic material. This degradation is done by microorganisms and can take place under aerobic and anaerobic 1 conditions. Under aerobic conditions, the reaction is ruled by the equation [13]: C6 H12 O6 + 6O2 ⇒ 6CO2 + 6H2 O

(2.1)

Otherwise, when anaerobic conditions are given, the carbon oxidation becomes a multistage process. Firstly, the complex molecules are transformed into acetic acid, lactic acid, ethanol and gases like carbon dioxide and hydrogen as the equations show [13]: C6 H12 O6 ⇒ 3CH3 COOH + H2

(2.2)

C6 H12 O6 ⇒ 2CH3 CHOHCOOH

(2.3)

C6 H12 O6 ⇒ 2CO2 + 2CH3 CH2 OH

(2.4)

On the second phase, the resultant molecules from these processes are also transformed depending on the substrate given. There are three different possibilities: methanogenesis, sulphate reduction and denitrification. Methanogenesis 4H2 + CO2 ⇒ CH4 + 2H2 O

(2.5)

CH3 COOH + 4H2 ⇒ 2CH4 + 2H2 O

(2.6)

Sulphate reduction 2CH3 CHOHCOOH + H2 SO4 ⇒ 2CH3 COOH + 2CO2 + 2H2 O + H2 S CH3 CHOHCOOH + H2 SO4 ⇒ 2CO2 + 2H2 O + H2 S

(2.7) (2.8)

1 anaerobic conditions does it mean that neither oxygen nor nitrate are present, in contrast to anoxic conditions where only absence of oxygen takes place

State of the art

7

Denitrification C6 H12 O6 + 4N O3− ⇒ 6H2 O + 6CO2 + 2N2 + 4e−

(2.9)

It must be taken under consideration that the most efficient process is given under aerobic conditions [13].

2.1.2

Nitrification

Nitrification is a two-step biological conversion of ammonia to nitrite and then to nitrate under aerobic conditions [9]. Two different groups of bacteria are needed for this process: Nitrosomonas and Nitrobacter. The stoichiometric equations that define those processes are N H4+ + 1.5O2 ⇒ N O2− + 2H + + H2 O

(2.10)

which expresses the oxidation of ammonium to nitrite done by Nitrosomonas and N O2− + 0.5O2 ⇒ N O3−

(2.11)

that expresses the oxidation of nitrite to nitrate by Nitrobacter. The overall equation of the nitrification process is then N H4+ + 2O2 ⇒ N O3− + 2H + + H2 O

2.1.3

(2.12)

Denitrification

On the other hand, denitrification takes place under anoxic conditions and is the process that reduces nitrate and nitrite into nitrogen gas. This reactions are developed by a wide range of heterotrophic bacterial species. [9] This can be accomplished with different carbon sources like methanol, ethanol, acetic acid, food processing organic waste materials or organics present in wastewater. If methanol is used, the stoichiometric equation is: 6N O3 + 5CH3 OH + H2 CO3 ⇒ 3N2 + 8H2 O + 6HCO3−

(2.13)

8

2.1.4

Mathematical modeling of heterogeneous biological wastewater processes

Phosphorus removal

Otherwise, phosphorus removal through biological processes needs the alternation of aerobic and anaerobic stages to favor biophosphorus or phosphorate-accumulating organisms (PAOs), which are the heterotrophic organisms responsible for the biological phosphorus removal. In the anaerobic stage, biophosphorus convert readily available organics material to energy-rich carbon polymers called polyhydroxyalkanoates without growing. But is in the aerobic stage where the PAOs can oxidize those polymers and use the energy and the carbons for their growth and maintenance requirements.[9]

2.2

Biochemical modeling

Modeling and simulation of wastewater treatment plants has been found important as these treatments are day by day more in use. Nowadays, there are different models of the biological treatment systems, normally developed not only for the treatment plant design and operation but also for the research applied on understanding these processes. But the best known mechanistic model are the ones developed by the International Water Association (IWA) [9]. A mechanistic model quantifies the biological processes using the actual or believed physical, chemical and microbiological expressions -in contrast of an empirical model, based upon a mathematical function that reasonably represents data from the system-. These models considers both kinetics and the stoichiometry of reactions. They try to represent the biochemical reactions within the activated sludge system and how they act on its different compounds. But modeling biological processes forces to make many simplifying assumptions. On the one hand, there are the different simplifications that appear on considering the microbial population of the wastewater, which comprises different individual organisms of multiple species, and its behavior as a unique specie whose behavior is an average of the whole. On the other hand, there are still some aspects of the non-linear dynamics and properties of the biological processes that need to be understood [6]. Those processes are not as the traditional mechanical and electrical systems, which equations

State of the art

9

are directly taken from physical laws. Moreover, models do not express all the processes given because of the large number of reactions that take place on them. Models might be made by those processes which are the most important but trying to keep the most simple expression that can allow a simple solution.

2.2.1

IWA models

In 1982 a Task Group on Mathematical Modeling for Design and Operation of Activated Sludge Processes was created in the bossom of the then called International Association on Water Pollution Research Control (IAWPRC). This Task Group had the aim of developing a model that could represent the reactions given on the activated sludge processes. Since then, four different models have been released. The first one, named Activated Sludge Model No. 1 (ASM1) was achieved on 1987. It introduced the matrix notation, where both, kinetics and stoichiometry are represented. There, carbon oxidation, nitrification and denitrification are modeled. But on mid-nineties the principles of the phosphorus removal on wastewater was being understood. For this reason on 1995 it was released a new model, the Activated Sludge Model No. 2 (ASM2), which also incorporated this phenomena. It was a few years later, as denitrifying phosphorus-accumulating organisms (PAOs) were needed and their processes understood, that the ASM2 was extended in the ASM2d, where those were included (1999). But because of the complexity that those models had taken, on 2000 a new modeling platform was developed in order to create a tool for the next generation of activated sludge models. It was the Activated Sludge Model No. 3. As complexity grows through all these different models, this master thesis is going to work around the first model, the ASM1. Employing this model may let, not only find out how might the equations be solved, but also observe how do they work. The posterior ampliation to the other IWA models should not bring much problems.

2.2.2

Activated Sludge Model No. 1

As it is said before, the ASM1 models the carbon oxidation, nitrification and denitrification. The model is expressed through a matrix where both, kinetics and stoi-

10

Mathematical modeling of heterogeneous biological wastewater processes

chiometry are represented. After taking a glance on the matrix (fig. 2.1), its elements and composition are explained.

Figure 2.1: Activated Sludge Model No.1 matrix. Source: IWA On top of the matrix, the symbols for the different components that take place on these processes can be found. It should be distinguished between those components represented by an X and those by an S. The first ones are those that are particulate while the others are soluble. Otherwise, on the bottom of the matrix, it can be read the name of the component as also units employed. On the leftmost column, the fundamental processes are written, while their rate expressions can be read on the rightmost. Those are eight processes: growth of biomass, separated on aerobic growth of heterotrophs, anoxic growth of heterotrophs and aerobic growth of autotrophs; decay of biomass, also separated on the decay of heterotrophs and autotrophs; ammonification of organic nitrogen and hydrolysis of

State of the art

11

particulate organics entrapped in the biofloc. The rate equations are denoted ρj , where j expresses the number of the process. Otherwise, within the matrix there are the stoichiometric coefficients, which express the mass relationships between the components on each process. Those coefficients are denoted as νi,j , where i expresses the component which is involved on the process, while j, as before, expresses which process is taken under consideration. On this matrix, the sign convention used is positive for production and negative for consumption (notice that oxygen is expressed as negative oxygen demand). Once the elements on the matrix are presented, lets explore a bit more the rate expressions ρj . These expressions refer to the kinetics of the process. The model on the kinetics employs the switching functions concept. A switching function is a function that allows turning on or off mathematically a process. This is really important on those processes that depend on the existence of the electron acceptor. The switch functions used are of the kind f (S) =

S K +S

(2.14)

where S is the concentration of the component which switches on or off the process, while K is the half-saturation coefficient. A lower value of K implies a more rapid change of the slope. On figure 2.2 the function f (S) is presented.

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

0.2

0.4

0.6

0.8

1

S

Figure 2.2: Representation of a switching function

Mathematical modeling of heterogeneous biological wastewater processes

12

This function, known also as Monod equation, can be rewritten as µ ¶ 1 f (S) = · S = α(S) · S K +S

(2.15)

where α(S) represents the velocity of reaction (production or destruction) of S. On figure 2.3 this velocity is represented.

0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0

0

0.2

0.4

0.6

0.8

1

S

Figure 2.3: Representation of the velocity of reaction α It may be seen that velocity of reaction has a maximum on S = 0. About zero concentration of S implies maximum production of S. Otherwise, it tends to zero when S goes to infinite, so for higher values of concentration, littler velocity of production. So the ways on reading this matrix are two: following the rows all the processes are being expressed, and attending to the columns it can be seen the evolution of the components. Our main interest is this second lecture. And an interesting point of the matrix is the easiness on representing this. It can be seen that the conversion rates of each component may be written as X ri = νi,j ρj (2.16) j

So, as explained before, the conversion rate of a component appearing on the matrix is the sum of every stoichiometric coefficient multiplied by the process rate to which it belongs. This way of expressing the different process is an easy and clear

State of the art

13

way to represent the reaction processes given. Otherwise, it must be remarked that the oxygen equation may have, in addition to the conversion rate ri , a term refereing to the oxygen transfer from gaseous air to solved oxygen. The change on the oxygen concentration may be presented by [7] KL a · (SO,sat − SO )

(2.17)

where KL a is known as Oxygen Transfer Rate (OTR). It can be described by the function ¡ ¢ KL a = k1 · 1 − ek2 ·Qair (2.18) where k1 and k2 are parameters which should be estimated [16]. As oxygen will be present in some parts of the domain of the problem, models will be solved in both conditions: allowing entrance of oxygen (so considering this term on oxygen’s equation) and not allowing its entrance (without adding this term on the matrix model). For the implementation of the numerical solution (which is going to be explained on next chapter), it has been thought to start with a reduced model before using the whole formulation of the ASM1. This was planned to facilitate the initial work, so that it would be easier to implement and to prove if this had been done satisfactorily. It has also been considered important to introduce a model with the oxygen as a component because, as it has been already explained, are the problems associated to oxygen the initial motivation of this work. Also, because oxygen presence in water is fundamental to define which processes take place. So it was decided to use a six equation model which represents an aerobic media where only COD removal and nitrification are given [1]. This involves four different process, which are the aerobic growth of heterotrophs, expressed as process j = 1 on the ASM1 matrix (see 2.2.2); the aerobic growth of autotrophs, which is process j = 3; the decay of heterotrophs, process j = 4; and the decay of autotrophs, denoted as j = 5. Otherwise, the components involved are the readily biodegradable substrate (SS ), active heterotrophic biomass (XB,H ), active autotrophic biomass (XB,A ), oxygen (SO ), nitrate and nitrite nitrogen (SN O ) and N H4+ + N H3 nitrogen (SN H ). These components are i = 2, 5, 6, 8, 9 and 10 respectively on the ASM1 notation.

Mathematical modeling of heterogeneous biological wastewater processes

14

2.3

Mass balances

Commonly, there are two patterns of flow considered in wastewater treatment: complete mixed and plug flow. If the hydraulic regime is complete mixed, it means that uniform conditions in terms of concentration of the components are given in the whole volume of water. Otherwise, the plug flow conditions represents the need of considering an elemental volume because each one has different characteristics. Plug flow

Complete mixed

Figure 2.4: Schemes of both hydraulic regime: plug flow and complete mixed Traditionally, mass balances have been used on the study of the components concentration in water. A mass balance can be represented as IN − OU T + GEN ERAT ION = ACCU M U LAT ION

(2.19)

where IN - OUT represents the net transport which is given through the volume, while GENERATION expresses the production or destruction given within it and ACCUMULATION is the amount left over.[5] Wastewater treatment wetlands are normally constructed as a large basin where the direction of the flow is in one direction and the concentration of its components changes along it, so it has not an homogeneous composition. It is a common practice on the study of these wetlands to divide the domain representing the wetland into different volumes, where water flows through all of them and each one is considered as a complete mixed basin. But due to the characteristics explained above, it seems more realistic to study its hydraulics as a plug flow problem. Then, if area and volumetric flow are assumed

State of the art

15

as constant, and flow is considered to be only along one direction, the mass balance will be µ ¶ ∂C ∂C QC − Q C + dx + rdV = dV (2.20) ∂x ∂t where Q is the volumetric flow rate, C is the concentration of the constituent, dV is the differential volume element and r is the reaction rate. Dividing the above equation for dV and assuming dV as Adx −

Q ∂C ∂C +r = A ∂x ∂t

(2.21)

So the behavior of the process can be expressed as a first order hyperbolic PDE. Effect of diffusion may be included through second derivative with respect to x. Then it will be written as −

2.4

Q ∂C ∂ 2C ∂C +ν 2 +r = A ∂x ∂x ∂t

(2.22)

Adsorption

Adsorption is a process that occurs when a gas or a liquid solute accumulates on the surface of a solid forming an atomic or molecular film. In the interior of the solid, all the atomic bonds are filled, while on the surface, these experiment a discontinuity. For those incomplete bonds it is energetically favorable to react with other atoms or molecules found on the surrounding medium. When adsorption is considered, the problem may change. Lets see its effects on the mass balance done before (equation 2.21) φ

∂ S˜ Q ∂C ∂C + (1 − φ) =− + r(c) ∂t ∂t A ∂x

(2.23)

where, while C is the concentration of solute, S˜ is the concentration adsorbed, and φ is the medium’s porosity. Two cases of adsorption can be given, one called equilibrium, which occurs when chemical reactions are fast, and the other one named kinetic, that occurs when they

Mathematical modeling of heterogeneous biological wastewater processes

16

are slow. In the first case, adsorption is considered as a function of the concentration in the fluid, S˜ = a(c) while when kinetic adsorption is given, it is considered that it satisfies S˜ ˜ = k(a(c) − S) dt where k is a rate parameter. From now on, only the equilibrium case is going to be considered. Different assumptions can be made to consider the function a, which receives the name of adsorption isotherm.2 One of them is the Langmuir isotherm (I. Langmiur, 1916), which is described as KC a(C) = N K>0 (2.24) 1 + KC where N is the saturation concentration of the adsorbed solute and K is an empirical constant related to the energy. When adsorption is achieved but on little quantities, the expression 2.24 can be simplified so that adsorption can be represented as linear (Henry’s law) a(C) = N KC

(2.25)

Under hypothesis of equation 2.25, equation 2.23 may be rewritten as [φ + (1 − φ)N K]

∂C Q ∂C =− + r(c) ∂t A ∂x

(2.26)

while K does not depend neither on time nor on concentration. If H = φ + (1 − φ)N K is defined, it can be seen that H

∂C Q ∂C =− + r(c) ∂t A ∂x

(2.27)

so, when adsorption is considered, the equation 2.21 is modified only by a term , H, multiplying the variation of concentration along the time. While it is a linear modification, which means not important modifications on the procedure of solving the problem, from now on, the adsorption problem has been neglected on implementing the solution. 2 It is called adsorption isotherm because the adsorption is studied considering constant temperature conditions.

Chapter 3 Mathematical model As it was exposed on section 2.3, the traditional way on studying the biochemical models on wastewater treatment is through mass balances, and there a particular and simplified case was developed. In this chapter, a more general point of view on mass balances is going to be developed, so that the equation that reigns the transport of particles studied can be found. Then, the reaction model will be mathematically presented. Moreover, the boundary conditions of the problem are going to be defined. Along this chapter, all the aspects to reach the complete problem statement will be defined, so that it can be presented at the end of it.

3.1

Continuum mechanical formulation

The goal of this work is implementing the wastewater treatment model from the continuum, considering wastewater as a continuous matter. Under this concept it is understood a infinite set of particles which is going to be studied macroscopically, without considering the possible discontinuities given at the microscopical level. That means that the mathematical description of this medium and its properties can be described by continuous functions [15]. At each point of this continuous matter there is supposed to be a unique value of the velocity, pressure, density and other so-called field variables. The continuous matter is then required to obey the conservation laws of mass, momentum and energy, which give rise to a set of differential equations governing the field variables.

Mathematical modeling of heterogeneous biological wastewater processes

18

Lets see some aspects that may help on defining the mathematical model.1

3.1.1

Eulerian and Lagrangian coordinates

As the continuum approach has been chosen, it should be then considered which reference frame may be used for the development of the formulation. There are basically two coordinate systems that might be employed, eulerian and lagrangian coordinates. On eulerian coordinates, also known as spatial coordinates, attention is focused on the matter which passes through a control volume that is fixed in space. So the matter inside the control volume at any instant in time will consist in different particles than the given at any other instant of time. In this case the independent variables are the spatial coordinates x, y and z and time t. Otherwise, in lagrangian or material approach, attention is fixed on a particular mass of the matter as it moves. Now a particle will be followed while it moves and changes the shape through the continuous matter. In this reference frame, independent variables are the coordinates x0 , y0 and z0 which the specific particle passed through at the instant of time t0 [3]. So the independent variables are x0 , y0 and z0 as also t.2 On the resolution of the problem, an eulerian coordinate system has been employed, so the mesh used on the finite element discretization will remain constant on space and along time.

3.1.2

Material derivative

Let φ be any field variable. From the eulerian point of view, φ may be considered as a function of x, y, z and t. But if the lagrangian approach is made, the specific particle observed during a period of time δt, is going to change its position an amount of δx, δy and δz while its value of φ will change δφ. During the δt the change of φ will be δφ =

∂φ ∂φ ∂φ ∂φ δt + δx + δy + δz ∂t ∂x ∂y ∂z

(3.1)

Equating this change in the lagrangian framework and dividing throughout δt 1 More information about the constitutive laws and other aspects of continuum mechanics may be found, between others, in [3] and [15]. 2 Note that some bibliography represent the independent variables on the lagrangian reference frame as X, Y , Z and t.

Mathematical model

19

gives

∂φ δx ∂φ δy ∂φ δx ∂φ δφ = + + + (3.2) δt ∂t δt ∂x δt ∂y δt ∂z The left-hand side of the expression represents the total change of φ as observed in the lagrangian framework during δt. In the limit, reducing δt to zero, the lefthand side represents the time derivative of φ, while the ratios δx/δt, δy/δt and δz/δt express the velocity, u, in each component. So it can be written as Dφ ∂φ = + (u · ∇) φ Dt ∂t

(3.3)

which is known as the material derivative. It might be remarked that the term (u · ∇) φ is known as the convective derivative.

3.1.3

Reynolds transport theorem

Let A being a property of the matter treated in the continuum, and φ(x, t) the quantity of this property per mass unit. Then, ρφ(x, t) is the quantity per volume unit. Consider then an arbitrary material volume within the continuum which in the instant t occupies a volume Vt ≡ V . The quantity of A given in this material volume Vt at instant t will be: Z Q(t) =

ρφdV

(3.4)

Vt ≡V

Then, the variation of this property A in the material volume Vt , that in t coincides with the volume V , will be given by the material derivation Z Z Z d ∂(ρφ) ρφdV = dV + ∇ · (ρφu)dV (3.5) dt Vt ≡V ∂t V V Using the Reynolds Lema 3 on the left side of the equation and the divergence theorem4 on the right one, it can be seen that Z Z Z dφ ∂(ρφ) ρ dV = dV + ρφu · ndS (3.6) dt ∂t V V ∂V where n is a unitary vector normal to the volume surface. It may be rewritten as Z Z Z dφ ∂ (ρφ)dV = ρ dV − ρφu · ndS (3.7) ∂t V dt V ∂V d The Reynolds Lema may be expressed as dt 4 The divergence theorem can be announced as

3

R

R ρφdV = V ρRdφ dt dV ∇ · ΦdV = Φ · ndS V ∂V

Vt ≡V R

Mathematical modeling of heterogeneous biological wastewater processes

20

The left-hand side of the expression can be seen as the variation of the quantity of property A in the control volume per time unit, while on the right-hand side, the first term expresses the variation of the quantity of A due to the internal variation and the second term expresses the variation given due to the convective flux of A through the surface ∂V . The equation 3.7 is commonly known as the Reynolds transport theorem. It can also be expressed on the local form as ∂(ρφ) dφ =ρ − ∇ · (ρφu) ∂t dt

∀x ∈ V

(3.8)

As it was said, ρφ(x, t) is the quantity per volume unit, that is the concentration of A. While working on an incompressible medium (and in this case, it is so), ρ may be considered constant. Otherwise, under this conditions, it might be considered that ∇·u=0 (3.9) So changing ρφ(x, t) for concentration (let it be c), and considering incompressible flow hypothesis, the equation will rest as ∂c dc = − u · ∇c ∂t dt

3.1.4

(3.10)

Diffusion

Diffusion is a transport procedure caused by the Brownian movement. The Brownian movement is caused by the crash of the molecules of the fluid with the particles studied. To represent this phenomena there are the Fick’s laws (Adolf Fick, 1885). In this case, second Fick’s law ought to be used, it describes diffusion on the threedimensional case as ∂c = ∇ · (ν∇c) (3.11) ∂t It may be remarked, that ν is a symmetric positive defined matrix, but when the medium may be considered as isotropous, then ν is a scalar.

3.1.5

Convection-diffusion-reaction

As diffusion processes are involved on particles fluid transport, it is necessary to add this term at the Reynolds transport theorem. So, from equations 3.10 and 3.11,

Mathematical model

21

transport of particles on a fluid media can be written as ∂c dc = ∇ · (ν∇c) − u · ∇c + ∂t dt

(3.12)

where on the right-hand side of the equation, the first term involves transport by the diffusion process, the second term expresses the convective contribution on the transport and the last term involves the reaction (or as it was said before, the production of a specie within the control volume). As it was presented on section 2.2.2, the reaction presented here as dc/dt is, in fact, the observed conversion rate defined by the ASM1. To characterize the relative importance between convection and diffusion there is the P´eclet number. It can be defined as Pe =

uh 2ν

(3.13)

which expresses the ratio of convective to diffusive transport. As seen before, u expresses convective velocity, while ν is the diffusion coefficient. And h expresses the length of transport.

3.2

Reaction model

As it was exposed on 2.2.2, the ASM1 model lets express the conversion rate of a component as the sum of each stoichiometric coefficient, related to that component, multiplied by the process rate to which it belongs. X ∂c = ∇ · (ν∇c) − u · ∇c + νi,j ρj ∂t j

(3.14)

where, as also was said before, j represents each process that has been described on the Activated Sludge M odel N o.1, while i denotes each compound involved on those processes.

3.3

Boundary conditions

To solve a differential equation, it is necessary to know the boundary conditions that define the problem. The most common types of boundary conditions are Dirichlet,

Mathematical modeling of heterogeneous biological wastewater processes

22 Neumann and Robin’s ones.

Let Ω be the domain of the problem. In order to define boundary conditions along ∂Ω it is necessary to equal the flux on both sides of a boundary. It can be considered that out of the domain, there will be only convective flux, while on the inside, convection and diffusion will take part on the transport process (reaction term may be neglected while defining the boundary conditions). So, if u− is the velocity and c− is the concentration on the outside of the boundary, and u+ and c+ are the velocity and the concentration on the inside, the expression of the equation will be, c− u− · n = c+ u+ · n − n · ν∇c (3.15) that it may be rewritten as ¡ ¢ n · ν∇c = c− u− − c+ u+ · n

(3.16)

This is the most general form on writing a boundary condition on a convectiondiffusion problem, which corresponds to the Robin boundary conditions type. On entrance boundary, it may be supposed that velocity on both sides of the boundary is the same (u− = u+ = u). Then equation 3.16 is expressed as ¡ ¢ n · ν∇c = c− − c+ u · n (3.17) Lets consider that diffusion is much more little than convection in the inside of the domain, then ¡ ¢ 0 = c− − c+ u · n (3.18) so that c− = c+

(3.19)

This hypothesi can be made while diffusion coefficient, ν, is really small. Although the value of n · ∇c is not known a priori, its importance on the entrance boundary condition has not such an influence [17] [18]. It can be seen that the condition along the entrance boundary, is a Dirichlet condition. On both lateral boundaries transport velocity is parallel to the boundary, so that the scalar product u · n will be equal zero. That leads to n · ν∇c = 0

(3.20)

Otherwise, the exit boundary has on both sides the same concentration as also the same velocity, c − u− = c + u+ (3.21)

Mathematical model

23

which entails that expression 3.16 can be rewritten as n · k∇u = 0

(3.22)

It may be seen that both, lateral and exit boundaries, can be written as Neumann conditions. Those conditions normally represent a impermeable boundary, but on convection-diffusion problems they have another meaning, while they can define exit or lateral boundaries. It is important to remark that oxygen should enter through a lateral boundary, but as it has been explained on section 2.2.2, oxygen will be considered with an expression added on its equation. Although it seems that this would not reproduce a realistic entrance of oxygen, it has the advantage that on subsurface flow wetlands oxygen does not only enter through a boundary, but also through the plants living in there. So that there is also oxygen entering from insider points of the domain.

3.4

Problem statement

After going through the sections of this chapter, the whole definition of the problem may be presented. Let Ω = [0, 1]x[0, 1] be the spatial domain associated to the problem, while the time domain is ]0, 2[. Then, the system of equations may be presented as ∂c = ∇ · (ν∇c) − u · ∇c + r(c) ∂t

in Ω x ]0, T [

(3.23)

As it was explained before, the reaction term, r(c) has two different forms. The first one involves six components from the whole formulation of the ASM1, so the system will have six equations. The unknowns of the system may be then SS , XB,H , XB,A , SO , SN O and SN H (see section 2.2.2) so the reaction term may be presented as 8 X νi,j ρj i=2,5,6,8,9 and 10. (3.24) j=1

The other formulation used as reaction term involves the whole ASM1 model. So thirteen equations may be given and the unknowns are all the components expressed on the first array of the ASM1 matrix. Then, the reaction term can be written as 8 X j=1

νi,j ρj

i = 1, ..., 13.

(3.25)

24

Mathematical modeling of heterogeneous biological wastewater processes

Attending to the spatial domain defined, boundary conditions may now be written as c(x, t) = cext on x = 0 cn (x, t) = 0 on x = 1 ∪ y = 0 ∪ y = 1

(3.26)

where x = 0 represents the entrance boundary, while y = 0 and y = 1 designate the lateral ones and x = 1 is the exit boundary. When the problem is solved on one dimension, as the domain of the problem is Ω = [0, 1], boundary conditions may be rewritten as c(x, t) = cext on x = 0 cx (x, t) = 0 on x = 1

(3.27)

Last but not least, initial condition may be also presented. It can be expressed as c(x, 0) = c0 (x)

on Ω

Chapter 4 Numerical solver In order to solve this unsteady convection-diffusion-reaction problem, it is going to be used a double discretization scheme involving linear C 0 finite elements in the spatial discretization and a multistage step method which involves only first derivatives in the time discretization. The finite element methods is nowadays a basic tool on the work related to civil engineering. The strongpoints of this methods are the consistent treatment of the differential boundary conditions, the easy modeling of complex continuum mechanical problems and domains as also the options given to be programmed in a flexible format[4]. A multistage step method is a integration method where every time-step, is subdivided into different stages in order to find the value reached after the time-step. The most traditional multistage methods are the Runge-Kutta methods and the Pad´e approximations. At the end of this chapter, stabilization techniques for the spatial discretization are going to be presented in order to avoid the problems associated to the convection term.

4.1

Spatial discretization

For the spatial discretization the Galerkin formulation of the method of weighted residuals is going to be used. When it is applied to problems governed by self-adjoint elliptic or parabolic differential equations, this method leads to symmetric stiffness matrices.

Mathematical modeling of heterogeneous biological wastewater processes

26

But on flow problems this advantage is not presented while the convection operator is not symmetric. This may lead to oscillations on the solution, but this will be studied on section 4.3. Lets see how is the Galerkin procedure implemented on the problem. Although it is a system of equations, its development will be made on an equation without the lose of veracity. All the processes explained on this section and the later ones may be then extended to all the equations of the problem. Rewriting equation ?? c˙ + u · ∇c − ∇ · (ν∇c) = r(c)

(4.1)

and using the differential operator L = u · ∇ − ∇ · (ν∇)

(4.2)

the problem can be rewritten as c˙ + Lc = r(c)

(4.3)

Let V be the space of weighting functions which satisfy the homogeneous boundary conditions on ΓD . V = {w(x) ∈ H1 (Ω) | w(x) = 0 f or x ∈ ΓD } It can be noticed that these w functions are not time-dependant. The solution c of the system lies in L (0, T ; H1 (Ω)). And the time dependency of the approximate solution c can be translated to the trial space St , which does vary on time, St := {c | c(·, t) ∈ H1 (Ω),

t ∈ [0, T ] and c(x, t) = cD f or x ∈ ΓD }

So the weak form of the initial boundary value problem is defined as (w, c) ˙ + (w, Lc) = (w, r(c)) where

(4.4)

Z (w, Lc) =

(w (u · ∇c) − ∇w · (ν∇c)) dΩ, Z (w, r(c)) = wr(c)dΩ

(4.5)





(4.6)

Numerical solver

27

The spatial discretization of the problem on attendance to the Galerkin formulation consists on defining two dimensional spaces Sth and V h , subsets of St and V, V h := {w ∈ H1 (Ω) | w|Ωe ∈ P 1 (Ωe ) ∀e and w = 0 on ΓD } (4.7) Sth := {c | c(·, t) ∈ H1 (Ω), c(·, t)|Ωe ∈ P 1 (Ωe )t ∈ [0, T ] ∀e and c = cD on ΓD } (4.8) 1 where P is the finite element interpolating space formed by polynomials of order 1. The semi-discrete Galerkin formulation is obtained by restricting the weak form to these finite dimensional spaces. So for any t ∈ [0, T ] find ch ∈ Sth such that for all wh ∈ V h , (wh , c˙h ) + (wh , Lch ) = (wh , r(ch )) (4.9)

The approximation ch can be written as X X NA (x)cD (xA , t), NA (x)cA (t) + ch (x, t) = A∈ηD

A∈η\ηD

where η is the set of global node numbers in the finite element mesh and ηD ⊂ η the subset of nodes belonging to the Dirichlet portion of the boundary. The test functions are defined as wh ∈ V h = spanA∈η\ηD {NA }. It shall be noticed that on ch (x, t), the shape functions NA (x) do not depend on time, but time dependence is accounted by the nodal values of the unknown. The assembly process delivers finally the semi-discrete system M c˙ + (C − K)c = f

(4.10)

where M , C and K are the consistent mass matrix, the convection matrix and the diffusion matrix respectively. These matrices are obtained by topological assembly of element contributions R e M = Ae M e Mab = Ωe Na Nb dΩ C = Ae C e K = Ae K e

e Cab = e Kab =

R Ωe

R Ωe

Na (u · ∇Nb )dΩ

∇Na · (ν∇Nb )dΩ

where A is the assembly operator, 1 ≤ a, b ≤ nen and nen is the number of element nodes.

Mathematical modeling of heterogeneous biological wastewater processes

28

Otherwise, the vector f is built up also with the help of the assembly operator: R f = Ae f e fae = Ωe Na r(c))dΩ It may be noticed that equation 4.10 can be rewritten as ¡ ¢ c˙ = M−1 (C − K) c + M−1 f = B · c + f¯(c)

(4.11)

Remark 4.1.1 At each time-step, equation 4.11 must be solved. It may be noticed that on this equation, the only time-dependant term is the vector f¯ -it depends on time due to the reaction term, which depends on the existing concentration at each point of Gauss, and this concentration changes along the time-, while B does not change in time. While f must be build up at each time step, matrices M , C and K should be calculated out of the time integration process, in order to reach an implementation that may be executed with less waste of time. Remark 4.1.2 As it has been seen before, the reaction term depends on the concentration, so due to the discretization it can be expressed as ³X ´ X r(c(x, t)) = r cA (t)NA (x) ' r(cA (t))NA (x) Then, while building up f , it can be seen that Z Z X X e e fa = Na r(c)dΩ = r(cA (t)) Na NA dΩ = r(cA ) · MaA Ωe

Ωe

So that now f are the reaction values evaluated on the nodes multiplied by the mass matrix M . This has been used on the one-dimension examples. Then equation 4.11 may be written as ¡ ¢ c˙ = M−1 (C − K) c + M−1 Mr(c) = B · c + r(c) (4.12)

4.2

Time discretization

Unsteady convection-diffusion-reaction problems are not easy to be solved by a highorder time-stepping method because of the second-order diffusion operator given in

Numerical solver

29

the equation. When C 0 finite elements are implemented, the time-stepping schema should only involve first order derivatives[4]. This leads, between others, to the Runge-Kutta methods. Runge-Kutta methods are multistage methods that only use the solution cn at time tn to compute the solution cn+1 at time tn+1 . This can be made by calculating intermediate values of the time derivative of the unknown c within the interval ∆t = tn+1 − tn . Equation 4.11 shall be rewritten as c˙ = F(t, c)

(4.13)

Then, time integration shall be applied in order to solve this ordinary differential equation. So integrating respect t on both sides of the equation, Z

tn+1 tn

∂c dt = ∂t

Z

tn+1

F(t, c(t))dt

(4.14)

tn

While the left-hand side of the equation may be represented as Z c

n+1

tn+1

n

−c =

F(t, c(t))dt

(4.15)

tn

the right-hand side shall be integrated using a numerical quadrature, where the integration points are expressed as ξ i , while the weights are denoted as ∆tbi . Then c

n+1

n

− c = ∆t

ntg X

bi F(ξ i , c(ξ i )) + E(∆t).

(4.16)

i=1

While to find out c(ξ i ) it can be seen Z i

n

ξi

c(ξ ) − c = tn

∂c dt = ∂t

Z

ξi

F(t, c(t))dt

(4.17)

tn

so it shall be used a numerical quadrature again, on the same integration points but with different weights, which may help to integer between [ti , ξ i ]. This leads to i

c(ξ ) − ci = ∆t

ntg X j=1

ai,j F(ξ j , c(ξ i )) + Ei (∆t).

(4.18)

Mathematical modeling of heterogeneous biological wastewater processes

30

Despising the errors E(∆t) and Ei (∆t) on both equations 4.16 and 4.18, changing the variable ξ i = tn + c˜i ∆t and using the notation li = F(ξ i , c(ξ i )), the ntg -stage Runge-Kutta method can be expressed as c

n+1

n

= c + ∆t

ntg X

bi li ,

(4.19)

i=1

where li is defined as n

n

li = F(t + c˜i ∆t, c + ∆t

ntg X

ai,j lj ),

i = 1, 2, ..., ntg .

(4.20)

j=1

The coefficients ai,j , bi and c˜i expressed on the former equations must verify the consistency conditions c˜i =

ntg X

ai,j

and

j=1

ntg X

bi = 1

i=1

and they are normally displayed on a table called the Butcher array, which is built up as Table 4.1: Generic Butcher array used on Runge-Kutta methods. c˜1 c˜2 .. .

a1,1 a2,1 .. .

a1,2 a2,2 .. .

... ... ...

a1,ntg a2,ntg ...

c˜ntg

antg ,1 b1

antg ,2 b2

... ...

antg ,ntg bntg

When ai,j = 0 ∀j ≥ i, a explicit Runge-Kutta method is presented, otherwise it would be implicit. A way on improving the resolution of the ODE is using an algorithm with a step size control, which means that on each step, the length of it, ∆t, may be different depending on the easiness of the function that may be integrated. One of this methods is the known as Runge-Kutta-Fehlberg 4-5 (E. Fehlberg, 1969). It defines the step-size that shall be used at each step on comparing the solution given by a forth-order Runge-Kutta method with the one found with a

Numerical solver

31

fifth-order one. And its development was planned on trying to reduce as much as possible the truncation error that may present while the integration step-size depends strongly on the magnitude of these errors[10]. It can also be presented on a Butcher array, but as it can be seen next, it has some changes in the structure due to the representation of two Runge-Kutta methods with different order. Table 4.2: Butcher array for the Runge-Kutta-Fehlberg 4-5. 1 4

1 4

3 8

3 32

9 32

12 13

1932 2197

−7200 2197

7296 2197

1

439 216

−8

3680 513

−845 4104

1 2

−8 27

2

−3544 2565

1859 4104

−11 40

25 216

0

1408 2565

2197 4101

−1 5

0

16 135

0

6656 12825

28561 56430

−9 50

2 55

1 360

0

−128 4275

−2197 4275

1 50

2 55

The sixth array is composed by the coefficients bi used to solve the ODE by a fourth order method, while the coefficients in the seventh array are used for the fifth order method. Finally, the eighth array is used to calculate the error committed (note that the values are those values of the seventh array minus the ones on the sixth). So the error committed by the use of a ∆t, E∆t , may be found as · ¸ 1 128 2197 1 2 E∆t ' ∆t k1 − k3 − k4 + k5 + k6 360 4275 75240 50 55 From the expression above, the step size can be found, when a tolerance is given. Limiting the error to that tolerance and solving the part in brackets, the ∆t is defined and the Runge-Kutta method of fourth order can be solved.

Mathematical modeling of heterogeneous biological wastewater processes

32

4.3

Stabilization

The convective term may bring instabilities on the space-discretization. In order to stabilize this term in a consistent manner, ensuring that the solution of the differential equation is also a solution of the weak form, there are different techniques. The most commonly used are: - The Streamline Upwind Petrov-Galerkin method (SUPG), - Galerkin/Least-squares (GLS), - Subgrid scale method (SGS), and - Least-squares. The three first ones have a similar structure. It consists on adding an extra term over the element interiors to the Galerkin weak form, which could be seen as a diffusion addition which ensures consistency. This term is a function of the residual of the equation. The residual of the equation 4.3 is R(c) = c˙ + Lc − r(c)

(4.21)

The general form of these stabilization techniques is (w, c) ˙ + (w, Lc) +

nel X

(P(w), τ R(c))Ωe = (w, r)

(4.22)

e=1

where P(w) is a certain operator applied to the test function, τ is the stabilization parameter and R(c) is the residual of the differential equation (equation 4.21). It may be noticed that R(c) is computed only for each element interior Ωe Those three stabilization techniques differ ones from the others by the definition of the operator P(w). After taking a glance over all of them, it shall be seen that, because of the use of C 0 elements on the Galerkin discretization, the three techniques apply the same operator P(w). Lets start the Streamline Upwind Petrov-Galerkin stabilization method, which may be characterized by P(w) = u · ∇w, (4.23) The restriction imposed on the weak form to the usual finite dimensional subspaces, defines the discrete problem that must be solved, ch ∈ S h has to be found such that (wh , c˙h ) + (wh , Lch ) + ¢ ¤ ¢ £ ¡ P R ¡ + e Ωe u · ∇wh τ c˙h + u · ∇ch − ∇ · ν∇ch − r dΩ = (wh , r) (4.24)

Numerical solver

33

Otherwise, stabilization may be also achieved using the Galerkin/Least-squares technique as it is shown next. On the use of GLS, the operator is defined as P(w) = L(w) = u · ∇w − ∇ · (ν∇w)

(4.25)

For an unsteady-convection-reaction problem, the weak form restricted to the standard finite dimensional subspaces may be written as (wh , c˙h ) + (wh , Lch ) + ¡ ¢¤ £ ¡ ¢ ¤ P R £ + e Ωe u · ∇wh − ∇ · ν∇wh τ c˙h + u · ∇ch − ∇ · ν∇ch − r dΩ = = (wh , r) (4.26) Another stabilization technique is the subgrid scale method (SGS), first introduced by T.J.R. Hughes (1995). Here the P(w) is quite similar than the one seen on GLS, but instead of using the operator L it shall be used its adjoint. That means P(w) = L(w) = u · ∇w + ∇ · (ν∇w)

(4.27)

Then, restricting the weak form to the standard dimensional subspaces the problem can be expressed as (wh , c˙h ) + (wh , Lch ) + ¡ ¢¤ £ ¡ ¢ ¤ P R £ + e Ωe u · ∇wh + ∇ · ν∇wh τ c˙h + u · ∇ch − ∇ · ν∇ch − r dΩ = = (wh , r) (4.28) But, as it has been previously announced, some changes may occur on the definition of those operators due to the use of linear elements. The expression for P(w) may be simplified on GLS and SGS methods while the second-order derivatives on linear elements are zero in the element interiors. So the operator P(w) becomes P(w) = u · ∇w

(4.29)

and the weak form restricted to the usual dimensional subspaces may be presented as (wh , c˙h ) + (wh , Lch ) + ¤ £ ¤ P R £ + e Ωe u · ∇wh τ c˙h + u · ∇ch − r dΩ = (wh , r)

(4.30)

So, as it can be seen, the three methods are the same on the resolution of the problem with the C 0 elements that have been chosen for the Galerkin formulation. Another method used for the stabilization is the least − squares (LS) method. The use of Galerkin finite element method is not a good procedure for a spatial

Mathematical modeling of heterogeneous biological wastewater processes

34

discretization on convection dominated problems [12]. The LS procedure does not add a stabilization term on the weak form, but uses another discretization procedure. It is remarkable that the classical procedure of the least squares formulation for convection-diffusion problems requires a work in H2 , unless a mixed LS formulation is used, while the minimization of the square of the residual of the governing equation that uses LS includes second spatial derivative. But there are multiple examples of those mixed LS formulation which lead to the use of H1 [4], [11] . Those mixed LS applied on transport problems generate a symmetric positive system, in contrast to the Galerkin method where the system is non-symmetric. The LS equation can be written as ¡ ¢ L(w), R(cn+1 ) = 0 (4.31) where R(c) is, as defined before, the residual of the equation 4.3. As it can be noticed on equation 4.31, on c time discretization has already been applied. That is because on the use of a LS procedure, firstly the time discretization must be done, and later the spatial one can be applied. Traditionally the mixed LS have been used with the known as Pad´e approximations for the time discretization. These methods are a multistage factorization of the Pad´e approximation to the exponential function and can be expressed in the incremental form as cn+βi

cn+β1 = cn n+β = cn + βi 4tct i−1 , i = 2, ..., ntg + 1

(4.32)

where βi can also be expressed by a Butcher array. But the whole implementation using Pad´e time-discretization and the later least squares spatial-discretization is laborious enough to decide to implement the Galerkin discretization with the Runge-Kutta-Fehlberg method introducing a parameter to reach stabilization.

Chapter 5 Validation Once the numerical solver is implemented, it is necessary to validate it. Although the mathematical model is not such a realistic model, the tendency of the solution may allow to see if the steps done before were correct or not. Firstly, the parameters chosen for the resolution will be presented. Those refer to some aspects of the solver as also of the mathematical model. Then, the model will be validated in one dimension. This means playing a bit with the formulation, restricting some components so that the difference on the results may help on establishing its goodness. Once this has been made, it is necessary to verify also the two dimensions one. For those models where diffusion is added in order to obtain a stabilized solution, a little study of the stabilization parameter and its order of magnitude will be also presented.

5.1

Parameters

Before solving the problem, there are some parameters on the equation which should be specified. There are two types of parameters, the ones needed exclusively for the implementation and the others which define the mathematical model. Those which refer to the numerical solver are the number of elements on the mesh and the tolerance admitted on the solution. The mesh defined to solve the problem is formed by 20 elements on the one dimensional case, while on the two dimensional one, the domain Ω is defined by 20x20 elements. In reference to tolerance admitted,

36

Mathematical modeling of heterogeneous biological wastewater processes

it has been chosen a tolerance of 0.5x10−4 . When looking to the mathematical model, the parameters that must be defined are the ones refereing to the constants on transport equation, those which appear on the ASM1 formulation, as also the initial and boundary conditions values. Constants which appear on the ASM1 may be presented on two different groups: the ones refereing to stoichiometric processes and the others related to the kinetics. On table 5.1 the first ones, which are given on the coefficients presented within the matrix, can be seen.

Table 5.1: Stoichiometric parameters appeared on the ASM1. Source:J. Bolmstedt and G. Olsson, 2002. Parameter Value Units Definition YH 0.67 g cell COD formed/(g yield of growth rate for hetCOD oxidized) erotrophic biomass YA 0.24 g cell COD formed/(g yield of growth rate for auCOD oxidized) totrophic biomass fp 0.08 fraction of biomass yielding particulate products iXB 0.086 g N/g COD in mass N/mass COD in biomass biomass iXP 0.06 g N/g COD in endoge- mass N/mass COD in products nous mass biomass The stoichiometric parameters, which refer to the ones presented on the process rates, are expressed on table 5.2. Values presented on both tables are realistic values that might be used on a wastewater treatment plant formed by different reactors when its operation is made under constant temperature of 20◦ C [1]. Although this has nothing to do with the initial motivation of the problem, its orders of magnitude seem to be similar to the real ones, and this may help on detecting possible troubles which could appear on the real problem. Then, on convection-diffusion equation, as diffusion coefficient, ν, the value ν = 10−4 has been employed (as it has been previously said, diffusion coefficient is an scalar when isotropy conditions are considered, as it is in this case). Otherwise convection velocity has been defined as u = (1, 0), so that flux is only in one direc-

Validation

37

Table 5.2: Kinetic parameters appeared on the ASM1. Source:J. Bolmstedt and G. Olsson, 2002. Parameter Value Units Definition −1 µ ˆH 6. h maximum specific growth rate for heterotrophic biomass 3 KS 20. g of COD/m saturation coefficient for heterotrophic biomass g of 02 /m3 oxygen saturation coefficient for KO,H 0.2 heterotrophic biomass 3 KN O 0.5 g N O3 − N/m nitrate hsc for denitrifying heterotrophs −1 bH 0.62 h decay rate for heterotrophic biomass −1 µ ˆA 0.8 h maximum specific growth rate for heterotrophic biomass KN H 1. g of N3 − N/m3 ammonium saturation coefficient for autotrophic biomass KO,A 0.4 g of 02 /m3 oxygen saturation coefficient for autotrophic biomass −1 h decay rate for autotrophic bA 0.2 biomass ηg 0.8 correction factor for anoxic growth of heterotrophs 3 ka 0.08 m / g COD day ammonification rate kh 3. g slowly biodeg.COD max. specific hydrolysis rate / g cell COD day KX 0.03 g slowly biodeg.COD hsc for hydrolysis of slowly / g cell COD biodeg. substrate ηh 0.4 correction factor for anoxic hydrolysis

38

Mathematical modeling of heterogeneous biological wastewater processes

tion. When refereing to the Dirichlet boundary condition, the value that has been chose is cext = 1 for all the different unknowns. Otherwise, as initial condition value, c0 , two different values have been chosen on the implementation. Those are c0 = 0 and c0 = 1, and they will be specified on each case when results are shown.

5.2

Validation of the system resolution in one dimension

After defining parameters, the system of equations may be solved. In this section, results taken by the model will be analyzed to validate its goodness. Equations have been solved in one dimension, without using any stabilization technique and its results have been found out with initial condition c0 = 1. The procedure employed in order to evaluate the model will consist on comparing the results given in different situations of components presence. Assuming that the main difference can be given by the presence of oxygen, the different casus will be seen on both conditions, when entrance of oxygen is given through the whole domain and when this entrance is not allowed. As it has been previously explained on section 2.2.2, entrance of oxygen means adding the term on equation 2.17 on the reaction term in oxygen’s equation (SO ). So, over this two different situations, the presence or absence of other influent components will be forced in order to observe how results vary.

5.2.1

Validation of the six equations modeling

Lets start by the analysis of the six equations model. As a reminder, this model analyzes six components (SS , XB,H , XB,A , SO , SN O and SN H ) which take part on four different processes (aerobic growth of heterotrophs (ρ1 ), aerobic growth of autotrophs (ρ3 ), decay of heterotrophs (ρ4 ) and decay of autotrophs (ρ6 )). On the six equations model the component which will be restricted is the N H4+ and N H3 nitrogen, (SN H ). This has been made so because, when other components are restricted, not such big differences can be observed.

Validation

39

No oxygen’s entrance 1.5

1.5

1

1

1

1

Ss 0.5

0

0.5 xnodo

0

1

0.5

0

0.5 xnodo

0

1

1.5

1

1

1

0.5

0

0.5 xnodo

0

1

2

0.5 xnodo

0

1

Snh 0

0.5 xnodo

1

0

0

0.5 xnodo

0.5 xnodo

1

0

0.5 xnodo

1

0

0.5 xnodo

1

5

0

1

1.5

1.5

0.5

0.5

0

10

2

1

1

0

0

1.5

1.5

0

1

0.5

Sno

0

0.5 xnodo

So

So 0.5

0

Xba

1.5

Xba

1.5

0.5

1 Snh

Ss 0

Xbh

1.5

Xbh

1.5

0.5

Sno

Oxygen’s entrance

1

0.5

0.5 0

0.5 xnodo

1

0

0

0.5 xnodo

1

0

Figure 5.1: Evolution of concentration on the six equation model without and with oxygen entrance

Mathematical modeling of heterogeneous biological wastewater processes

40

First of all, lets observe results with and without entrance of oxygen maintaining other components not limited (they may be seen on figure 5.1). There concentration along the domain Ω = [0, 1] is presented when solution reaches the steady state. This may be seen as the concentration of each component along the subsurface flow wetland when their stationary flow is achieved. The main differences between both situations are given by SN O and SN H . When presence of oxygen is assured, production of SN O is given while in its absence (on the second half of the domain on the N o oxygen0 s entrance’s graphics), there is neither production nor destruction. That occurs so because the process rate on its reaction term goes to zero when absence of SO is given. A similar situation is given on SN H destruction, as when SO is not present on the medium, its decreasing is stopped, this time because both processes rates, ρ1 and ρ3 , become zero. From the biological point of view, as denitrification processes are not considered on this model (only aerobic processes are taken under consideration), destruction of N H4+ can not occur without oxygen presence, as it may be seen on nitrification equilibrium on section 2.1.2. Therefore nitrate and nitrite can neither be produced. On the other results, similar behavior may be observed, but its effects are not that spectacular as on the ones above commented. Results produced when N H4+ and N H3 nitrogen reaction term is forced to be zero can be seen on figure 5.2. In this case, when there is not an additional entrance of oxygen, results observed are quite the same as on figure 5.1. But bigger differences appear when oxygen enters through the whole domain. As it may be seen, production of SN O increases more than the one seen on figure 5.1 while SN H is only transported through the domain. This occurs because its process rate depends on SN H , and now its presence on the domain does not decrease. Also XB,A increases more in this case, and the reason is also the no-decreasing of SN H . It can be also appreciated that oxygen presence in water is not that high as when SN H is not restricted. Here, when Oxygen Transmission Rate 1 is stabilized, SO almost stops growing because the aerobic growth of heterotrophs rate, which reduces SO presence, is bigger than before due to its dependance on SN H . Also this results have a certain biological sense. When the concentration of N H4+ 1

see section 2.2.2.

Validation

41

No oxygen’s entrance

Oxygen’s entrance 1.5

1.5

1

1

1

1

0

Ss

Ss 0.5

0.5

0

0.5 xnodo

0

1

Xbh

1.5

Xbh

1.5

0.5

0

0.5 xnodo

0

1

1.5

1

1

1

0

0.5

0

0.5 xnodo

0

1

0.5 xnodo

0

1

0

0.5 xnodo

0

1

0

0.5 xnodo

0

1

2

1

0.5

0.5 xnodo

1

0

0

0.5 xnodo

1

0

0.5 xnodo

1

Snh

1 Sno

2 Snh

1.5

Sno

3

0

1

5

1.5

0

0.5 xnodo

0.5

3

1

0

10

So

So 0.5

0

Xba

1.5

Xba

1.5

0.5

1

0

0.5 xnodo

1

0

0.5

0

0.5 xnodo

1

0

Figure 5.2: Evolution of concentration on the six equation model without and with oxygen entrance when SN H = 0

42

Mathematical modeling of heterogeneous biological wastewater processes

and N H3 nitrogen is maintained constant, then nitrification will not stop and so production of nitrate and nitrite nitrogen should be higher (equilibrium seen on expressions 2.10 and 2.12). On the oxygen case, process is similar to this one, because as nitrification takes place, waste of oxygen is given by the process.

5.2.2

Validation of the ASM1 modeling

Now the changes on the whole ASM1 will be studied. Although it does not represent all the processes given on wastewater treatments, this model is quite more complete as the other one. So from the biological point of view, results may have much more sense. Following the procedure made on last section, lets start looking at results given when oxygen’s entrance is given through the whole domain, as also when it is not. The main difference given between the six equations model (figure 5.1) and the whole ASM1 (figure 5.3), is that here, SS is produced instead of destructed. That occurs because hydrolysis of entrapped organics rate (ρ7 ) is bigger than both of heterotrophs growth (ρ1 and ρ1 ). As it may be seen, it continues growing also when oxygen is not present, but with littler velocity due to that one term of this ρ7 rate is proportional to SO while the other one is not. When analyzing both solutions, letting entrance of oxygen and not, main differences are presented on SN O , SN H , SS , XS , SN D and XN D concentrations. Differences on the first two components are the same as on the six equation model and they have been already described on last subsection. On XS and XN D , its destruction is widely increased when oxygen is presented mainly because of its hydrolysis of entrapped organics rates. On the other hand, SS and SN D production is increased. From the biological point of view, the littler destruction of XS is because under anoxic conditions (when nitrate is the only terminal electron acceptor) conversion of slowly biodegradable material into readily biodegradable one is lower. Between particulate organic nitrogen and soluble organic one, conversion follows the same patron. Under anaerobic conditions (neither oxygen nor nitrate are given), those conversion rates should stop. On figure 5.4, SN O has been forced to be zero in the whole domain, so that anaerobic conditions are partially presented. As it may be seen, when oxygen is already given, both SS and SN D increase their concentration, while when anaerobic conditions are given, SS remains constant while SN D concentra-

Validation

43

No oxygen’s entrance 1.5

2.5

1.5

1.5

1.5

1

1

1

1

0.5

Xs

1.5 Xi

Si

Ss

1

Xbh

2

0.5

0.5

0.5

0.5 0

0

0.5 xnodo

0

1

1.5

0

0.5 xnodo

0

1

1.5

0 0

0.5 xnodo

1

0

1.5

0.5 xnodo

0

1

2.5

0

0.5 xnodo

1

0

0.5 xnodo

1

0

0.5 xnodo

1

0

0.5 xnodo

1

1.5

2

0.5

0.5

0.5

1

1.5

Snh

Sno

Xp

1 So

1

Xba

1

1

0.5

0.5 0

0

0.5 xnodo

0

1

0

0.5 xnodo

0

1

2.5

0

0.5 xnodo

0

1

1.5

1.5

1

1

0

0.5 xnodo

1

0

0.5 xnodo

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

0.5 xnodo

0

1

0

0.5 xnodo

0

1

Oxygen’s entrance 1.5

2.5

1.5

1.5

1.5

1

1

1

1

0.5

Xs

1.5 Xi

Si

Ss

1

Xbh

2

0.5

0.5

0.5

0.5 0

1

1.5

0

0.5 xnodo

0

1

1.5

1 Xp

Xba

1

0.5

0 0

0.5 xnodo

1

10

2.5

8

2

6

1.5

4

0.5

0

2 0

0

0.5 xnodo

0

1

0

0.5 xnodo

0

1

2.5

0.5 xnodo

0

1

1.5

1 Snh

0.5 xnodo

Sno

0

So

0

1

0.5

0.5 0

0.5 xnodo

0

1

1.5

1.5

1

1

0

0.5 xnodo

1

0

0.5 xnodo

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

0.5 xnodo

1

0

0

0.5 xnodo

1

0

Figure 5.3: Evolution of concentration on the complete ASM1 equation model with oxygen entrance

Mathematical modeling of heterogeneous biological wastewater processes

44

1.5

2.5

1.5

1.5

1.5

1

1

1

1

0.5

Xs

1.5 Xi

Si

Ss

1

Xbh

2

0.5

0.5

0.5

0.5 0

0

0.5 xnodo

0

1

1.5

0

0.5 xnodo

0

1

1.5

0 0

0.5 xnodo

1

0

1.5

0.5 xnodo

0

1

2.5

0

0.5 xnodo

1

0

0.5 xnodo

1

1.5

2

0.5

0.5

0.5

1

1.5

Snh

Sno

Xp

1 So

1

Xba

1

1

0.5

0.5 0

0

0.5 xnodo

0

1

0

0.5 xnodo

0

1

2.5

0

0.5 xnodo

0

1

1.5

1.5

1

1

0

0.5 xnodo

1

0

0.5 xnodo

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0 0

0.5 xnodo

1

0

0.5 xnodo

1

0

Figure 5.4: Evolution of concentration on the complete ASM1 equation model with neither oxygen entrance nor nitrite presence

Validation

45

tion even decreases (due to ammonification). And XS , when anaerobic conditions are given, is increased while its conversion into readily biodegradable material is stopped but decay of autotrophs and heterotrophs is not. After observing solutions on those different situations, both on the six equations model and on the whole ASM1 formulation, it may be assumed that the numerical modeling has been properly implemented.

5.3

Verification of the two dimensional model

As the one dimensional case has already been proved, now is the turn of the two dimensional model to be verified. It’s not necessary to validate it as it has been done on the one dimensional case, but just observing if results are similar to those ones.

No oxygen’s entrance

Oxygen’s entrance 1.5

1.5

Xbh

Ss

1.5

1.5

1 0.5

1 Xbh

0 0.5

0

0

0

0.2

0.4

0.6

0.8

1

1.5

1

1

0.2

0.4

0.6

0.8

0.5

0

1

1

0

0.5

1

0

0.5

1

0

0.5

1

10

1

5

0.5

So

1.5

0

0

1.5 Xba

0.5

So

1

1 0.5

0

0.5

0

0.2

0.4

0.6

0.8

0

1

2.5

0 0

0.2

0.4

0.6

0.8

0

1

2 Sno

1 Snh

1

0.5

1.5

1.5

2 1.5

0

1

Snh

0.5

1

1 0.5

0.5

0.5 0

0

0.2

0.4

0.6

0.8

1

0

0 0

0.2

0.4

0.6

0.8

1

0

0.5

1

0

Figure 5.5: Evolution of concentration on the six equation model without and with oxygen entrance on two dimensions On figures 5.5 and 5.6 solutions have been presented showing only x-axe and concentration, so that it is easier to compare results from both models. Results have been calculated using as initial condition c0 = 1. It may be seen that results on both cases are quite similar, so they can be accepted and the two dimensional model can be assumed as correct also.

Mathematical modeling of heterogeneous biological wastewater processes

46 No oxygen’s entrance 1.5

1.5

1.5

1

1

1

1

1

0

0.5

0

0.5

0

1

1.5

Xs

Xi

Si 0.5

0.5

0.5

0

0.5

0

1

1.5

Xbh

1.5

Ss

1.5

0.5

0 0

0.5

1

0

1.5

0.5

0

1

2.5

0

0.5

1

0

0.5

1

0

0.5

1

0

0.5

1

1.5

2

0.5

0.5

0.5

1

1.5

Snh

Sno

Xp

1 So

1

Xba

1

1

0.5

0.5 0

0

0.5

0

1

0

0.5

0

1

2.5

0

0.5

0

1

1.5

1.5

1

1

0

0.5

1

0

0.5

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

0.5

0

1

0

0.5

0

1

Oxygen’s entrance 1.5

1.5

1.5

1

1

1

1

1

0.5

0

1

1.5

0

0.5

1.5

1 Xp

Xba

1

0.5

Xs 0

1

0 0

0.5

1

0

10

2.5

8

2

6

1.5

4

0.5

2 0

0

0.5

0

1

0

0.5

0

1

2.5

0.5

0.5

0

1

1.5

1 Snh

0

0.5

0.5

Sno

0.5

So

0

Xi

Si 0.5

Xbh

1.5

Ss

1.5

1

0.5

0.5 0

0.5

0

1

1.5

1.5

1

1

0

0.5

1

0

0.5

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

0.5

1

0

0

0.5

1

0

Figure 5.6: Evolution of concentration on the complete ASM1 equation model with oxygen entrance

Validation

5.4

47

Defining τ

It is not that easy to see instabilities on results produced with a initial condition equal one, but when it is defined as zero, those are quite important. In figure 5.7, where initial condition c0 = 0 has been used, they can be seen. Those instabilities can be better appreciated on the time-concentration graphic, so on this figure concentration along the time on the exit boundary is presented.

1.5

1

1

Ss

Xbh

1.5

0.5 0

0.5

0

1 t

0

2

1.5

0

1 t

2

0

1 t

2

0

1 t

2

10

So

Xba

1 5

0.5 0

0

1 t

0

2

2.5

1.5 1

1.5

Snh

Sno

2

1

0.5

0.5 0

0

1 t

2

0

Figure 5.7: Results of the six equation model using initial condition c0 = 0 Then, the need of a stabilization technique is evident. And so a stabilization parameter τ shall be defined. There is plenty of different formulations for the wide range of equations and its different conditions. Just for unsteady convection-diffusionreaction problems several of them can be found. Some examples can be found in [4]. When steady state is the one which concerns in the solution, next expressions can be used. ¶ µ h 1 (5.1) τ= coth(Pe ) − 2a Pe

Mathematical modeling of heterogeneous biological wastewater processes

48

µ ¶−1 h 1 h τ= 1+ + σ 2a Pe 2a à µ ¶2 !−1/2 h 9 h τ= 1+ 2 + σ 2a Pe 2a

(5.2)

(5.3)

Due to nature of the reaction term, which is not constant, it may not be considered on the stabilization parameter, so expressions can be rewritten as µ ¶ 1 h coth(Pe ) − (5.4) τ= 2a Pe µ ¶−1 h 1 τ= 1+ 2a Pe µ ¶−1/2 h 9 τ= 1+ 2 2a Pe

(5.5) (5.6)

On figure 5.8, three different graphics of SN O concentration can be seen, which are the result of the six equations model. Each one of them compares the resolution of the Galerkin solution without stabilization versus the use of a stabilization technique where each τ is used. There, τa corresponds to the use of expression 5.4, τb is expressed on formulae 5.5 while τc is the one on expression 5.6. τa

τb

τc 2.5

2

2

2

1.5

1.5

1.5 Sno

2.5

Sno

2.5

1

1

1

0.5

0.5

0.5

0

0

0.5

1 t

1.5

2

0

0

0.5

1 t

1.5

2

0

0

0.5

1 t

1.5

2

Figure 5.8: Results of SN O on six equations model using initial condition c0 = 0 calculated with stabilization technique using different expressions from τ As it may be seen, stabilization is achieved, and on the three, results are almost the same. But it adds a lot of diffusion, which could be reduced. That has a lot of sine, because those expressions are not made for this particular problem. So a

Validation

49

readjustment should be made, which may be seen on figure 5.9. To obtain these results τ has been divided by different numbers (5, 10 and 25) in order to obtain a better value, that means trying to find a stabilization parameter which does not give instabilities but also does not add too much diffusion. Results have been also presented on SN O solution obtained by the six equations model.

2.5

2

Sno

1.5

1

0.5

0

0

0.5

1 t

1.5

2

Figure 5.9: Comparison on SN O results of solutions reached without stabilization and using τ /5, τ /10 and τ /25 It may be seen that solution presented on red does not have any perturbation, and its added diffusion is littler than the one in green. This one corresponds to the τ /10 solution. Then the stabilization parameter which shall be used in this particular problem under the conditions before described is ¶ µ 1 h 1 (5.7) τ= coth(Pe ) − 10 2a Pe But from a more generic point of view, τ should be written as ¶ µ h 1 τ =α coth(Pe ) − 2a Pe where α could be another value when some aspects of the problem change.

(5.8)

50

Mathematical modeling of heterogeneous biological wastewater processes

Chapter 6 Results Once the model has been implemented and the goodness of it has been proved, it is time to present results. On the first section of this chapter, results presented are those of the flow of water on an homogeneous domain, either where oxygen’s entrance is allowed or it is not. Later, results presented are quite more realistic because oxygen’s entrance is allowed only on the upper part of the domain. In this case, new oscillations have been found, and ways to reduce its presence on the solution tried out.

6.1

Results on a homogeneous domain

Two different models have been implemented, and different domains have been used (one dimension and two dimensions). But all this steps have been made only to reach a model from the whole Activated Sludge Model No. 1 in a two dimensions domain. So here results presented will refer just to this case. Results on modeling the ASM1 on Ω = [0, 1]x[0, 1] as it was defined on chapter 3 are shown next. Both cases, letting oxygen enter through the whole domain and avoiding its entrance, are considered. Firstly results when oxygen entrance is allowed on the whole domain are presented. Next, results shown are the ones when the only oxygen which gives on the wetland is the one that enters solved on the wastewater. It can be seen that results do not have oscillations thanks to the stabilization parameter chosen before.

Mathematical modeling of heterogeneous biological wastewater processes

52

1.5

1.5

1.5

1

1

1.5

0 0

1 0.5

1

0.5

0 0

0.5 1 0

1.5

1.5

1

1

Xs

Xi

Ss

Si 0.5

1 0.5

0.5

0 0

0.5

Xbh

2 1

1 0

1 0.5

0 0

0.5 1 0

1 0.5

1 0.5 0 0

0.5 1 0

10

1 0.5

0.5

1 0

0.5

1 0

1.5

5

0.5

0 0

1 0.5

0 0

0.5 1 0

1 0.5

0 0

0.5 1 0

1 0.5

Snh

Sno

0.5

So

Xp

Xba

2 1 0 0

0.5 1 0

1.5

1

0 0

0.5

0.5

1 0

0.5

1 0

1 0.5 1 0.5

1.5

1 0 0

1 0.5

1

Salk

Xnd

Snd

2 0.5 0 0

0.5 1 0

1 0.5

1 0.5 0 0

0.5 1 0

1 0.5

Figure 6.1: Evolution of concentration on the ASM1 model with oxygen entrance on the whole section of the wetland

1.5

2.5

1.5

1.5

1.5

1

1

1

1

0.5

Xs

1.5 Xi

Si

Ss

1

Xbh

2

0.5

0.5

0.5

0.5 0.5

0

1

1

1

0.5

0

0.5

0

1

10

Xp 0.5

0

1

So

1.5

Xba

1.5

0

8

2

6

1.5

2 0

0

0.5

0

1

0

0.5

0

1

2.5

0.5

0

1

2.5

4

0.5

0

0

0.5

1

0

0.5

1

1.5

1 Snh

0

Sno

0

1

0.5

0.5 0

0.5

0

1

1.5

1.5

1

1

0

0.5

1

0

0.5

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

0.5

1

0

0

0.5

1

0

Figure 6.2: Evolution of concentration on the ASM1 model with oxygen entrance along the wetland

Results

53

1.5

2.5

1.5

1.5

1.5

1

1

1

Xs

1.5 Xi

Si

Ss

1

1

0.5

Xbh

2

0.5

0.5

0.5

0.5 1 t

0

2

1.5

0

1 t

0

2

1.5

Xp 0.5

So

1

Xba

1

0

1 t

10

2.5

8

2

6

1.5

4

0.5

0

2

2 0

0

1 t

0

2

0

1 t

0

2

2.5

0

1 t

0

2

0

1 t

2

0

1 t

2

1.5

1 Snh

0

Sno

0

1

0.5

0.5 0

1 t

0

2

1.5

1.5

1

1

0

1 t

2

0

1 t

2

0

Salk

Xnd

Snd

2 1.5 1

0.5

0.5

0.5 0

0

1 t

0

2

0

1 t

0

2

Figure 6.3: Evolution of concentration of each component on the exit boundary when oxygen entrance is given on the whole domain along the time

1.5

1.5

1.5

1

1

1.5

0 0

1 0.5

1 0

Xs

1 0 0

0.5

Xi

Ss

Si 0.5

0.5 1 0.5

0 0

0.5

Xbh

2 1

1 0

1.5

1.5

1.5

1

1

1

0.5 1 0.5

0 0

0.5 1 0

1 0.5

1 0.5 0 0

0.5 1 0

1 0.5

0.5

1 0

0.5

1 0

1.5

0 0

1

0 0

0.5 1 0

0.5 1 0.5

0 0

0.5 1 0

1 0.5

1 0 0

0.5 1 0

1.5

Snh

0.5

0.5

Sno

0.5

So

Xp

Xba

2

1 0.5

0.5

1 0

0.5

1 0

1 0.5 0 0

1 0.5

1.5

1 0 0

1 0.5

0.5 1 0

1

Salk

Xnd

Snd

2 0.5 0 0

1 0.5

0.5 1 0

1 0.5 0 0

1 0.5

Figure 6.4: Evolution of concentration on the ASM1 model with oxygen entrance on the whole section of the wetland

Mathematical modeling of heterogeneous biological wastewater processes

54

1.5

2.5

1.5

1.5

1.5

1

1

1

1

0.5

Xs

1.5 Xi

Si

Ss

1

Xbh

2

0.5

0.5

0.5

0.5 0

0

0.5

0

1

0

0.5

0

1

1.5

1.5

1.5

1

1

1

0

0.5

0

1

0

0.5

0

1

2.5

0

0.5

1

0

0.5

1

1.5

0.5

0.5

0.5

1 Snh

Sno

So

Xp

Xba

2 1.5 1

0.5

0.5 0

0

0.5

0

1

0

0.5

0

1

2.5

0

0.5

0

1

1.5

1.5

1

1

0

0.5

1

0

0.5

1

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

0.5

0

1

0

0.5

0

1

Figure 6.5: Evolution of concentration on the ASM1 model with oxygen entrance along the wetland

1.5

2.5

1.5

1.5

1.5

1

1

1

1

0.5

Xs

1.5 Xi

Si

Ss

1

Xbh

2

0.5

0.5

0.5

0.5 0

0

1 t

0

2

0

1 t

0

2

1.5

1.5

1.5

1

1

1

0

1 t

0

2

0

1 t

0

2

2.5

0

1 t

2

0

1 t

2

1.5

0.5

0.5

0.5

1 Snh

Sno

So

Xp

Xba

2 1.5 1

0.5

0.5 0

0

1 t

0

2

0

1 t

0

2

2.5

0

1 t

0

2

1.5

1.5

1

1

0

1 t

2

0

1 t

2

0

Salk

1.5

Xnd

Snd

2

1

0.5

0.5

0.5 0

0

1 t

2

0

0

1 t

2

0

Figure 6.6: Evolution of concentration of each component on the exit boundary when oxygen entrance is given on the whole domain along the time

Results

6.2

55

Results on a heterogeneous domain (or epilogue)

Models presented until now which separate both situations, when oxygen enters through the whole domain and when it does not, may be interesting to observe how do the implementations work. But this models have not much to do with wetlands reality. There, as briefly explained before, oxygen enters through the superior boundary, as also from roots of the existing plants on the upper part of the domain. So a quite more realistic implementation has been made. It consists on having different reaction terms depending on the position of the point P of Gauss within the domain. While ones use the term presented on the ASM1 as j νij ρj , others add to this term the Oxygen Transmission Rate (equation 2.17).

6.2.1

Discontinuous OTR reaction term

On the first case, the border has been defined on the upper 40% of the domain, where the reaction term presents the OTR, while the other 60% does not. On this 60% of the domain, this discontinuous reaction term is imposed with values equal zero, while on the others, its values are equal one. As it may be seen on next figures, this option presents some computational problems, and its computational cost is so high that just till final time 0.2 has been computed. On figure 6.7 it can be seen that the discontinuity on oxygen entrance by the reaction term produces a big oscillation on several unknowns of the problem, where several variables take negative values. Those negative values could be filtered using a modification of the time integration scheme imposing at each time step that all negative values equal zero. However this would affect the convergence of the time integrator, and, moreover, the tendency to present negative values will remain and the overall reaction term will be affected.

Mathematical modeling of heterogeneous biological wastewater processes

56

0 0

5

0 1 0 0.5 0 0.5 0 0.5 0.5 1 0 mlss autotrophic mlss heterotrophic 1

0.5

2

1 0

1.5

1.5

1

1

1

0.5

0 0

1 0 0.5 0 1 0 Sno

0.5

1 0 0.5 0 1 0 bod

0.5 1 0 Xs

0 1 0

0.5

0.5 1 0 Xp

1 0.5

1 0 0

0.5

1 0.5 1 0 Snh

0.5

1 0 0 0.5 1 0 Snd

0.5

1 0.5 1 0 Xnd

0.5

0.5 1 0 Xi

1.5

1.5 1

1

0.5

0.5

0.5

1 0 0.5 0 1 0 Salk

1 0 0

1

1

0.5

0.5 1 0 do

1.5

0.5

0.5

0.5

2

1.5

0 0

0.5

0.5

0.5

0.5

1

1

1

1

1.5

10

1.5

1.5

1.5

0.5

0.5 1 0 Si

1

0 0

1

Figure 6.7: Evolution of concentration of each component when oxygen entrance is given on the upper 40% of the domain (final dimensionless time = 0.2)

Results

6.2.2

57

Continuous OTR reaction term

Both, because of the problems of the discontinuous reaction term and because actually a transition between both zones, the one with oxygen and the one without, may be expected, a second modeling option has been tested. A linear transition from zero to one in reaction term is imposed on y direction from 0.6 to 0.8. Provided that a vertical discretization on 20 elements is used, the transition includes four elements. A reaction term with value one is imposed on the 20% upper domain. Results may be seen on figure 6.8, where computation was executed until final dimensionless time 0.5, and on figure 6.9, where final dimensionless time was 2.

1.5

1.5

1

1

0.5

0.5

0 0

0.5

1 0 0

0 1 0 0.5 0.5 1 0 mlss autotrophic

1.5

1

1

0.5

1 0.5 1 0 bod

0 0

0.5

1 0.5 1 0 do

0.5

0 0

0.5

1 0.5 1 0 Xp

0 0

1.5 1

1

0.5

0.5

0.5

1 0.5 1 0 Salk

0 0

1 0.5 1 0 Snh

0.5

1 0.5 1 0 Xnd

1

1

0.5

0.5

1.5

0.5

0.5 1 0.5 1 0 Xs

0 0

2

1.5

0 0

0.5

1

0.5 0 0

5

0.5

1.5

1 0.5 1 0 Sno

1

1

0 1 0 0.5 0.5 1 0 mlss heterotrophic

2

1.5

10

1.5

0.5

1 0.5 1 0 Snd

0.5

1 0.5 1 0 Xi

0 0

1.5

0.5

1 0.5 1 0 Si

0 0

Figure 6.8: Evolution when oxygen entrance is given on the upper 20% of the domain, while on the lowest 60% it is not, and on the 20% in between there is a linear transition (final dimensionless time = 0.5). As it should be expected, results on the upper part of the domain are as those obtained with O2 presence, and in the lower part as those without it.

Mathematical modeling of heterogeneous biological wastewater processes

58

1.5

1.5

1.5

1

1

1

0 0

0.5

2

0 0 1 0.51 0 0.5 0.5 0 0.5 1 0 mlss autotrophic 1 0 mlss heterotrophic

1.5

1.5

1

1

1

0.5

0 0

0 1 0.5 0 1 0 Sno

0.5

1 5 1 0.5 1 0 bod

0.5

0 0

1.5

1

1

0.5

0.5

1 0.5 1 0 Xp

0.5

0 0

0.5

0.5 1 0 do

0 0

1

0.5

1 0.5 1 0 Snh

2

1

1

0.5

0 0

0.5

1 0.5 1 0 Snd

0.5

0.5 1 0 Xi

0 0

0.5

1 0.5 1 0 Xnd

1.5 1 0.5

0.5 1 0.5 1 0 Salk

0 0

1.5

0.5 1 0.5 1 0 Xs

1.5

0 0

0.5

0.5

0.5

0.5

1.5

10

0.5

0.5 1 0 Si

1

0 0

1

Figure 6.9: Evolution when oxygen entrance is given on the upper 20% of the domain, while on the lowest 60% it is not, while on the 20% in between there is a linear transition (final dimensionless time = 2.0).

Results

59

It is remarkable that the small oscillations presented during the evolution of the process (figure 6.8) are also presented on the stationary state (figure 6.9), although there they are concentrated around sharp profiles. On figure 6.10 a finer discretization along the y direction can be seen, where 40 elements have been used. So now transition is made through eight different elements.

1.5

1.5

1.5

1

1

1

0.5

0.5

0.5

0 0

0.5

1

1.5

1.5

1

1

1 0.5 1 0 Sno

0 0

1 0 0.5 0 1 0 bod

0.5

0.5 1 0 do

1 0 0

0.5

1 0.5 1 0 Snh

0.5

1 0.5 1 0 Xnd

1.5 1 1

0.5 0.5

1 0.5 1 0 Xs

0 0

0.5

0.5 1 0 Xp

1 0 0

1.5

1

1

1

0.5

0.5

0.5

0.5

1 0 0 0.5 1 0 Salk

0.5 0.5

1 0 0.5 1 0 Snd 0

0.5

1 0.5 1 0 Xi

1.5

1.5

0 0

1 0.5

2

0.5 0.5

1.5

5

0 1 0 1 0.5 0 0.5 0 0.5 0.5 1 0 mlss autotrophic 1 0 mlss heterotrophic

2

0 0

10

0.5

0.5 1 0 Si

1 0 0

Figure 6.10: Evolution when oxygen entrance is given on the upper 20% of the domain, while on the lowest 60% it is not, and on the 20% in between there is a transition, while velocity is increased by a factor 5. Here oscillation on most of the variables on the transition between both zones,the one with O2 presence and the one without, does not decrease but increase. This results confirm that the oscillation is produced due to the sharp profile of oxygen on the y direction. As the element size has been decreased (halved on on y direction), the oxygen profile is better defined. However, the profile in this direction increases suddenly with the presence of O2 , and the first element affected by it concentrates the main

Mathematical modeling of heterogeneous biological wastewater processes

60

part of this increase. The smaller the element is, the sharper this variation will be, and this produces the oscillations on all the other variables. Oxygen discontinuity can be directly related to different relative velocities between oxygen source/reaction term and convective transport. This can be confirmed on figure 6.11, where horizontal velocity has been increased by a factor 5.

1.5 1.5

1.5

1

1

0.5

0.5

10

1.5

1

0 0

0.5

1 0

0.5

1 0 0 0.5 0.5 mlss heterotrophic

2.5

1 5

1 0

0.5

1 0 0 0.5 0.5 mlss autotrophic

1.5

1 0.5 1 0 bod

1.5

0.5

1 0.5 1 0 do

2.5

2

0 0

0.5

1 0.5 1 0 Snh

0.5

1 0.5 1 0 Xnd

1.5

2 1

1.5 1

1

0.5

1

1.5 1

0.5

0.5 0 0

0 0

0.5

0.5 0.5

1 0.5 1 0 Sno

0 0

0.5

1 0.5 1 0 Xs

0 0

0.5

1 0.5 1 0 Xp

0 0

1.5

1.5

1.5

1

1

1

0.5 0 0

0.5

0.5

1 0.5 1 0 Salk

0 0

0.5

1 0.5 1 0 Snd

0.5

1 0.5 1 0 Xi

0 0

0.5

0.5

1 0.5 1 0 Si

0 0

Figure 6.11: Evolution of concentration of each component when oxygen entrance is given on the upper 20% of the domain, on the lowest 60% no oxygen entrance is allowed, while on the 20% in between there is a transition (final dimensionless time = 2.0). Now oscillations have been widely decreased.

Chapter 7 Conclusions and future work 7.1

Conclusions

In the present master thesis, a new way on the mathematical modeling of wastewater treatment processes which occur in subsurface flow constructed wetlands has started to be developed. Mathematical models which already had been used in other contexts of wastewater treatment, are not much useful for those treatment plants. From the biological part of the modeling, the Activated Sludge Models from IWA have been seen as the more adequate ones to express rates. Not only due to its international acceptance and reputation, but also to its presentation, on a clear and workable matrix. This allows an easily formulation of the components rates. Also the use of continuous functions on its process rates helps on the resolution of the equations system. Although it is a really good model, not all the aspects may help. The worst of them is that its process rates are coupled. This leads to more slow resolutions of the system. Obviously, as nature processes are coupled, this could not be simplified. When refereing strictly to the mathematical part of the model, here the finite element method has been applied, in contrast from the cases found on bibliography, where finite differences where normally used [19] [8]. There, normally domain’s discretization was made doing partitions in only one axe. In this work, a two dimensional section of the whole domain has been used, and its division has been done on both directions. This may allow a better approximation on the differentiation of aerobic and anaerobic zones on the wetlands, which is really useful to understand process there given.

Mathematical modeling of heterogeneous biological wastewater processes

62

For the resolution of the spatial discretization, the Galerkin formulation has appeared to be a reasonable method to solve it, although a stabilization technique is necessary. The three presented on this work (Stream Upwind Petrov-Galerkin, Galerkin-Least Squares and Subgrid Scale) apply same results as they end up with the same expression due to the election of linear elements for the discretization. So the three of them are adequate for the resolution. Otherwise, on the spatial discretization the Runge-Kutta-Fehlberg method appears to be a good method due to its variable stepsize. The implementation made here using fourth and fifth order methods should probably be reduced to lower order ones. This would lead to a resolution with minor computational cost, although damages accuracy should be studied. Observing first results acquired, it seems that the model solves the problem in a proper manner. So it can be concluded that the mathematical modeling of the ASM1 on a homogeneous domain works properly. Posterior adaptations to get a more realistic situation, where oxygen’s entrance is only allowed in a part of the domain, have shown spurious oscillations on the discontinuity produced due to this difference. These can be reduced implementing a linear transition between both zones, although they don’t disappear. Also increasing the horizontal velocity of the flow reduces their presence.

7.2

Future work

Through all this work, one and two dimension models have been used. Extending this to a three dimensional one, although it may bring a bit more of information, does not seem much important. Firstly because its computational cost would increase widely, but mainly because much more adaptations should be done before in order to obtain more realistic values. Lets see some of them that have been outlined along the present master thesis.

• Posterior IWA models could be used, either the ASM2, the ASM2d or the ASM3. • Other geometries of the domain, as also flow models, should be considered, attending to more realistic situations.

• Adsorption could be considered, not only the linear expression but also a more complex one which could represent much better this process. • Other spatial discretization techniques (instead of the GLS, SUPG and SGS) should be proved, in order to avoid instabilities. The Least-squares formulation, which has been briefly described on chapter 4 may be a good one to start with. • Order of the time discretization schema could be reduced, starting probably with a Runge-Kutta-Fehlberg 2-3. • Stabilization of discontinuities produced due to the entrance of oxygen on the heterogeneous domain shall be studied. • The upper boundary condition for some of the components should be modified to represent the fact that some of them escape on a gaseous manner. • Finally, some more aspects relating to these wetland operation would be necessary. Between others hydraulic parameters and deposition of the particulate components. Although the list seems quite short, there is a lot of work to do, but it seems not only possible but also useful the use of finite elements method to study this black box that subsurface flow constructed wetlands already are.

Acknowledgements I beg your pardon, but on this part I will use catalan, my language and the one I use with the people whom I would like to thank. It is not that I don’t want others to understand this, but just that due to my (poor) English level, these acknowledgements would be colder than they should. Primerament, agrair a l’Agust´ı P´erez Foguet, el tutor, l’esfor¸c i les hores de dedicaci´o que ha posat en aquesta tesina, aix´ı com la seva paci`encia en front dels meus dubtes i desconcerts. Ha estat un llarg per´ıode de treball en que sempre que ho he necessitat m’ha ajudat. Moltes gr`acies per la teva gran disponibilitat, pel teu treball i per la teva confian¸ca. Tamb´e agrair a l’Esther Sala Lardies la seva paci`encia en atendre’m cada cop que l’Agust´ı no hi era i algun dubte em sorgia. Moltes gr`acies per rebre’m sempre amb un somriure i amb molta paci`encia. ´ de just´ıcia agrair tamb´e en Roger Vilaseca i Cabo, per la seva ajuda inicial Es en els meus dubtes sobre els EF, i el cop de m`a que m’ha donat amb el LATEX 2ε . Gr`acies, tot i que el meu agra¨ıment m´es sincer per a tu b´e tamb´e despr´es. Cal fer un especial agra¨ıment a tots els grans amics que m’han acompanyat en el (llarg) cam´ı que han estat aquests estudis, no nom´es els darrers dos anys amb la tesina, sin´o els gaireb´e 8 anys de carrera. Primerament als amics de ”la casa”, per les moltes hores compartides a l’Escola, pels moments de nervis (sin´o hist`eria) en `epoques d’ex`amens, per les hores al bar i les hores de festa, per les esquiades, les hores de biblioteca, les converses intimistes... Gr`acies per ser com sou i aguantar-me tal com s´oc, amb les meves desaparicions i el meu magn´ıf ic humor. Per`o tamb´e als de PACCS que, tot i que de forma m´es distant, tamb´e han patit la meva carrera. I en especial a en Xavi, la Mireia i la Guime per compartir l’estudi

(i uns quants caf`es) al llarg d’aquests vuit anys. Gr`acies a tots per haver-me acompanyat en el cam´ı i haver-me deixat ser el vostre intendent. I com no, a aquells qui han patit m´es tant aquests dos anys de tesina com els vuit d’estudis, la meva fam´ılia. A la tia, pels seus `anims. A les meves germanes, l’Ester i la S´ılvia, i als meus cunyats, en Marco i en David, pel suport que m’han donat sempre, per la seva voluntarietat i per intentar treure ferro als patiments dels meus pares. I sobretot als meus pares, que m’han ajudat en tot el que han pogut i m´es, i han patit cada examen i cada nota tant o m´es que jo. Moltes gr`acies per ser-hi sempre i per aguantar totes les meves idees tot i que les trob´essiu f´ora de lloc (per dir-ho finament). A tots vosaltres, moltes gr`acies.

Bibliography [1] J. Bolmstedt and G. Olsson. Exercise 1. control of biological wastewater treatment, 2002. [2] F.L. Burton and G. Tchobanoglous. Wastewater Engineering. Treatment, disposal, and reuse. McGraw-Hill, 1972. [3] I.G. Currie. Fundamental Mechanics of Fluids. Marcel Dekker Inc., New York, third edn edition, 2003. [4] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. John Wiley & sons., West Sussex, England, 2003. [5] R.L. Droste. Theory and Practice of Water and Wastewater Treatment. John Wiley & sons., 1997. [6] H. Vanhooren et al. West: modelling biological wastewater treatment. Journal of Hydroinformayics, pages 27–50, January 2003. [7] M. Henze et al. Activated Sludge Models ASM1, ASM2, ASM2d and ASM3. IWA Publishing, 2000. [8] Y. Zeng et al. Hydraulic modelling and axial dispersion analysis of uasb reactor. Biochemical Engineering Journal, Vol. 25:113–123, 2005. [9] Water Environment Federation. Biological Nutrient Removal (BNR) Operation in Wastewater Treatment Plants. McGraw-Hill, 2005. [10] E. Fehlberg. Low-order classical runge-kutta formulas with stepsize control and their application on some heat transfer problems. 1969. [11] Y. Shen G.F. Carey and R.T. McLay. Parallel conjugate gradient performance for least-squares finite elements and transport problems. International Journal for Numerical Methods in Fluids, pages 1421–1440, 1998.

[12] B. Roig J. Donea and A. Huerta. High-order accurate time-stepping schemes for convection-diffusion problems. Comput. Methods Appl. Mech. Engrg., pages 249–275, 2000. [13] J. Morat´o J. Garc´ıa and J.M. Bayona. Nuevos criterios para el dise˜ no y operaci´ on de humedales construidos. Ediciones CPET, 2004. [14] A. Hern´andez Muoz. Depuraci´ on de aguas residuales. Paraninfo, S.A., 1990. anica de medios continuos para [15] X. Oliver and C. Agelet de Sarac´ıbar. Mec´ ingenieros. Edicions UPC, Barcelona, second edn edition, 2000. [16] G. Olsson and B. Newell. Wastewater Treatment Systems. Modelling, diagnosis and control. IWA Publishing, London, 1999. [17] F.Y. Wang T.T. Lee and R.B. Newell. On the evaluation of the exit boundary condition for the axial dispersion bioreactor system. Chemical Engineering and Technology, Vol. 21. [18] F.Y. Wang T.T. Lee and R.B. Newell. Dynamic modelling and simulation of activated sludge process using orthogonal collocation approach. Water Research, Vol. 33, No. 1:73–86, 1999. [19] F.Y. Wang T.T. Lee and R.B. Newell. Dynamic simulation of bioreactor systems using orthogonal collocation on finite elements. Computers and Chemical Engineering, Vol. 23:1247–1262, 1999.