Next Article in Journal
Assessing Urban Flooding Extent of the Baunia Khal Watershed in Dhaka, Bangladesh
Previous Article in Journal
Comparative Study for Daily Streamflow Simulation with Different Machine Learning Methods
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

E-Learning Proposal for 3D Modeling and Numerical Simulation with FreeFem++ for the Study of the Discontinuous Dynamics of Biological and Anaerobic Digesters

by
Saulo Brito-Espino
1,*,†,
Tania García-Ramírez
1,
Federico Leon-Zerpa
1,
Carlos Mendieta-Pino
2,
Juan J. Santana
2 and
Alejandro Ramos-Martín
2
1
Institute for Environmental Studies and Natural Resources (i-UNAT), University of Las Palmas de Gran Canaria (ULPGC), 35017 Las Palmas de Gran Canaria, Spain
2
Department of Process Engineering, Universidad de Las Palmas de Gran Canaria, 35017 Las Palmas de Gran Canaria, Spain
*
Author to whom correspondence should be addressed.
Current address: Sagasta, 86, 5°, 35008 Las Palms, Spain.
Water 2023, 15(6), 1181; https://doi.org/10.3390/w15061181
Submission received: 1 March 2023 / Revised: 14 March 2023 / Accepted: 14 March 2023 / Published: 18 March 2023
(This article belongs to the Section Wastewater Treatment and Reuse)

Abstract

:
This work presents an original 3D code in FreeFem++ to recreate the behavior of anaerobic microorganisms in non-stirred anaerobic reactors with an intermittent feed. The physical and biochemical phenomena have been considered using a mathematical model based on a set of partial differential equations: Stokes, advection–diffusion, and diffusion–reaction. The description of the anaerobic metabolism was carried out by implementing the structured AMD1 model. The Galerkin finite element method has been used to solve the partial differential equations defined in the model. Finally, the methodology and procedures are presented by means of a concrete example. Thanks to the inclusion of this e-learning tool for use in virtual laboratories, it is possible to improve the understanding of engineering students on the functioning of the metabolism that takes place inside non-stirred anaerobic reactors that are fed discontinuously. This proposal reinforces to students, in a transversal way, both environmental sensitivity and awareness of the circular economy focused on the implementation of natural wastewater treatment systems in rural areas.

1. Introduction

Anaerobic digestion is arousing great interest in the agricultural and livestock sectors. The use of eco-technologies based on anaerobic processes for wastewater treatment is presented as a sustainable alternative to conventional systems due to their low energy demand and their capacity to produce biogas. Non-Stirred Anaerobic Reactors (NSARs) appear to be ideal for isolated farms where the accesses to sewage network and electricity supply are extremely difficult or nearly impossible. In this context and in order to carry out the commitment of engineering education to environmental protection, it is crucial to equip engineers with the necessary knowledge for the design of resilient environmental technologies with low greenhouse gas emissions [1].
In educational settings, the use of information and communication technologies as a pedagogical tool allows students, on the one hand, to access information and, on the other, to communicate and collaborate both with other students and with the teaching professional [2]. An application of all this can be seen in the field of engineering, where a high degree of abstraction is required to learn important concepts, properties and parameters, all with a certain degree of complexity. For these cases, Virtual Laboratories (VLs) are introduced as an effective tool to be used as an alternative and/or a complement to real laboratories within the educational system, since real laboratories have a few restrictions in terms of equipment, cost and number of experiments [3]. Through VLs, inquiry-based learning is encouraged where students must find solutions to given problems through a process of investigation [4].
Learning about anaerobic digestion processes in engineering education is justified since it is a renewable energy resource with an important contribution to make in the effort to mitigate climate change [5]. However, it also allows us to find tailor-made solutions for wastewater treatment in agricultural and livestock farms, thus contributing to sustainable development and the circular economy. The development of VLs, based on mathematical models, capable of simulating anaerobic metabolisms inside bioreactors, allows students to acquire rich experience in terms of design, control and optimization, enabling significant learning about anaerobic digestion processes and the operation of biological reactors [6].
The application of mathematical models for the description of anaerobic reactors was settled as a trend since the end of the last century. Initially, studies were limited by the lack of knowledge of the complex microbial ecology and the limited computer tools available at the time [7,8,9]. However, over the years, new models, applied to several research works, have made it possible to advance in the description of the problem [8,10]. It is important to highlight the publication, in 2002, of the “Anaerobic Digestion Model No. 1 (ADM1)” [11], which has since become a universal platform for modeling anaerobic digestion processes [12]. In recent years, many of the initial barriers have been overcome thanks to new technologies, instrumentation, and software applied to Computational Fluid Dynamics (CFD) modeling [8,13]. In educational settings [14,15], many of these models have been used to promote inquiry-based learning. The use of simulation software, such as Comsol Multiphysic, OpenFOAM and C++ libraries, has been fundamental for the visualization and experimentation of the model results [16]. The vast majority have been developed for stirred anaerobic reactors in which substrate homogeneity and the uniformity of the microbial biomass are obtained by mixing the entire fluid within the vessel. However, the study of NSARs requires a much more complex analysis, with non-linear parameters and coefficients that vary both in space and time.
The objective of this work is to provide a new tool based on CFD for the design and the simulation of NSARs in order to be used for studying the processes of the breakdown of organic matter that take place inside it. The proposed reactor type in this model is one of those reactors with a low level of mechanization, which is suitable for environmentally friendly facilities and whose activity does not involve an external energy supply. The objective is to develop a tool for educational purposes to provide a critical overview of metabolic processes and microbial interactions along the flow. This study is accompanied by the algorithm that has been designed in the FreeFem++ software. This, like the other software used in this study, is an open-source platform, which offers possibilities for both students and teachers to access the source code, distribute it, and to use it for different purposes. The computational algorithm, included in this manuscript, contains the governing equations of the problem represented in the weak forms. The Galerkin finite element method has been considered to approximate the governing equations, second-order partial differential equations, defined in the model.
Physical phenomena, advection–diffusion, and biochemical processes involved in each phase of anaerobic digestion are taken into account within the model. The proposed methodology allows the exploration of the evolution, in time and space, of the substrate and microorganisms inside the NSARs, taking into account that the results will be strongly influenced by both the reactor geometry and the initial and boundary conditions.
This methodology is presented as a suitable digital tool to promote meaningful learning for students in the continuous and discontinuous dynamics of anaerobic biodigesters. This proposal aims to improve the competences and skills of the teaching team to design active learning strategies in teaching laboratories within the field of engineering, allowing them to learn new methods and theories, and providing them with the versatility to adapt to different scenarios within the framework of respect for the environment and the circular economy.

