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

: 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 ﬁnite element method has been used to solve the partial differential equations deﬁned 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.


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.

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.

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;

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.

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].

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.

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].

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.

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.

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. Table 1. Physical and kinetic parameters for: 1 sugar, 2 propionic acid, 3 acetic acid.

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.
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; they are accompanied by the following boundary conditions.
The weak formulations and boundary conditions for ADE (Equation (6)) and DRE (Equation (7)) are shown below.
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.

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.

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. 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.