While computational modelling has increasingly supported wastewater bioreactor engineering, novel numerical techniques have been developed such as the lattice-Boltzmann method (LBM). With vinasse treatment as case study, this work is a first step towards a comprehensive LBM simulator of a continuous-flow anaerobic packed-bed reactor. Extensions from typical models comprise one-dimensional (besides time) dependence, species transport via convection and diffusion, and imposition of either Dirichlet or Danckwerts condition at inlet. The LBM simulator proved to be operational when simulating the bioreactor at different hydraulic retention times (HRTs). Simulated profiles show that stepwise feeding concentrations are smoothed as they are transported towards the bioreactor exit while concentrations increase or decrease in response to generation or degradation kinetics. Good fitting was observed for concentrations of acetic acid (2.1 kg-COD/m3 for HRT = 24 h) and butyric acid (1.3 kg-COD/m3 for HRT = 16 h) at the exit whereas other concentrations were numerically simulated at proper order of magnitude.
Computational modelling can provide insights for both fundamental and applied research, particularly when acquiring real-operation or experimental data becomes awkward or even unfeasible (de Souza-Santos 2010). Virtualisation allows one to rapidly examine exploratory scenarios and search for optimal conditions. Accordingly, the shift from empirical (i.e. black-box or data-driven) to physics-based (i.e. white-box or phenomenological) models is vital to benefit from new sensor technology, internet of things and big data (Saguy 2016).
Comprehensive knowledge of bioreactors and wastewater treatment may rely on equations invoking numerical solution methods (Głuszcz et al. 2011; Naessens et al. 2012). Besides variability in chemical kinetics and composition, mathematical hurdles refer to the mutual interference between fluid velocity and chemical species concentrations (Parco et al. 2007). Computational fluid dynamics (CFD) arises as a helpful tool (Brannock et al. 2010) and the importance of computational modelling to bioprocess engineering has been recognised while novel numerical methods have been developed and applied (Datta & Sablani 2007).
The lattice-Boltzmann method (LBM) was envisaged in McNamara & Zanetti (1988) as an offspring of another mesoscale method, namely lattice gas cellular automata. LBM has become a powerful technique to numerically simulate bioprocesses and biosystems at different length scales (van der Sman 2007), including bioreactor design (Štumbauer et al. 2013). Within the scope of bioreactor simulation, one may apply LBM to computationally model biofilm formation and detachment (Picioreanu et al. 2001) as well as biogas bubbles generation and transport (Chen 2010).
From a modelling standpoint, LBM does not directly solve Navier–Stokes equations (Succi 2001) so that it can quite straightforwardly deal with fluid flow, transport phenomena, moving boundaries and multicomponent systems, whether or not coupled to chemical reactions (Sukop & Thorne 2006). On the computational side, LBM is suitable for parallel processing while it renders relatively simpler codes (Mohamad 2011). Those features sound attractive when programming in-house CFD simulators via classic discretisation methods such as finite differences (FDM), finite elements (FEM) or finite volumes (FVM)
LBM avoids the numerical stiffness related to either Galerkin formulations in FEM (Bathe 1996) or pressure-velocity coupling in FVM (Patankar 1980). Discretisation of convective terms (e.g. upwind scheme) may yield false (i.e. numerical) diffusion (Patankar 1980), which is not the case of LBM (van der Sman 2007; Mohamad 2011). False diffusion effects arise in FDM simulations of convective-dominant transport phenomena while they are mitigated by means of LBM (Rabi 2012).
Inspired by the IWA Anaerobic Digestion Model No. 1 (ADM1) (Batstone et al. 2002), this work is a first step in the development of a comprehensive in-house LBM simulator of continuous-flow bioreactors for wastewater treatment. While there is no doubt about the efficiency of classic numerical methods and off-the-shelf software, the present work points to LBM as an alternative route to computationally model bioreactors for wastewater treatment.
The physics-based model proposed in this work extends the habitual time-dependent (i.e. dynamic) behaviour to include one-dimensional (1-D) spatial dependence. As far as species transport is concerned, the model considers convection as well as diffusion. The latter has been sometimes ignored as convective-dominant differential equations remain parabolic with respect to spatial dependence (Wu & Hou 2001; Sanchez et al. 2004). Consequently, one can dismiss one boundary condition (usually at the exit) while being able to use marching numerical methods such as the Runge–Kutta method (Kreyszig 1993).
Although governing differential equations become elliptic, diffusive transport should be modelled not only when comprehensive knowledge is necessary but also because it is influential in large equipment (Reis-Vasco et al. 2000; Gaspar et al. 2003). In view of that, convective–diffusive transport has been considered since early versions of our in-house LBM simulators (Durán et al. 2015; Okiyama et al. 2015; Rabi & Kamimura 2016; Rosa et al. 2016). In the present work, the null Neumann condition was imposed at the bioreactor exit while either Dirichlet or Danckwerts boundary conditions were numerically implemented and tested at the bioreactor inlet.
Physics-based model of continuous-flow bioreactor for vinasse treatment
A cylindrical laboratory-scale anaerobic packed-bed reactor (APBR) for sugarcane vinasse treatment has been taken as case study (Ferraz 2013; Ferraz et al. 2014). Vinasse is an effluent from the sugar-ethanol industry whose large-scale use has pointed to ferti-irrigation in sugarcane crops after it has undergone anaerobic treatment (Vlissidis & Zouboulis 1993). Figure 1(a) depicts the experimental APBR assembled within acrylic tubes. Each APBR comprises a feeding module (FM), bed section (BS), effluent collection (EC) area and biogas collection (BC) area, as shown in Figure 1(b).
Low-density polyethylene cylinders (approximate dimensions: length ∼5.0 mm, diameter ∼4.5 mm) were used as supporting medium in the BS module. Those small cylinders randomly filled up the BS module so that bed porosity ɛ was observed to range from 0.47 to 0.54 (Ferraz 2013; Ferraz et al. 2014). Table 1 summarises length, inner diameter, geometric volume and volume occupied by the liquid phase at each APBR compartment. As the BC module is solely occupied by the gas phase, the volume occupied by the liquid phase could range from Vliq,min = 2.187 L (for lower porosity) to Vliq,max = 2.363 L (for higher porosity).
APBR is thermophilic (at 55 °C) and four different volumetric flow rates of vinasse were tested in Ferraz (2013) and Ferraz et al. (2014). Flow rates were experimentally set so as to yield hydraulic retention times (HRTs) of 24 h, 16 h, 12 h and 8 h, assessed as . In numerical simulations, corresponding interstitial fluid velocities were evaluated in terms of lower porosity ɛmin and bed diameter d as . Measured data available for comparisons with numerical simulations (ahead in ‘Results and discussion’ section) refer to species concentrations in the EC outflow, Figure 1(b).
This work proposes a dynamic 1-D model in which coordinate axis z is oriented upwards along with vinasse flow. As sketched in Figure 1(c), stratification is assumed so that species concentrations yi are functions of time t and coordinate z, i.e. yi = yi(z,t). The 1-D solution domain spans from z = 0 (inlet) to z = L (exit), where length L is numerically set to yield the same HRT from experiments, i.e. L = v×HRT. Table 2 shows the values of APBR operational parameters HRT, , v and L used in LBM simulations.
Comprehensive versions of the LBM simulator intend to invoke as many chemical species as those listed in Table 3. For the sake of successful implementation of a preliminary version, this work numerically solves Equation (1) for selected species, namely y1 to y5 in Table 3 together with single (say ‘unified’) biomass species y7 (i.e. rather than invoking species y71 to y74).
Underlying concepts of LBM
Inspired by the kinetic gas theory, LBM considers any medium (whether solid or fluid) as comprised of fictitious constituent particles following a sequence of synchronised streaming and collision steps in a discrete space (namely a fictitious lattice). During streaming, particles displace from one site to another through links defined by the fictitious lattice assigned to the true medium. As particles simultaneously arrive at lattice sites, they mutually collide so that their discrete velocities are rearranged for subsequent streaming–collision steps. By imposing conservation principles (e.g. mass and momentum) to such repetitive particle dynamics, one may numerically simulate macroscopic medium behaviour (Succi 2001).
The basic mathematical entity in LBM is the particle distribution function giving, at time t and about position , the number of particles per unit volume with velocities between and . By taking suitable moments of function f, one is able to retrieve macroscopic (i.e. observable) properties such as chemical species concentration, bulk flow velocity, and temperature (Succi 2001; van der Sman 2007; Mohamad 2011).
LBM is implemented to numerically solve Equation (2) as written in terms of the fictitious lattice, when it becomes known as the lattice-Boltzmann equation (LBE). Particle distribution functions fk are then assigned to each streaming velocity ck in the lattice. LBM lattices are identified as DnQm, where n is the problem dimension (e.g. n = 1 for 1-D problem) while m refers to the number of functions fk to be solved for each macroscopic property in the true medium.
The present work proposes a dynamic 1-D model and Figure 2 sketches the basic (and repetitive) linear structure of (a) D1Q2 and (b) D1Q3 lattices. It comprises a central lattice site linked to two neighbouring sites, one in forward (k = 1) and the other in backward (k = 2) direction, thus leading to streaming velocities c1 = +c and c2 = −c. The difference between D1Q2 and D1Q3 relies on the fact that the latter lattice allows some particles to remain at rest (i.e. null velocity c0 = 0) after collisions. Usual 2-D and 3-D lattices in LBM are sketched elsewhere (Succi 2001; van der Sman 2007; Mohamad 2011).
NUMERICAL METHOD: LBM SIMULATION OF THE CONTINUOUS-FLOW APBR
As far as LBM simulation of species transport in multicomponent systems is concerned, one may follow two approaches. Species can be simulated either (i) as an active component whose local concentration increase leads to local decrease in background fluid concentration or (ii) as a passive solute carried by the flowing solvent (Sukop & Thorne 2006). In view of LBM simulators already implemented in our research, the present work follows the second approach.
As D1Q2 lattice is used, streaming was implemented in forward (Δz1 = +Δz) and backward (Δz2 = −Δz) directions. Streaming can be suppressed (with no loss of functionality) for any modelled property explicitly lacking derivatives with respect to spatial coordinates in its governing equation (Durán et al. 2015; Okiyama et al. 2015; Rabi & Kamimura 2016; Rosa et al. 2016).
RESULTS AND DISCUSSION
For the chemical species preliminary invoked, Table 4 also shows initial conditions yi(z,0) for Equation (12) as well as feeding values yi,in for either Danckwerts or Dirichlet inlet conditions, Equation (9) or (10). The latter were measured in Ferraz (2013) and Ferraz et al. (2014) and they refer to typical features of sugarcane vinasse. Null Neumann boundary condition, Equation (11), was imposed to all species at exit.
Besides the operational parameters from Table 2, LBM simulations were performed by setting kd1 = 0.02 day−1, Ysu = 0.1, fC2,su = 0.41, fC3,su = 0.27, fC4,su = 0.13, fH2,su = 0.19, and Di = 0.02 m2/day. The latter value was a small diffusivity arbitrarily set for all chemical species and a sensitivity analysis is left for future work. LBM parameters were set as Δz = 0.01 m and Δt = 0.0005 day in order to render small lattice-based Mach number (namely Ma = 0.1) for numerical stability purposes (Mohamad 2011).
Available experimental data refer to species concentrations y1 to y4 measured at EC exit, Figure 1(b), after 1 day of operation (Ferraz et al. 2014). For each operation (i.e. HRT) scenario in Table 2, species concentrations numerically simulated at z = L and t = 1 day are compared with aforesaid experimental data in Table 5. Good fitting was observed for some species concentrations (e.g. acetic acid for HRT = 24 h and butyric acid for HRT = 16 h) and other concentrations were numerically simulated at appropriate order of magnitude.
Differences between experimental values and LBM simulations might refer to the fact that chemical parameters have not been fine-tuned. Also, the model must be extended to include further chemical species and corresponding governing equations (together with additional parameters). Those computational modelling tasks are among future developments of the LBM simulator and they fall beyond the scope of this preliminary work.
Species concentrations in Table 5 are in fact a small sample of the numerical results one is able to obtain. The LBM simulator can indeed provide species concentrations yi at any axial position z and time t. While those numerical values can be contrasted with measured data, comparisons in this work are limited to experimental values available in Table 5. Yet, one may grasp bioreactor behaviour from simulations and, in what follows, concentration profiles are shown for LBM simulations performed for APBR operated at HRT = 12 h.
For each chemical species invoked in the current preliminary model, Figure 3 shows time-dependent concentrations yi(L,t) numerically simulated at the exit. Concentrations simulated under Danckwerts condition at the inlet, Equation (9), are compared to Dirichlet counterparts, Equation (10). Only minor differences are observed when one inlet boundary condition is replaced by the other. In effect, differences are prone to become evident for relatively high diffusivities Di (Rabi & Kamimura 2016). LBM simulations suggest that a time lag of approximately 8 h (∼1/3 day) is necessary for concentration changes to become traceable at the bioreactor exit.
While Figure 3 shows how species concentrations may vary in time at a specific position, the LBM simulator can also depict how concentrations may either spread or cluster within the bioreactor during operation. Accordingly, Figure 4 shows species concentration profiles yi(z,t) numerically simulated over the solution domain (i.e. for 0 ≤ z ≤ L) at the following instants: t = 0.25 day (=tfinal/4), t = 0.5 day (=tfinal/2), t = 0.75 day (=3tfinal/4), and t = 1.0 day (=tfinal).
From COD profiles simulated at t = 0.25 day and t = 0.50 day in Figure 4(a), one may note that the stepwise feeding concentration is smoothed while it is transported towards the bioreactor exit. As a result, there is a concentration increase at this position, which is evidenced in Figure 3(a) for t < 0.6 day. Subsequently, COD concentration at z = L is reduced due to its degradation throughout the bioreactor.
As far as acid concentrations are concerned, Figure 4(b)–4(d) suggest that stepwise feeding fronts are equally smoothed as they are transported to the exit. Conversely, those species concentrations are augmented due to acidogenesis over the bioreactor. Similar increasing behaviour was numerically simulated for dissolved H2 and biomass as Figure 4(e) and 4(f) respectively suggest, with degradation effects observed in the latter species.
All species concentration profiles in Figure 4 were numerically simulated by imposing null Neumann condition at the exit and Danckwerts condition at the inlet. LBM simulations under Dirichlet inlet condition yielded similar concentration profiles over the solution domain and they are not shown here for brevity. Yet, in Figure 4(f) one may spot the effect of imposing Danckwerts condition as the inlet concentration y7(0,t) slightly increases over time. This occurs because Danckwerts condition comes from a mass balance accounting for species transport by both convection and diffusion at the inlet (Danckwerts 1958). Conversely, in the case of imposing Dirichlet condition at inlet, concentration y7(0,t) remains fixed at the feeding value, i.e. y7,in = 0.1 kg-COD/m3 (as given in Table 4).
Short-term developments of the LBM simulator will refer to (i) including further species in Table 3 (together with governing equations) and (ii) coupling optimisation routines to best-fit model parameters (e.g. fine-tuning diffusivities Di) against available experimental data. Long-term improvements will include LBM simulation of fluid (i.e. vinasse) flow while extending the solution domain to two (and eventually three) dimensions.
Advances in dynamic modelling of bioreactors were successfully implemented in the LBM simulator. Specifically, this work extends bioreactor simulation to include 1-D spatial (besides time) dependence, species transport through convection as well as diffusion, and imposition of either Dirichlet or Danckwerts condition at the inlet. At its preliminary development stage, the LBM simulator proved to be operational when simulating the continuous-flow APBR at different HRT.
There was good fitting of concentrations of acetic acid (2.1 kg-COD/m3 for HRT = 24 h) and butyric acid (1.3 kg-COD/m3 for HRT = 16 h) while other species concentrations were numerically simulated at proper order of magnitude. LBM simulations were able to depict in detail how species concentrations change throughout the bioreactor length at specific time instants as well as how they vary in time at the bioreactor exit. Those numerical results provide insights to how chemical species may either spread or cluster within the bioreactor.