2. Materials and Methods

The model is designed to simulate NSARs. The complexity of the model lies in the fact that it is a distributed parameter model in which the state variables vary both in time and space. The model is characterized by the use of partial differential equations.

2.1. Governing Equations

2.1.1. Stokes Equations

Stokes equations define the three-dimensional motion of the viscous fluid inside the SNARs. These equations are all derived from the linear approximations of the steady-state Navier–Stokes equations [17,18]. The Stokes equations are used for scenarios where the inertial advective forces are small compared to the viscous forces. The strong form read: for a domain Ω , find the velocity u and the pressure p satisfying;
ν Δ u + p = F f o r x , y , z ϵ Ω u = 0 f o r x , y , z ϵ Ω

2.1.2. Advection–Diffusion Equation (ADE)

This equation has been used to describe the mixing behavior in the reactor feed process. The strong ADR formulation is depicted below.
ϕ t D Δ ϕ + u ϕ = F f o r x , y , z ϵ Ω

2.1.3. Diffusion Reaction Equation (DRE)

DRE (Equation (3)) is designed to describe the behavior of the mixture in the period between reactor loading/unloading. This is a much longer time interval than the ADE [19].
ϕ t D Δ ϕ + f ( ϕ ) = F f o r x , y , z ϵ Ω

2.2. Boundary and Initial Conditions

Dirichlet and Neumann boundary conditions (Figure 1) were imposed for the three governing Equations (1)–(3).
The initial concentrations of cells and substrates were set to zero (kg·m−3) in the entire domain.

2.3. Implementation of ADM1

The integration of the biochemical reactions into the CFD mode has been completed through the term reactant ( f ( ϕ ) ) of DRE (Equation (3)). Following the ADM1 [11] criteria, the microbial growth rate expressions, defined in the model, were based on Monod’s kinetic. (Equation (4)) [20,21].
f ( S i ) = S i t = μ α i X i Y i ; f ( X i ) = X i t = μ X i K d X i μ = μ m a x i S i K s i + S i ;

2.4. Numerical Method, Finite Element

The Galerkin Finite Element Method (GFEM) was used to solve the governing equations defined in the model (Equations (1)–(3)). The procedure for converting the partial differential equations into a sum of simple polynomials was as follows:
1.
Discretization of the domain into unstructured meshes.
2.
The differential equation is first multiplied by a weight function w(x) and then integrated over the domain.
3.
Choose the order of interpolation (e.g., linear, quadratic, etc.) and corresponding shape functions N i = 1 m , with trial function;
p = p ˜ ( x ) = i = 1 m N i ( x ) p i .
4.
Evaluate all integrals over each element in order to set up a system of equations in the unknown p i s .
5.
Solve the system of equations for the p i s .

Tools

The calculation was performed using FreeFem++/FreeFem-cs (FF). This is a programming language designed to solve partial differential equations by the finite element method. One of the most important aspects to be taken into account in the design of the algorithm was the need for the calculation of all the biochemical reactions to be carried out simultaneously, as occurs in nature itself. For this purpose, the algorithm structure was based on parallel computation, contemplated in the FF message passing interface (MPI), to be executed on PCs with parallel computing architectures. The use of MPI in this methodology allows for improved large-scale computation, considering the linkages between different processes defined in the ADM1, maximizing CPU performance, and reducing computation time. The other software used for the calculation were GMSH v4.11, for the generation of finite element meshes, and PARAVIEW v5.11, for the post-processing. All of these, along with FreeFem++ v4.12, are free and open-source software.

2.5. Case Study

In this section, we have been considered an NSAR to be implemented in an isolated farm and with limited power supply resources. This semicontinuous batch reactor has cylindrical geometry: 3 m in diameter and 3 m high (Figure 2a). It has an attached dome from which the biogas collected inside is extracted through a gas outlet pipe located at the top of the dome. The influent is loaded into the reactor from the lower part of the vessel, while the effluent is unloaded from the upper part. This process takes place in an instant and causes an upward movement of the entire liquid. An inlet flow Q = 7.07 (m· s−1) and a retention time of 1 day have been considered.
The anaerobic degradation of sugars from glucose will be studied considering the three stages: hydrolysis, acidogenesis and methanogenesis. The physical and kinetic parameters are as shown in Table 1.

FreeFem++ Code

The FreeFem++ code for solving this problem is given in Appendix A. Five main blocks can be highlighted in the proposed algorithm.
1.
Definitions of parameters, boundaries and finite element spaces.
2.
Description of macros and functions.
3.
Problem definition.
4.
Loading the loop and calculation of different processes.
5.
Saving the results for manipulation in the post-processes
Figure 2b shows the mesh generated by the discretization with a total of 401,133 tetrahedrons, 30,118 triangles, and 431,252 elements. The 3D mesh is loaded with the command “gmshload3”.
The description of the reactor feed will be made in a vertical and upward direction. Thus, an attempt has been made to reproduce an instantaneous loading into the reactor of the influent. Therefore, the loop in which the calculation is described is outside the time path defined in the flow diagram. The Stokes equation in the weak formulation is as follows;
μ Ω u · v Ω ( d i v v ) p = Ω f v f o r a l l v ϵ H 0 1 ( Ω ) 3
Ω ( d i v u ) q = 0 f o r a l l q ϵ L 3 ( Ω ) ,
they are accompanied by the following boundary conditions.
u ( x , y , z ) = h D ( x , y , z ) f o r x , y , z ϵ Γ D u · n = h N ( x , y , z ) f o r x , y , z ϵ Γ N
The weak formulations and boundary conditions for ADE (Equation (6)) and DRE (Equation (7)) are shown below.
Ω ϕ t · Y i + D Ω ϕ · Y i + Ω u ϕ · Y i = F ( x i ) · Y i f o r a l l Y i ϵ H 0 1 ( Ω ) ϕ ( x i , t ) = g D ( x i ) f o r x i ϵ Γ D Ω , t > 0 ϕ ( x i , t ) n n ( x i ) = g N ( x i ) f o r x i ϵ Γ N Ω , t > 0
Ω ϕ t · Y i + D Ω ϕ · Y i + Ω f ( ϕ ) · Y i = F ( x i ) · Y i f o r a l l Y i ϵ H 0 1 ( Ω ) ϕ ( x i , t ) = g D ( x i ) f o r x i ϵ Γ D Ω , t > 0 ϕ ( x i , t ) n n ( x i ) = g N ( x i ) f o r x i ϵ Γ N Ω , t > 0
These equations are described through macros that are applied to both the substrates and the microorganisms for each of the processes: acidogenesis, acetogenesis and methanogenesis.

3. Results and Discussion

This section contains the results obtained for the example proposed in Section 2.5.
Figure 3 shows the steady-state velocity field as a function of the established boundary conditions. An upward flow is observable where the values are higher in the central zone, reaching values of u = 0.3 (m−1 s−1), as opposed to the extremes, in the upper zone with values close to zero. Since we are in the presence of a viscous fluid, the frictional forces between the fluid and the reactor walls force the velocities to be zero at these points—hence the tendency of the upward flow toward the central axis (Figure 3).
The velocity field map shown will change as soon as the boundary conditions are modified, both in input values and in the location of inputs and outputs.
Figure 4 shows a snapshot of the simulation, t = 19 days, for propionic acid acetogenesis. The substrate concentrations are projected onto a vertical plane that cuts the vessel in half (Figure 4a). Here, it is noted how the maximum values are located at the bottom of the reactor. There is also a narrow band in the central zone, forming an arc, which reflects higher concentrations in contrast to the surrounding zone.
Figure 4a shows the distributions of substrate and microorganisms for the acetogenesis of propionic acid at 19 days. The first, plotted on the cylindrical surface of the reactor vessel, illustrates the variation of concentrations along the flow. This is understood when the concentrations of microorganisms are added to the graph (Figure 4b). The growth of biomass in the lower zone leads to a reduction of substrate.
By studying the two graphs through overlapping, it is possible to analyze the system in detail. This facilitates the development of skills for the acquisition of knowledge from experience.
Figure 5 shows the concentration map for acidogenesis, acetogenesis and methanogenesis, at instant t = 26 days. In acidogenesis (Figure 5a), the activity is very intense. Both substrate and microorganism concentrations are localized in the lower part of the vessel. These are supplied at the moment of digester feeding. However, they do not transcend to higher levels as cell metabolism occurs rapidly. In the case of acetogenesis and methanogenesis (Figure 5b,c), the substrates, the product of the metabolism of the previous stage, are generated at higher levels. The acid-degrading bacteria, propionic and acetic, will have to wait until they reach the higher levels where the substrate sources are located.
Thanks to the potential of the provided algorithm, it is feasible to improve the biomass yield. For this purpose, the student can modify the input data in the code to establish different strategies related to the workings of the reactor, such as the extension of the residence time or the reduction of the input flow rate.
Figure 6 shows the projections of the cell concentration distributions in the xz and xy planes for t = 26 (d).
With this methodology, the student can both examine and plan his or her own work as well as identify the successes and difficulties of the results. It is also possible to use different study strategies depending on the situation, varying the boundary conditions and parameters related to the micro-organisms and substrates.
Figure 6 illustrates the biomass concentration map in two planes, horizontal and vertical. The maximum values are observed to be located on the central axis and on the lowest part of the reactor vessel.
Figure 6 is a plot of the projections, in the xz and xy planes, of the cell concentrations. It can be seen how the maximum values are concentrated in the central axis and in the lowest part of the reactor vessel. With this type of representation, it is possible to closely analyze any part of the system, thus becoming a tool to generate and enhance an enriching learning environment based on all the information obtained in the simulation.
Line graphs can be applied with this tool. In Figure 7, a line graph has been plotted along the axis of the cylinder representing the reactor vessel (a). The diagrams presented in this figure show the concentrations of substrates in COD along the plotted axis, the abscissa axis, for sugar (a), propionic acid (b), and acetic acid (c). The results indicate, as in the graphs in the previous figures, that the highest microbial activity occurs from the reactor floor to the 1.5 m level. Above this level, residual COD is observed. The COD of propionic acid produced in the previous step, acidogenesis, is presumably higher than that consumed in acetogenesis itself (c). In the case of acetic acid, the concentrations generated in acidogenesis and acetogenesis relative to those consumed in methanogenesis are higher in the upper part of the vessel (d). It is likely that due to the decrease in biomass as a consequence of substrate reduction along the flow pathway, and also due to the time lag between substrate production and consumption, organic growth rates are likely to be increased at the top of the reactor. Although the linear graph shown has been plotted at a given instant, the provided information can also be interpreted as the behavior of the biomass along the entire path from the bottom of the vessel to the top where the effluent outlet is located.
In this case, there is consistency with the graphs seen previously. The progression of the process, reflected in the entire upward trajectory of the flow, describes a growth of biomass and a reduction of substrates. From very small initial concentrations, Xi = 0.05 (kg·m−2 COD), the exponential growth of biomass requires some time before substrate removal can be clearly perceived.

4. Conclusions

This paper presents a 3D code based on CFD to be implemented in virtual laboratories for studying anaerobic metabolism inside NSARs. The inclusion of NSARs studies in the training of engineers is justified since the application of these natural systems, suitable for isolated livestock facilities in the treatment of their waste, contributes to the mitigation of climate change and the strengthening of the circular economy.
The problem described consisted of modeling a semi-batch reactor with heterogeneous mixing; i.e., the tank is not equipped with mixers. Therefore, the challenge of this work was to develop a mathematical model and a resolution methodology for distributed parameter systems whose variables are defined in terms of space and time. Governing equations were based on the Stokes, advection–diffusion and diffusion–reaction equations. Even the structured model ADM1 was implemented for the description of anaerobic metabolism. The finite element method was used to solve these equations via the open access FF software.
As a result of this innovative proposal, an original methodology has been developed for inclusion in the educational community, for the active/interactive participation of students in training activities of great value for meaningful learning, within the university environment. Through the code provided, it is feasible that students can adapt it to the different examples that may arise, allowing them to complete the experiments through virtual laboratory simulations.

Author Contributions

Conceptualization, S.B.-E. and C.M.-P.; Methodology, S.B.-E. and A.R.-M.; Formal analysis, S.B.-E. and A.R.-M.; Investigation, S.B.-E., T.G.-R. and A.R.-M.; Writing—original draft, S.B.-E.; Writing—review & editing, S.B.-E.; Supervision, J.J.S. and F.L.-Z.; Project administration, J.J.S. and F.L.-Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been partially funded by the The Next Generation EU (NGEU) fund under “Real Decreto 641/2021. de 27 de julio. por el que se regula la concesión directa de subvenciones a universidades públicas españolas para la modernización y digitalización del sistema universitario español en el marco del plan de recuperación. transformación y resiliencia (UNIDIGITAL)—Proyectos de Innovación Educativa para la Formación lnterdisciplinar (PIEFI)—Línea 3. Contenidos y programas de formación” in the scope of the Teaching Innovation Project “Aplicación de técnicas de aprendizaje autónomo y colaborativo para la mejora de los resultados del aprendizaje en entornos de simulación virtual en el ámbito de la ingeniería (PIE 2021-60)”. This research was co-funded as well by the INTERREG V-A Cooperation, Spain-Portugal MAC (Madeira-Azores-Canarias) 2014–2020 program, MITIMAC project (MAC2/1.1a/263).

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

α Stoichiometric parameters
Δ 2 x i ; x i = x , y , z
μ (d−1)Specific growth rate
μ m a x i (d−1)Maximum specific growth rate
x i ; x i = x , y , z
ν viscosity of the fluid
ϕ represent both the concentration of substrates and cells (kg·m−3)
fsource function. Reactive term
K S 1 (kg m−3)Saturation constant of the substrate
K d Cell decay rate
K s Monod substrate saturation constant
K d (d−1)Cell decay rate
S i (kg COD m−3)Substrate concentration
X i (kg COD m−3)Cells concentration
Y 1 Substrate yield coefficient
D Diffusion coefficient
μ m a x Maximum specific growth rate
u velocity vector
Fexternal force applied to the fluid
ppressure
S i e n t Boundary condition for the substrate
X i e n t Boundary condition for the microorganism

Appendix A. FreeFem++ Code

//      3D model (CFD) for a non-stirred anaerobic digester        
//       Autor :Saulo M. Brito                 
if ( mpisize != 3) {    // total number of processes 
cout << ” sorry, number of processors !=3” << endl;
exit(1);}
load ”medit”
load ”iovtk”
load ”gmsh”
load ”msh3”
load ”tetgen”
load ”MUMPS_FreeFem”
include ’’MPIGMRESmacro.idp’’
int[int] ffordert=1;
int[int] ffordert1=[1,1,1];
int worldrank;
mpiComm comm(mpiCommWorld, 0, 0);
int MPICommSize=mpiSize(comm);
int MPIRank=mpiRank(comm);
int intlet=6,sidedown=7,sideup=8,sky=9,inter=11;//labels
//load the mesh archive(import digester geometry)
mesh3 th1= gmshload3(”cilinder1.msh”);
int numvertices=th1.nv;
int numtriangle=th1.nt;
plot(th1, wait=true);
//definition of kinetic constants
real Numax1=2.85,Ks1=0.05,Kd1=0.9;     //for sugar
real Numax4=1.8,Ks4=1.145,Kd4=0.04;   //for propionic acid
real Numax8=3.9,Ks8=0.93,Kd8=0.037;   //for acetic acid
real XIent=0.05,SIent=100,Hcont= 1 × 10−7;  //values of reactor inlet concentrations (kg/m3 COD)
real Df=8.64 × 10−3;      //Diffusion coefficient m2/d
real Q=7.0686;      //flow (m3/s)
int conteo=0, cpu=clock();
//   time period   
real T=0.01,dt=1.5 × 102;      // T (d), dt(dias)
fespace Vh(th1,P13d);
fespace Nh(th1,[P13d,P13d,P13d,P13d]);
Nh [u1,u2,u3,pp1],[v1,v2,v3,pp2];
Vh dU,VH,U;
Vh EpH1,EpH2,EpH3,IpHant=0,pHll=6,pHul=8.5,pHcat,IpHcal;
//   definition of variables   
Vh si1,xi1,yi1,qi1,si2,xi2,yi2,qi2,si3,xi3,yi3,qi3,TCO2,TNH3=0.1;
Vh NuT,NuT1,NuT2,NuT3,TSiPac,TSiPb,TSiPva,TSiPp,IpH=0.6,INH3,KINH3=7.5;
real tmax1=45.8,tmin1=11,topt1=39.3,tmax2=47.3,tmin2=5.6,topt2=40.3,
tmax3=46.3,tmin3=11.1,topt3=34.1;
macro Grad(u) [dx(u),dy(u),dz(u)]//
macro div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3))//
macro grad(si) (dx(si)+dy(si)+dz(si))//
macro BAD(xi1,yi1) (int3d(th1)(1/dt*(xi1*yi1))
+int3d(th1)(Df/86400*(Grad(xi1)’*Grad(yi1)))  
+int3d(th1)(-1/dt*convect([u1,u2,u3],-dt,xiold[mpirank])*yi1)
) //
macro BDR(xi2,yi2) (int3d(th1)(1/dt*(xi2*yi2))
+int3d(th1)(Df*(Grad(xi2)’*Grad(yi2)))  
-int3d(th1)((Nu[mpirank]*siold[mpirank]/(Ks[mpirank]+siold[mpirank])*IpH- 
Kd[mpirank])*xi2*yi2)
-int3d(th1)(1/dt*xiold[mpirank]*yi2)
) //
macro SAD(si1,qi1) (int3d(th1)(1/dt*(si1*qi1))
+int3d(th1)(Df/86400*(Grad(si1)’*Grad(qi1)))  
+int3d(th1)(-1/dt*convect([u1,u2,u3],-dt,siold[mpirank])*qi1)
) //
macro SDR(si2,qi2) (int3d(th1)(1/dt*(si2*qi2))
+int3d(th1)(Df*(Grad(si2)’*Grad(qi2)))  
+int3d(th1)(Nu[mpirank]*siold[mpirank]/(Ks[mpirank]+siold[mpirank])*IpH* 
Alpha[mpirank]*xiold[mpirank]*qi2)
-int3d(th1)(1/dt*siold[mpirank]*qi2)
) //
func vel=Q/(6*7.0686); // function of velocity (Q/m2)[m/seg]
solve stokes([u1,u2,u3,pp1],[v1,v2,v3,pp2]) = int3d(th1,qforder=3)(Grad(u1)’*Grad(v1)
 +Grad(u2)’*Grad(v2)+Grad(u3)’*Grad(v3)
+0.000000001*pp1*pp2
- div(u1,u2,u3)*pp2 - div(v1,v2,v3)*pp1 + 1e-10*pp2*pp1 )
+on(sideup,u1=0,u2=0,u3=0) + on(intlet,inter,u1=0,u2=0,u3=vel);
savevtk("velocity.vtk",th1,[u1,u2,u3],order=ffordert1, dataname="ui");
//                                        
//********************Definition of the MPI workspace**********************
//                                        
Vh[int] Nu(mpisize),Ks(mpisize),siold(mpisize),xiold(mpisize),SiPac(mpisize),
SiPp(mpisize),SiPb(mpisize),SiPva(mpisize),SiPcl(mpisize),
CO2(mpisize),NH3(mpisize),sifont(mpisize),sifontb(mpisize),siV(mpisize),
siB(mpisize),CC1(mpisize),tco21
(mpisize),tnh31(mpisize),CC2(mpisize),tco22(mpisize),tnh32(mpisize);
real[int] Kd(mpisize),Alpha(mpisize),XiaciZ(mpisize),
SiaciZ(mpisize),XiaceZ(mpisize),SiaceZ(mpisize),XimethZ(mpisize),SimethZ
(mpisize),SIcont(mpisize),XIcont(mpisize);        // source coeficient
if (mpirank==0){
Kd[0]=Kd1;Alpha[0]=10.76;xiold[mpirank]=0;siold[mpirank]=0;      //SUGAR
SiPac[0];SiPp[0];CO2[0];NH3[0];Nu[mpirank]=Numax1;Ks[mpirank]=Ks1;};
if (mpirank==1){
Kd[mpirank]=Kd4;Alpha[mpirank]=11.29943503;xiold[mpirank]=0;
siold[mpirank]=0; // PROPIONIC ACID
SiPac[mpirank];Nu[mpirank]=Numax4;Ks[mpirank]=Ks4;sifont[mpirank];}
if (mpirank==2){
Kd[mpirank]=Kd8;Alpha[mpirank]=18.1812;xiold[mpirank]=0;
siold[mpirank]=0;  // ACETIC ACID
sifont[mpirank];Nu[mpirank]=Numax8;Ks[mpirank]=Ks8;} // end mpirank
//        Calculation             
int nn=T/dt+2,ii,in=1;  // loop
for(real tt=0; tt<=T;tt+=dt)
{ii=(tt+dt)/dt;
if(ii==in){                // reactor feeding
for(int mm=0;mm<=6; mm+=1)
{
if (mpirank==0){
solve Bich1(xi1,yi1, solver=GMRES)= BAD(xi1,yi1)
+on(intlet,xi1=XIent)+on(sideup,sidedown,xi1=0);
solve Sust1(si1,qi1, solver=GMRES)= SAD(si1,qi1)
+on(intlet,si1=SIent)+on(sideup,sidedown,si1=0);
xiold[mpirank]=xi1;siold[mpirank]=si1;
};
if (mpirank==1){
solve Bich2(xi2,yi2, solver=GMRES)= BAD(xi2,yi2)
+on(intlet,xi2=XIent)+on(sideup,sidedown,xi2=0);
solve Sust2(si2,qi2, solver=GMRES)= SAD(si2,qi2);
xiold[mpirank]=xi2;siold[mpirank]=si2;
};// end if(mpirank…)
if (mpirank==2){
solve Bich3(xi3,yi3, solver=GMRES)= BAD(xi3,yi3)
+on(intlet,xi3=XIent)+on(sideup,sidedown,xi3=0);
solve Sust3(si3,qi3, solver=GMRES)=SAD(si3,qi3);
xiold[mpirank]=xi3;siold[mpirank]=si3;
};// end if(mpirank…)
};// end for
in=in+33;
};// end if
if (mpirank==0){solve Bich1(xi1,yi1, solver=GMRES)= BDR(xi1,yi1);
solve Sust1(si1,qi1, solver=GMRES)= SDR(si1,qi1);
xiold[mpirank]=xi1;siold[mpirank]=si1;
// **************************product***************************************
SiPac[mpirank]=siold[mpirank]*2.669537137/10.76426265;SiPp[mpirank]=
siold[mpirank]*3.139935414/10.76426265;
//        Choice of the smallest value between Si y Xi       
for (int i=0;i<Vh.ndof;i+=1){
if(xiold[mpirank][](i)<=0.1e-9){xiold[mpirank][](i)=0;}
if(siold[mpirank][](i)>=SIent){siold[mpirank][](i)=SIent;}
if(siold[mpirank][](i)< 0.1e-9){siold[mpirank][](i)=0;SiPac[mpirank][](i)=0;
SiPp[mpirank][](i)=0; }
if(siold[mpirank][](i)/10.76426265>=xiold[mpirank][](i)){SiPac[mpirank][](i)=
2.669537137*xiold[mpirank][](i);
SiPp[mpirank][](i)=3.139935414*xiold[mpirank][](i);  ;}
}//end for
processor(1) << SiPp[mpirank][] ;processor(2) << SiPac[mpirank][];
savevtk("Acidsugar_xi"+ii+".vtk",th1,xiold[mpirank],order=ffordert, dataname="Xi_hidrol");
savevtk("Acidsugar_si"+ii+".vtk",th1,siold[mpirank],order=ffordert, dataname="Si_hidrol");
};// end if(mpirank…)
if (mpirank==1){processor(0) » SiPp[0][];
sifont[mpirank]= SiPp[0];
solve Bich2(xi2,yi2, solver=GMRES)= BDR(xi2,yi2);
solve Sust2(si2,qi2,solver=GMRES)= SDR(si2,qi2)
-int3d(th1)(sifont[mpirank]*qi2);
xiold[mpirank]=xi2;siold[mpirank]=si2;
// **************************product***************************************
SiPac[mpirank]=siold[mpirank]*6.0307/11.294; Choice of the smallest value between Si y Xi
for (int i=0;i<Vh.ndof;i+=1){
if(xiold[mpirank][](i)<=0.1e-9){xiold[mpirank][](i)=0;}
if(siold[mpirank][](i)< 0.1e-9){siold[mpirank][](i)=0;
SiPac[mpirank][](i)=0;SiPp[mpirank][](i)=0; }
if(siold[mpirank][](i)/11.294>=xiold[mpirank][](i)){SiPac[mpirank][](i)=
6.0307*xiold[mpirank][](i);};
}//end for
processor(2) << SiPac[mpirank][];
savevtk("Propionic_xi"+ii+".vtk",th1,xiold[mpirank],order=ffordert, dataname="Xi_propionic");
savevtk("Propionic_si"+ii+".vtk",th1,siold[mpirank],order=ffordert, dataname="Si_propionic");
};// end if(mpirank…)
if (mpirank==2){processor(0) » SiPac[0][];processor(1) » SiPac[1][];
sifont[mpirank]= SiPac[0]+SiPac[1];
solve Bich3(xi3,yi3, solver=GMRES)= BDR(xi3,yi3);
solve Sust3(si3,qi3, solver=GMRES)= SDR(si3,qi3)
-int3d(th1)(sifont[mpirank]*qi3);
xiold[mpirank]=xi3;siold[mpirank]=si3;
for (int i=0;i<Vh.ndof;i+=1){
if(xiold[mpirank][](i)<=0.1e-9){xiold[mpirank][](i)=0;}
if(siold[mpirank][](i)< 0.1e-9){siold[mpirank][](i)=0;}
}//end for
savevtk("Acetico_xi"+ii+".vtk",th1,xiold[mpirank],order=ffordert, dataname="Xi_acetic");
savevtk("Acetico_si"+ii+".vtk",th1,siold[mpirank],order=ffordert, dataname="Si_acetic");
};// end if(mpirank…)
};// end for

References

  1. States, N.; Stone, E.; Cole, R. Creating Meaningful Learning Opportunities through Incorporating Local Research into Chemistry Classroom Activities. Educ. Sci. 2023, 13, 192. [Google Scholar] [CrossRef]
  2. Camilleri, M.A.; Camilleri, A.C. Digital Learning Resources and Ubiquitous Technologies in Education. Technol. Knowl. Learn. 2017, 22, 65–82. [Google Scholar] [CrossRef]
  3. Liu, D.; Valdiviezo-Díaz, P.; Riofrio, G.; Sun, Y.M.; Barba, R. Integration of Virtual Labs into Science E-learning. Procedia Comput. Sci. 2015, 75, 95–102. [Google Scholar] [CrossRef] [Green Version]
  4. Watts, J.; Crippen, K.J.; Payne, C.; Imperial, L.; Veige, M. The varied experience of undergraduate students during the transition to mandatory online chem lab during the initial lockdown of the COVID-19 pandemic. Discip. Interdiscip. Sci. Educ. Res. 2022, 4, 14. [Google Scholar] [CrossRef]
  5. Appels, L.; Lauwers, J.; Degrève, J.; Helsen, L.; Lievens, B.; Willems, K.; Van Impe, J.; Dewil, R. Anaerobic digestion in global bio-energy production: Potential and research challenges. Renew. Sustain. Energy Rev. 2011, 15, 4295–4301. [Google Scholar] [CrossRef]
  6. Guzmán, J.L.; Joseph, B. Web-Based Virtual Lab for Learning Design, Operation, Control, and Optimization of an Anaerobic Digestion Process. J. Sci. Educ. Technol. 2021, 30, 319–330. [Google Scholar] [CrossRef]
  7. Husain, A. Mathematical models of the kinetics of anaerobic digestion—A selected review. Biomass Bioenergy 1998, 14, 561–571. [Google Scholar] [CrossRef]
  8. Wade, M.J. Not Just Numbers: Mathematical Modelling and Its Contribution to Anaerobic Digestion Processes. Processes 2020, 8, 888. [Google Scholar] [CrossRef]
  9. Batstone, D.J.; Puyol, D.; Flores-Alsina, X.; Rodríguez, J. Mathematical modelling of anaerobic digestion processes: Applications and future needs. Rev. Environ. Sci. Bio/Technol. 2015, 14, 595–613. [Google Scholar] [CrossRef]
  10. Borisov, M.; Dimitrova, N.; Simeonov, I. Mathematical Modelling of Anaerobic Digestion with Hydrogen and Methane Production. IFAC PapersOnLine 2016, 49, 231–238. [Google Scholar] [CrossRef]
  11. Batstone, D.; Keller, J.; Angelidaki, I.; Kalyuzhnyi, S.; Pavlostathis, S.; Rozzi, A.; Sanders, W.; Siegrist, H.; Vavilin, V. The IWA Anaerobic Digestion Model No 1 (ADM1). Water Sci. Technol. 2002, 45, 65–73. [Google Scholar] [CrossRef] [PubMed]
  12. Kleerebezem. Critical analysis of some concepts proposed in ADM1. Water Sci. Technol. 2006, 54, 51–57. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Li, L.; Wang, K.; Zhao, Q.; Gao, Q.; Zhou, H.; Jiang, J.; Mei, W. A critical review of experimental and CFD techniques to characterize the mixing performance of anaerobic digesters for biogas production. Rev. Environ. Sci. Bio/Technol. 2022, 21, 665–689. [Google Scholar] [CrossRef]
  14. Hoff, C.; Aurandt, J.; O’Toole, M.; Davis, G. Motivating Student Learning Using Biofuel-Based Activities. In Proceedings of the 2013 ASEE Annual Conference & Exposition (Online), Atlanta, GA, USA, 22–26 June 2020; Available online: https://peer.asee.org/19019 (accessed on 24 January 2023).
  15. Nanou, P.; Kuchta, D.I.K. Knowledge Transfer Networks in the Field of Biological Waste Treatment: Current Status; Literature Review; Hamburg University of Technology: Hamburg, Germany, 2015. [Google Scholar]
  16. Sadino-Riquelme, C.; Hayes, R.E.; Jeison, D.; Donoso-Bravo, A. Computational fluid dynamic (CFD) modelling in anaerobic digestion: General application and recent advances. Crit. Rev. Environ. Sci. Technol. 2018, 48, 39–76. [Google Scholar] [CrossRef]
  17. Song, L.; Li, P.W.; Gu, Y.; Fan, C.M. Generalized finite difference method for solving stationary 2D and 3D Stokes equations with a mixed boundary condition. Comput. Math. Appl. 2020, 80, 1726–1743. [Google Scholar] [CrossRef]
  18. Ukai, S. A solution formula for the Stokes equation in Rn+. Commun. Pur. Appl. Math. 1987, 40, 611–621. [Google Scholar] [CrossRef] [Green Version]
  19. Parshotam, A.; Mcnabb, A.; Wake, G. Comparison Theorems for Coupled Reaction-Diffusion Equations in Chemical Reactor Analysis. J. Math. Anal. Appl. 1993, 178, 196–220. [Google Scholar] [CrossRef] [Green Version]
  20. Muloiwa, M.; Nyende-Byakika, S.; Dinka, M. Comparison of unstructured kinetic bacterial growth models. S. Afr. J. Chem. Eng. 2020, 33, 141–150. [Google Scholar] [CrossRef]
  21. Monod, J. The Growth of Bacterial Cultures. Annu. Rev. Microbiol. 1949, 3, 371–394. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Red: Neumann boundary value problem; Blue: Dirichlet boundary value problem; boundary = Ω = Γ 1 Γ 2 ; ϕ and g are given function defined on Γ 1 and Γ 2 ; n = unit normal vector.
Figure 1. Red: Neumann boundary value problem; Blue: Dirichlet boundary value problem; boundary = Ω = Γ 1 Γ 2 ; ϕ and g are given function defined on Γ 1 and Γ 2 ; n = unit normal vector.
Water 15 01181 g001
Figure 2. (a) Reactor geometry. (b) Discretization of the volume of the liquid portion.
Figure 2. (a) Reactor geometry. (b) Discretization of the volume of the liquid portion.
Water 15 01181 g002
Figure 3. Perspective and top views of velocity field. The arrows represent the projection of the normalized velocity vectors, i.e., the behavior of the non-normal velocity components. The magnitude of the velocity vector is indicated by the size of the arrow and the intensity of the color in the key of the figure.
Figure 3. Perspective and top views of velocity field. The arrows represent the projection of the normalized velocity vectors, i.e., the behavior of the non-normal velocity components. The magnitude of the velocity vector is indicated by the size of the arrow and the intensity of the color in the key of the figure.
Water 15 01181 g003
Figure 4. Substrate (a) and cells (b) at time t = 19 days. The keys at the right of the figure indicate microbial and substrates concentrations in kg (COD).
Figure 4. Substrate (a) and cells (b) at time t = 19 days. The keys at the right of the figure indicate microbial and substrates concentrations in kg (COD).
Water 15 01181 g004
Figure 5. COD distribution at 26 days for biomass (top) and substrates (bottom): (a) acidogenesis, (b) acetogenesis, and (c) methanogenesis. The keys at the right of the figure indicate microbial and substrates concentrations in kg (COD).
Figure 5. COD distribution at 26 days for biomass (top) and substrates (bottom): (a) acidogenesis, (b) acetogenesis, and (c) methanogenesis. The keys at the right of the figure indicate microbial and substrates concentrations in kg (COD).
Water 15 01181 g005
Figure 6. Representation of cell concentrations in the planes: (a) longitudinal, parallel to the xz plane and cutting the cylinder into two equal parts, and (b) transverse, parallel to the axis of the xy plane and located at a vertical height of 0.8 m.
Figure 6. Representation of cell concentrations in the planes: (a) longitudinal, parallel to the xz plane and cutting the cylinder into two equal parts, and (b) transverse, parallel to the axis of the xy plane and located at a vertical height of 0.8 m.
Water 15 01181 g006
Figure 7. Evolution of COD concentrations in the NSAR along the y-axis (a), for the acidogenesis of sugars (b), the acetogenesis of propionic acid (c) and the methanogenesis of acetic acid (d).
Figure 7. Evolution of COD concentrations in the NSAR along the y-axis (a), for the acidogenesis of sugars (b), the acetogenesis of propionic acid (c) and the methanogenesis of acetic acid (d).
Water 15 01181 g007
Table 1. Physical and kinetic parameters for: 1 sugar, 2 propionic acid, 3 acetic acid.
Table 1. Physical and kinetic parameters for: 1 sugar, 2 propionic acid, 3 acetic acid.
Process μ max D Si ent Xi ent K d K s α
d−1m2·d−1kg· m−3kg· m−3d−1kg· m−3kg · kg−1
12.85 8 · 10 3 1000.050.90.0510.76
21.8 8 · 10 3 -0.050.041.14511.29
33.9 8 · 10 3 -0.050.0370.9318.18
Note: D —Diffusion coefficient. μ max —Maximum specific growth rate. Si ent —Boundary condition for the substrate. Xi ent —Boundary condition for the microorganism. K s —Monod substrate saturation constant. K d —Cell decay rate. α —Stoichiometric parameters.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Brito-Espino, S.; García-Ramírez, T.; Leon-Zerpa, F.; Mendieta-Pino, C.; Santana, J.J.; Ramos-Martín, A. E-Learning Proposal for 3D Modeling and Numerical Simulation with FreeFem++ for the Study of the Discontinuous Dynamics of Biological and Anaerobic Digesters. Water 2023, 15, 1181. https://doi.org/10.3390/w15061181

AMA Style

Brito-Espino S, García-Ramírez T, Leon-Zerpa F, Mendieta-Pino C, Santana JJ, Ramos-Martín A. E-Learning Proposal for 3D Modeling and Numerical Simulation with FreeFem++ for the Study of the Discontinuous Dynamics of Biological and Anaerobic Digesters. Water. 2023; 15(6):1181. https://doi.org/10.3390/w15061181

Chicago/Turabian Style

Brito-Espino, Saulo, Tania García-Ramírez, Federico Leon-Zerpa, Carlos Mendieta-Pino, Juan J. Santana, and Alejandro Ramos-Martín. 2023. "E-Learning Proposal for 3D Modeling and Numerical Simulation with FreeFem++ for the Study of the Discontinuous Dynamics of Biological and Anaerobic Digesters" Water 15, no. 6: 1181. https://doi.org/10.3390/w15061181

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop