Skip to content

Advertisement

  • Research
  • Open Access

Advanced integration of fluid dynamics and photosynthetic reaction kinetics for microalgae culture systems

BMC Systems Biology201812 (Suppl 5) :93

https://doi.org/10.1186/s12918-018-0611-9

  • Published:

Abstract

Background

Photosynthetic microalgae have been in the spotlight of biotechnological production (biofuels, lipids, etc), however, current barriers in mass cultivation of microalgae are limiting its successful industrialization. Therefore, a mathematical model integrating both the biological and hydrodynamical parts of the cultivation process may improve our understanding of relevant phenomena, leading to further optimization of the microalgae cultivation.

Results

We introduce a unified multidisciplinary simulation tool for microalgae culture systems, particularly the photobioreactors. Our approach describes changes of cell growth determined by dynamics of heterogeneous environmental conditions such as irradiation and mixing of the culture. Presented framework consists of (i) a simplified model of microalgae growth in a culture system (the advection-diffusion-reaction system within a phenomenological model of photosynthesis and photoinhibition), (ii) the fluid dynamics (Navier-Stokes equations), and (iii) the irradiance field description (Beer-Lambert law). To validate the method, a simple case study leading to hydrodynamically induced fluctuating light conditions was chosen. The integration of computational fluid dynamics (ANSYS Fluent) revealed the inner property of the system, the flashing light enhancement phenomenon, known from experiments.

Conclusion

Our physically accurate model of microalgae culture naturally exhibits features of real system, can be applied to any geometry of microalgae mass cultivation and thus is suitable for biotechnological applications.

Keywords

  • Microalgae
  • Mathematical modeling
  • Photosynthesis
  • CFD
  • Microalgae culture systems
  • Flashing light enhancement

Background

The global warming, gradual oil fuel depletion, higher demands for energy and food consumptions are some of the current challenges our society is facing. Thus, looking for alternative resources is a hot topic. After the failure of the first generation biofuels based on corn and soya, which created a food shortage in the third world, the scientific community focused on simple photosynthetic organisms, cyanobacteria and microalgae. These organisms have high photosynthetic efficiency, biomass growth and lipids content [1].

Cyanobacteria and microalgae are very versatile microorganisms which have attracted a significant amount of scientific interest in the last decades [24]. These unicellular oxygenic organisms absorb photons of specific wavelengths to fix inorganic carbon (CO2). Carbon molecules are then utilized in the Calvin-Benson cycle, producing complex organic molecules in dark photosynthetic reactions and thus reducing the concentration of carbon dioxide in the environment. In addition, cyanobacteria and microalgae consume inorganic phosphorus and nitrogen which could be used in the processes of waste-water treatment [5]. The current scientific focus is to employ microalgae and cyanobacteria in biotechnology field as a potential biofuel source, for example a genetic modification increasing lipid content [6]. The other applications, e.g., pharmaceutical (bioactive metabolites), nutritional (animal feed and human food and supplements), are less important compared to the renewable energy business.

However, there are technological and knowledge-based barriers preventing optimized mass cultivation of photosynthetic microorganisms. In particular, reliable approaches allowing qualitative predictions for simulation of microbial biomass production in microalgae culture systems (MCS), including the photobioreactors (PBR), are rather slowly emerging than being well established [79]. The cause of apparent lack of accurate models of microalgae growth is in the complexity of three-dimensional multiphase (liquid-gas-solid) flow dynamics, irradiance distribution within MCS and multi-level processes of cellular functions. A reliable model has to integrate these phenomena which interact across various time-scales. Therefore, it is necessary to solve theoretical (multi-time-scale phenomena) and practical (high computational requirements) problems. As a result, the complex interactions within the MCS are often integrated in inadequate way, leading to over-simplified models [10]. Especially, the correct integration of a computational fluid dynamics (CFD) code and photosynthetic reaction kinetics is essential for meaningful solution [11]. Finally, current models, describing general MCS or specific PBR performance analysis, have usually no connection to the bio-production itself, e.g., [4, 12], and references within.

In this study, we present a unified modeling framework of fluid dynamics and photosynthetic reaction kinetics for microalgae culture system of any geometry. The framework is validated by a numerical simulation of microalgae growth in a 2-D square cavity with one moving wall, a hypothetical system triggering a vortex flow. This system is able to reproduce the hydrodynamically induced high frequency light-dark cycles regime causing the so-called flashing light enhancement [1315].

Methods

We present a simulation tool for MCS (photobioreactor, open or raceway ponds) to connect the key aspect of real system, i.e., microalgae cell growth under changing hydrodynamical and optical conditions. This framework consists of (i) a state system – mass balance equations in form of advection-diffusion-reaction partial differential equations – PDEs (1-3), (ii) the fluid flow equations, i.e., Navier-Stokes equations (4), and (iii) the irradiance distribution inside MCS according to Beer-Lambert’s law (5), see [8, 16].

All mentioned parts of the introduced MCS model are interconnected, i.e., the mass balance equations for the state variables (species characteristics, nutrients, gases, etc) are solved simultaneously with the fluid dynamics (momentum balances, continuity equations). Nevertheless, we assume that the stationary flow field inside the MCS is not affected by mass transfer and reactions [17]. This assumption allows the detachment of biological and environmental states and thus biological and environmental parts (models) can be analyzed with different numerical methods and with different spatial-temporal discretizations, dramatically reducing the computational demands. This detachment could be complete (for the stationary flow field in a continuous system) or stepwise, e.g., reflecting some sequence of quasi-steady states in a production system operated in batch mode.

State model

The biological part of modeling framework is based on the mass balance equations for the state variables describing the transport and reaction among the species or compounds [18]:
$$ \begin{aligned} \frac{\partial c_{i}}{\partial t} + \nabla \cdot ({v}c_{i}) - \nabla \cdot (D_{e}\nabla c_{i}) = R(c_{i}) + S(c_{i})~,\\ ~ t \in [\!t_{0},T],~ i=1,\ldots,m, \end{aligned} $$
(1)

where ci=ci(x,t), is a conservative quantity (concentration or cell density), v is a velocity of flow field ruled by the fluid-dynamic model, cf. (4), and xΩR3 stands for a position vector in a coordinate system (e.g., Cartesian or cylindrical). The dispersion coefficient De is a second-order tensor, which corresponds to the diffusion coefficient in microstructure description. De is an empirical parameter describing mixing in the system and it is influenced by the molecular diffusion and velocity profile, i.e., it is not strictly a material constant.

The reaction kinetics and temporal changes of the reacting species are described by the reaction term R(ci) and source terms S(ci), respectively. The source term, e.g., the external load of nutrients into MCS, is usually modeled as a corresponding boundary condition. However, in order to (i) simplify the analytic study of the optimal solution existence, see [19], and (ii) respect that the location of discharge of some material could be inside the domain Ω, we employed the particular form of (1). The initial condition and boundary condition (impermeability of the domain boundary, e.g., PBR walls) are following:
$$ ~ c_{i_{0}} = c_{i}(x,t_{0}), ~x \in \Omega \subset R^{3}, ~ i=1,\ldots,m, $$
(2)
$$ ~\nabla c_{i}(x, t) = 0,~ x \in \partial \Omega, ~ t \in [t_{0},T],~ i=1,\ldots,m. $$
(3)

Fluid-dynamic model

Microalgal cells are considered to be solid particles and fixed CO2 and released O2 are gases, therefore the analyzed system can be described as a multiphase flow and transport. However, it is possible to consider our multiphase system as a suspension by neglecting the gaseous phase. Then, the microalgal cells represent the dispersed phase of a suspension. Also, considering the cell density about 10 kg m −3 for the dry weight of biomass, i.e., 1% of mass content, allows us to classify our system as a single-phase system based on Newtonian viscosity relationship within the employed computational software ANSYS Fluent [20]. Moreover, we assume that the mean diameter of a spherical microalgae cell is approximately ten micrometers [5] which supports this simplification.

Mass density of the suspension is determined as ρ=ρw (1−k)+ρs k, where ρw is the mass density of the medium and ρs represents the cell mass density, and k is the volume fraction. However, one can assume ρ=ρw because of a homogeneous distribution of algal cells (no aggregation) which are actually floating in the medium. Furthermore, the inter-particle distances in our case of dilute suspension are sufficiently large to calculate flow field over each particle or cell [18], i.e., particles do not interfere with the flow field. This idea has been proposed by J. Pruvost et al. [21] who showed that the smallest eddy size, based on Kolmogorov scale, is approximately ten times bigger than the cell size, i.e., neglecting the cell-vortex interaction.

Therefore, we model the incompressible liquid phase (suspension of water, nutrients and microalgae), thus the classical system of Navier-Stokes equations and the continuity equation is employed as:
$$ \frac{\partial v}{\partial t} + (v \cdot \nabla) {v} = f - \frac{1}{\rho}\nabla p + \nu \nabla^{2} v, ~~~ \nabla \cdot v = 0, $$
(4)

in [ t0,T],×Ω, with suitable boundary conditions on [ t0,T],×Ω and initial conditions in Ω. v,p,f,ρ and ν denote the fluid velocity, the pressure, the body forces, fluid density and kinematic viscosity, respectively.

Reaction model and its time-scale analysis

The photosynthetic reactions depend on many variables such as irradiance, accessibility of both inorganic (Fe, P, etc) and organic nutrients (e.g., carbohydrates), concentration of gases, osmolality, and many others. In this work, we focus only on the key variable determining the reaction kinetics, the fluctuations of irradiance; the other variables are assumed to be externally controlled to not limit the cellular growth. Generally, both light attenuation and scattering should be considered to describe the irradiance distribution inside MCS [16]. Moreover, for complicated MCS geometries and multiple light sources distribution, the irradiance field has to be determined experimentally. Here, for the sake of simplicity and manageable computational costs, we suppose the Beer-Lambert law (5) is valid, i.e., the incident irradiance I0 is exponentially decreasing with the increasing light path s:
$$ I(s) = I_{0} ~e^{-\Lambda s}, $$
(5)

where Λ is an attenuation coefficient describing the attenuation of irradiance by the suspension of microbial cells in the liquid medium (unit: m −1).

Based on the decades of experimental research of photosynthesis, the reliable model of reaction kinetics in MCS should implement no less than three time scales across the photosynthetic reactions, see [22]. The following systems are integrated in the reaction term R(ci) in (1): (i) activation of the light harvesting complex (photosynthetic unit – PSU) in light reactions, (ii) biomass output, and (iii) photoinhibition, i.e., damage of pigment-protein complexes by excessive irradiance. Photoacclimation (light adaptation) is another important process which, however, can be observed in the tens of hours-days time scale. Therefore, it is possible to model it as a (quasi) steady state mechanism, e.g., implicitly via a reaction model parameter when another differential equation describes dynamics of this parameter [23].

Comparison of the characteristic times of microalgae growth ( hours) with the mixing due to dispersion – turbulent diffusion ( seconds) indicates a huge gap between both processes. There is a possibility to adopt a very simple growth modeling approach (e.g., Monod or Haldane kinetics) leading to two alternatives. The first one is governed by the time scale of chosen steady state kinetics and has to neglect details concerning mixing phenomena. The second alternative, based on the time step in order of seconds, can observe transitions because of the hydrodynamic mixing but has to neglect the changes in photosynthetic growth (in form of a growth kinetic model). However, both alternatives entirely miss the connection across the transport and reaction processes.

As a suitable candidate for the reaction model, we adopted a mechanistic three-state model of photosynthetic factory (PSF) [2429]. PSF model considers 3 states in which microalgae cells may exist: activated – A, inhibited – B and rested – R. Although the fluid-dynamical properties of cells in each of three states are exactly the same, these states are impacted differently by dynamic changes of the environment, i.e., spatio-temporal changes of states concentrations determine the biomass production, see (6).

Let be the concentrations of respective components cA, cB, and cR (identical units as for the microalgae cell density cx in whole PBR – generally 106 cell ml −1 as in [28]). Then the following relation holds (for t[t0,t], and xΩ): cA(x,t)+cB(x,t)+cR(x,t)=cx(x,t). The dimensionless scalars yA=cA/cx, yB=cB/cx and yR=cR/cx (molar fractions) are respective states of the PSF model, see the next section.

We note that the spatio-temporal averaged rate of photosynthetic production (specific growth rate \(\mu := \dot {c_{x}}/c_{x}\)) is related (constant κ) to the activated state fraction [28, 30]:
$$ \dot c_{x} =\frac{\kappa \gamma}{\text{meas}(\Omega) ~T} \int_{0}^{T} \int_{\Omega}^{}\left [{y_{A}(I(x),t) - Me}\right ] c_{x} ~\mathrm{d}x~ \mathrm{d}t~, $$
(6)

where term Me represents a cellular maintenance; Me is considered constant although it may vary if hydrodynamical shear stress is considered [29].

Decimal value 10−4 [s −1] of term κγ in (6) allows the shift from the time scale of light fluctuation (caused by fluid dynamics or the light source via the irradiance distribution I(x) in PSF model) to the time scale of biomass growth (macro-scale) with no impact on the accuracy. Then, the integral in (6) is evaluated exactly and factor κγ in (6) provides the value of a specific growth rate. Thus, the employed PSF model satisfies the requirement for the time scales of the model candidate. Three of its five parameters can be seen as time constants for three different time scales: i) the light and dark reactions ( seconds), ii) the photoinhibition ( minutes) and iii) the microalgae growth ( hours). The remaining two parameters represent the optimal irradiance and the conversion from dimensionless growth to the real values, see [27, 31] for more details.

Model of photosynthetic factory – PSF model

Here we characterize the multi time-scale three-state model of photosynthetic factory (PSF model) proposed by Eilers and Peeters [24, 25] and further developed by Wu and Merchuk [28, 29] and Celikovsky, Papacek and Rehak [26, 27, 3032]. PSF model (Fig. 1), is used for derivation of the reaction term R(ci), see the transport Eq. (1).
Fig. 1
Fig. 1

Scheme of Photosynthetic factory model (PSF). Four parameters (α, β, γ, δ) describe transitions among respective states (resting R, activated A, and B inhibited)

PSF model (7) incorporates the dynamics of three fundamental states propagating in different orders of time scales: (i) cell growth (including the shear stress effect), (ii) photoinhibition, and (iii) photosynthetic light and dark reactions. The state vector of the PSF model is three dimensional, y=(yR,yA,yB), where yR indicates the probability that PSF is in the resting state R, yA the probability that PSF is in the activated state A, and yB the probability that PSF is in the inhibited state B.
$$ \begin{array}{c} \left[ \begin{array}{c} \dot{y_{R}} \\ \dot {y_{A}} \\ \dot {y_{B}} \end{array} \right] \,=\, \left[ \begin{array}{rrr} 0 & \gamma & \delta \\ 0 & - \gamma & 0 \\ 0 & 0 & - \delta \end{array} \right] \left[ \begin{array}{c} y_{R} \\ {y_{A}} \\ {y_{B}} \end{array} \right]\\ \qquad \qquad\qquad \,\,\,\,+ u(t)\left[ \begin{array}{rrr} - \alpha & 0 & 0 \\ \alpha & - \beta & 0 \\ 0 & \beta & 0 \end{array} \right] \left[ \begin{array}{c} y_{R} \\ {y_{A}} \\ {y_{B}} \end{array} \right] \end{array} $$
(7)

The values of PSF model parameters, i.e., α, β, γ, δ, see (7), and κ, see (6), are taken from original study [28], where Wu and Merchuk identified these values of PSF model parameters for the microalga Porphyridium sp.: α = 1.935 × 10−3 μE−1 m2, β = 5.785 × 10−7 μE−1 m2, γ = 1.460 × 10−1 s −1, δ = 4.796 × 10 −4 s −1, κ = 3.647 10 −3 and Me = 0.059 h −1. For details regarding the experimental setup, for parameter estimation and the identifiability study, see [31, 32].

In order to make the PSF model parameter estimation more robust, the following re-parametrization was introduced and the singular perturbation method was used in [31]: \(q_{1}: = \sqrt {\frac { \gamma \delta }{ \alpha \beta }}, ~ q_{2}: = \sqrt {\frac {\alpha \beta \gamma }{\delta {(\alpha +\beta)}^{2}}},~ q_{3} := \kappa \gamma \sqrt {\frac { \alpha \delta }{\beta \gamma }},~ q_{4}:=\alpha q_{1},~ q_{5} := \beta /\alpha.\) Consequently, the PSF model acquires the following form:
$$ \dot{y} = \left[ {\mathcal{A}}+u(t){\mathcal{B}} \right] y, $$
(8)
$$ {\mathcal{A}}=q_{4} \left[ \begin{array}{ccc} 0 & {q_{2}(1+{q_{5}})} & \frac{q_{5}}{q_{2}(1+{q_{5}})} \\ 0 & - {q_{2}(1+{q_{5}})} & 0 \\ 0 & 0 & - \frac{q_{5}}{q_{2}(1+{q_{5}})} \end{array} \right], $$
(9)
$$ ~~ {\mathcal B}=q_{4} \left[ \begin{array}{ccc} - 1 & 0 & 0 \\ 1 & - q_{5} & 0 \\ 0 & q_{5} & 0 \end{array} \right]. $$
(10)
For the given input variable, i.e., the irradiance u(t), the ODE system (7) can be solved either by numerical methods or by asymptotic methods. However, the irradiance level of 250 μE m −2 s −1 (the value of the parameter q1 which maximizes the steady-state growth [31]), provides two negative eigenvalues of the (7) system matrix: λ1 = -0.63, λ2 = -0.59 10 −3. The ODE system (7) is stiff due to value 103 of the ratio \(\frac {\lambda _{1}}{\lambda _{2}}\). This fact points to the existence of two processes: (i) photosynthetic reactions, and (ii) photoinhibition, widely separated from the point of view of their characteristic times. Further, without significant loss of accuracy, we use the so-called fast reduction [27], i.e., the behavior of the system (7) is characterized by the only one following ODE
$$ \begin{array}{c} \frac{\mathrm{d} y_{A}}{\mathrm{d}t} = - q_{4} \left(u_{av} + q_{2}\right)\left[\frac{u(z) + q_{2}}{u_{av} + q_{2}}y_{A} - \frac{u(z)}{u_{av}} {y_{A}}_{ss}(u_{av})~\right], \end{array} $$
(11)
and the “slow” variable yB can be regarded as a constant depending on the averaged value u=uav. According to our works [27, 33] it holds
$$ \begin{array}{c} {{y_{A}}_{ss}(u)} = \frac{u/q_{2}}{ u^{2}+u/q_{2}+1}, \\ ~{{y_{B}}_{ss}(u)} = \frac{u^{2}}{ u^{2}+u/q_{2}+1} ~. \end{array} $$
(12)

Both the non-reduced (7) and the reduced order PSF model Eq. (11) have been incorporated as a User-Defined Function (UDF) in ANSYS Fluent, see the following section.

Numerical aspects of PDEs implementation and solution

The key issue is in the incorporation of a CFD code and photosynthetic reaction kinetics in one modeling framework. A proper analysis of characteristic times of microalgae growth and of mixing due to the convective motion is especially important. Two appealing approaches, Lattice Boltzmann method [34] and hybrid multicompartment/CFD approach [35], were previously used for bioreactor modeling [36, 37]. Both methods seem to be computationally efficient, however, the discussion of their advantages and disadvantages is far of the scope of this paper. Our approach is based on the implementation of a User-Defined Function (UDF) within the commercial CFD code ANSYS Fluent, which provides the possibility to define an arbitrary reaction term, see below.

As denoted in the previous section, ODE system (7) can be simplified to only one differential Eq. (11). The right-hand part of this equation indicates the rate of change of activated state yA which can be used in the definition of UDF function in ANSYS Fluent. Macro DEFINE_VR_RATE provides the possibility to define the arbitrary reaction term in (1), including its dependency on the spatial coordinate as represented by (5). Macro C_CENTROID can retrieve the corresponding coordinates of the current mesh cell. Three states of the microalgae culture were represented as individual species with the same molar weights in ANSYS Fluent. Therefore, mass and molar fractions should be identical. Example of the user-defined reaction rate definition employed in our simulations for the reduced system (11) follows:

We performed a grid independence study for three mesh sizes of the given geometry (2-D square cavity with moving wall), see Fig. 2. Grid Convergence Index (GCI) describing an average velocity of the flow field for the finest mesh was evaluated as 1.14% according to Celik et al. [38]. All simulations of the 2-D geometry were performed on a grid with approx. 10 thousand mesh elements. Refined mesh elements near walls were created so that dimensionless thickness Y+ was around 1 to capture properly the viscous sublayer in the turbulent flow regime. Time step 0.025 s was used and the accuracy of the solution at the end of time range 0 – 5 s was evaluated as 0.19%.
Fig. 2
Fig. 2

Contours of steady-state velocity magnitude of the microalgae suspension within the idealized 2-D square cavity. The vortex flow is imposed by the moving upper wall (wall coordinate is L and its velocity is vL in x-direction, i.e., from left to right), corresponding Reynolds number is Re=1000 in this case (laminar flow regime)

Every simulation consisted of two steps. In the first step, the Navier-Stokes Eq. (4) were solved iteratively to get a converged solution of the steady-state flow field for the corresponding mixing rate (Reynolds number). In the second step, the flow field was assumed to be fully developed and only Eq. (1) describing the transport and reaction of individual species (states of the algae in our case) were solved. This approach is more efficient but in cases when the flow field exhibits a transient behaviour, the Navier-Stokes equations must be solved simultaneously with the species transport equations. For example, the Coutte-Taylor device shows an inherently transient behaviour at some specific rotations rates [39]. Periodic oscillations can be observed even in the square cavity [40, 41] which we neglected in this study.

Results and discussion

Our early works [27, 30] theoretically confirmed the phenomenon of flashing light enhancement, known from experiments [1315]. However, these models were applicable only for specific PBR geometries, being a Couette-Taylor bioreactor in [30] and flat panel PBR in [27], and suffered from over-simplification (Navier-Stokes equations were not solved there within a fluid dynamic model). Here, for the first time, we introduce a general approach applicable to any geometry of MCS and thus ready for industrial applications. For the validation of our integrated model, we have chosen the special geometry consisting of a 2-D square cavity with one moving wall (Fig. 2), representing an idealized part of a Couette-Taylor device, see Fig. 3. We note that microalgae cells are forced to travel between dark and light sides of the Couette-Taylor device (see Taylor vortex in [39]) which is an analogy to the flashing light regime and consequently the flashing light enhancement [15]. Indeed, very high cell densities (about 150 g/l) were reported in this type of PBR by E. A. Davis in [42]. The author claims that the growth enhancement is due to an ordered mixing (E. A. Davis misuses the term “turbulence”) produced by the Couette-Taylor device. Nevertheless, the mixing enhanced light and mass transfer is only one part of the whole picture. Based on the fact that the standing issue of microalgae mass production in PBR is the cell fragility [43], we see also the other enhancing effect - low shear stress transmitted to microbial cells in the laminar Taylor vortex flow regime [44].
Fig. 3
Fig. 3

Result of CFD simulation of one particle trajectory by CFD code ANSYS Fluent [20]. The particle trajectory in the Couette-Taylor device cross-section is shown on the left side. The right-side picture describes the time course of the microalgae particle radial position R. The “irradiance history” of an individual cell can be calculated by concatenation of trajectory and the irradiance field I=f(R,t) within the device

Photosynthesis in fluctuating light regime

The growth of photosynthetic microorganisms could be described by the steady-state photosynthetic reaction kinetics called P–I (photosynthesis-irradiance) curve representing the photosynthetic response in small cultivation systems with a homogeneous light distribution [17]. However, MCS models based on P–I curves only do not take into account the effect of hydrodynamic mixing (or mixing via static mixers) which has a different characteristic time scale compared to the time scale of the steady-state reaction kinetics, see “Reaction model and its time-scale analysis” section and [27, 31] for more details. In order to overcome this problem, Terry in [15] employed “artificial” corrections to the steady-state kinetic model reflecting the enhancement in the photosynthetic efficiency during fluctuating light conditions present in the dynamic model.

In our work, we aim to build a general predictive model, independent of specific experimental data. Such model must be built on universal principle that time-scales of all relevant processes are considered; see details in section “Reaction model and its time-scale analysis”.

To determine if our modeling framework is able to describe the flashing light enhancement phenomenon, we have chosen a 2-D square cavity with the moving top wall as our hypothetical MCS, see Fig. 2. This MCS can be viewed as the representation of a Couette-Taylor device, i.e., one cell from its axial cross-section when the so-called Taylor vortex flow regime takes place. This occurs when the Taylor number exceeds the first critical value [39], see Fig. 4. The cultivation system is illuminated from the bottom side only by the incident irradiance I0 and the irradiance level is exponentially decreasing with the increasing z-coordinate until the top (z=L) according to the Beer-Lambert law (5). It is convenient to work further with the normalized irradiance u=I/Iopt, where Iopt reaches the value of 250 μE m −2 s −1 [28]. Our analysis showed that the value of the averaged irradiance within MCS uav=1, see Table 1, provides an optimal condition for the microalgae culture growth. This optimization problem depends mainly on mixing intensity and it is discussed elsewhere [33]. Only for the lumped parameter system, the value u=1 is the optimal one. The relations for the irradiance level in depth z and for the average (absorbed) irradiance in the whole MCS are then:
$$ u(z) = u_{0} ~e^{-\Lambda z}~,~ u_{av} = u_{0}\frac{1 - e^{-\Lambda L}}{\Lambda L}~, $$
(13)
Fig. 4
Fig. 4

Fluid flow velocity profile in the axial section of the laboratory Couette-Taylor device calculated by CFD code ANSYS Fluent 15.0 (laminar model, inner cylinder angular frequency ω = 2.4 rad s −1, Re=2000). One rectangular cell used for the 2D case study is marked by the bold rectangle

Table 1

Basic parameters specifying the growth rate in MCS system used in this study and CFD simulations: uav represents the value of average irradiance in the culture, L stands for MCS characteristic length, velocity of the upper MCS wall is given by Reynolds number Re, ΛL stands for optical thickness, q2, q4, yR(t0), yA(t0) and yBss(uav=1) represent 2 parameters, 2 initial conditions for rested and activated state of the PSF model and (quasi) steady state level of the inhibited state, respectively (see section “Model of photosynthetic factory – PSF model”)

u av

L [m]

Re =vL L/ν

Λ L

q 2

q4 [s−1]

yR(t0)

yA(t0)

yBss(uav)

1

0.02

0 – 100000

8 ln(2)

0.3

0.5

1

0

0.1875

where u0=I0/Iopt is the incident irradiance on the bottom (z=0) and L is the MCS thickness in direction of light gradient. It is further convenient to define a dimensionless “optical thickness” of the system Ot>0 as follows: Ot:=ΛL. It holds \(\frac {r_{1/2}}{L}=\frac {\ln (2)}{Ot}\), where r1/2 is the length interval (unit: m) in which the light intensity diminishes to one half. Then, ratio \(\frac {u(z)}{u_{av}}\) for optically thin culture system (Ot≈5.5 in our case) has
$$\frac{u(z)}{u_{av}}\approx {\Lambda L}~{e^{-\Lambda L \frac{z}{L}}}. $$
(14)

We suggest that the “irradiance history” of an individual cell within our system in Fig. 2 corresponds to the hydrodynamically induced light/dark cycles shown in Fig. 3 (produced by the Couette-Taylor device). The irradiance history can be calculated by linking the trajectory of an individual microalgae cell and the irradiance field within device.

The final step is the implementation of the transport-reaction model based on the PDEs (1-4). The key role plays the reaction term R(ci) in (1). The R(ci) term is “local” in sense that it depends on the actual transport and reaction rates, and it contains all relevant time scales, thus, it perfectly fits to our requirements. The technical simplicity of PDE based approach, which takes advantage of the sophisticated CFD codes, e.g., ANSYS Fluent [20], is obvious when compared to the alternative Lagrangian approach. While the Eulerian approach reaches its governing equations from the mass balance in an infinitesimal control volume, the Lagrangian approach is based on the above mentioned irradiance history of an individual cell. Although conceptually simple, the application of Lagrangian approach to MCS modeling presents serious problems. First, either an experimental technique, e.g., computer-automated radioactive particle tracking (CARPT) [45] or a CFD code is needed for the cell trajectory description, see Fig. 3. Second, using an ensemble of cell trajectories across the irradiance field inside MCS, the probabilistic description of “irradiance history” of individual cells is reached. Finally, the problem reduces to a reaction kinetics model represented by an ordinary differential equation with a stochastic input.

The employed Eulerian approach is more convenient, mainly because it benefits from CFD codes and the implementation of an appropriate dynamic model of microbial growth (with lumped parameters) into a distributed parameter system is merely technical.

Case study: Biological performance of MCS in flashing-light regime

Having the modeling framework assembled, we aimed to reproduce the hydrodynamic regime where high frequency light-dark cycles are induced and trigger the flashing light enhancement [14, 15]. The values of parameters employed in the simulation of MCS performance are summarized in Table 1.

The PSF model consists of the ODE system (7), nevertheless it can be reduced by using the singular perturbation method [33]. Then, the behavior of the system (7) is characterized by only one ODE (11). Furthermore, for the local description of PSF model, we need the expression for u(z) in the case of average irradiance uav=1; based on (14): \(u(z)\approx \Lambda L~{e^{-\Lambda L \frac {z}{L}}}\).

For the initial condition and other parameters in Table 1, the case study was resolved using CFD code ANSYS Fluent: PDE (1-4) embedded and the reaction kinetics implemented as a special UDF (User-Defined n) according to the description in section “Numerical aspects of PDEs implementation and solution”. The resulting steady-state solution for molar fraction yA is shown in Fig. 5. The spatial distribution of molar fraction yB steady-state solution is shown in Fig. 6. The spatio-temporal average for molar fraction yA is shown in Fig. 7 and, finally, the temporal evolution of yB steady-state is depicted in Fig. 8 (right). Notice nice patterns of yB for low mixing conditions in Fig. 6 (left and center).
Fig. 5
Fig. 5

The steady-state spatial distribution of the activated state molar fraction yA for three different Reynolds numbers, Re = 0, 1000, 100000, from left to right. The irradiation, governed by the Lambert-Beer law (13), is maximal on the bottom side. Contour plots of the state yA were calculated by the CFD code ANSYS Fluent 15.0 using the non-reduced (full) PSF model (7). The dynamics of yA state (and yB as well) is depicted in Fig. 8

Fig. 6
Fig. 6

Contour plots of the inhibited state yB calculated by the CFD code ANSYS Fluent 15.0 using the non-reduced model (7) for different mixing rates. The steady-state spatial distribution of yB shows the clear pattern for low mixing conditions corresponding to no mixing (Left: Re = 0) and low mixing (Center: Re = 1000), while for an almost ideally mixed case (Right: Re = 100000) the yB concentration is evenly distributed within the reaction domain. The numerical values are clearly consistent with the results shown on the right part of Fig. 8 (non-reduced model)

Fig. 7
Fig. 7

Simulated and experimental results of the flashing light enhancement phenomenon. Left picture shows experimental data reprinted from [14]. Two curves of Oxygen evolution rate vs. L-D cycles period are tagged with the ratio of light-to-dark interval L:D. Right: the normalized performance index J vs. Da, for larger mixing rate (i.e., lower Da - bigger Re), we get better performance J approaching its maximal theoretical value corresponding to the growth in averaged continuous light, cf. (12). Red circles represents simulation results of the whole non-reduced system (7) or (8) gained for time range 0 - 20000 s, with time step 0.1 s

Fig. 8
Fig. 8

Time dependency of the activated state yA for three different mixing intensities, Re = 0, 1000, 100000, with and without fast reduction, respectively. On the left, the time dependency of activated state yA of the reduced system (11) is shown for time range 0 - 40 s. On the right, the time dependency of activated state yA and inhibited state yB of the non-reduced system (7) is shown for time range 0 - 20000 s. The dynamics of activated state yA has two phases, fast phase (the transition from the initial value yA(t0)=0 to the maximal value occurring in seconds) and slow phase (imposed by the gradually growing inhibited state yB). In the case of applied fast reduction (left graph), the slow dynamics does not occur due to constant

It is clear that the presented simulations are in the qualitative match with the previously published experimental data [14], see Fig. 7 on the left. Also, the time scale of yA state dynamics is shown in Fig. 8, describing the time course of the spatial average of yA for three different velocities vL, corresponding to Reynolds number of 0, 1000, and 100000, respectively. The mixing intensity determines the frequency of hydrodynamically induced light/dark cycles. Similarly to our previous simplified studies [27, 30], we can see the hydrodynamically induced growth enhancement due to the flashing light. The major difference and improvement to previous studies is that presented modeling framework can be used to any MCS setup and not only to special cases.

Analysis of our numerical model indicates that the process is characterized by the Damköhler number Da:=ttr/tr which describes the ratio between the characteristic time of transport (\(t_{tr}=\frac {L}{v_{L}}\)) and reaction time (tr=1/[q4(uav+q2)]). In chemical engineering, the Damköhler numbers Da are frequently used to describe the ratio between the chemical reaction time scale (reaction rate) and the flow time scale (convective mass transport rate) occurring in a system. For our specific case, when the characteristic time tr is fixed, it holds:
$${\text{Da}}:=1/t_{r} \frac{L^{2}}{\nu}\frac{1}{\text{Re}}=\frac{260}{\text{Re}}. $$

It is known that extremely high or low values of the characteristic numbers leads to numerical difficulties during solving PDEs, e.g., (1). Employing Monod or Haldane kinetics (steady-state) might lead to such difficulties and that is why a multi time-scale model would be more appropriate.

In order to quantify the impact of mixing velocity on the cellular growth, we defined an objective function JMCS as a volumetric productivity (usually in grams per liter per day):
$$ J_{MCS} = \frac{1}{\text{meas}(\Omega) ~T} \int_{0}^{T} \int_{\Omega}^{}\left[{y_{A}(u,z,t)}\right] c_{x} ~\mathrm{d}z \mathrm{d}t~. $$
(15)
For a usual case of steady state operation, it is only the activated state fraction yA(u,z,t) in JMCS which has to be integrated over the MCS domain. Thus, in Fig. 7 we show the dependence of the normalized performance index
$$\frac{J}{c_{x}}=\frac{1}{L^{2}}\int_{\Omega}^{}{y_{A}(u,z,t)} ~\mathrm{d}z ~. $$

We note that the empirical data have an illustrative and testing purpose only and for quantitative simulations one has to employ parameters for specific MCS and particular microalgae. However, the result of biological performance for our MCS compared to experimental results shown in Fig. 7 (left) clearly shows that the positive impact of mixing on the growth induced by the convective motion, i.e., the simulated flashing light enhancement is comparable to experiments. Analyzing simulations in Fig. 7 (right) implies that lower Da (bigger v, better mixing) leads to better performance J (approaching its maximal theoretical value 0.625). Nevertheless, the harmful impact of the hydrodynamically induced shear stress originally considered in Me term [29], see (6), is not taken into account and thus in reality certain mixing intensity would start damaging the cellular integrity.

Conclusions

In order to reach commercially sustainable production for algae biofuels, the empirical data indicates the need for accurate and reliable modeling approach with both predictive and optimizing capabilities. In this work, we presented the general unified modeling framework for microalgae culture systems (MCS) which can be applied to an arbitrary geometry of the production system and any microbial strain in general. This allows evaluating the performance of various MCS architectures, operating conditions, and microalgae strains as well. All parts of the multidisciplinary framework, i.e., the growth-rate model, the fluid-dynamic model, and the irradiance-distribution model, are integrated within CFD package ANSYS Fluent. In the illustrative case study, we have shown that a three-state model of photosynthetic factory (PSF model) well behaves under hydrodynamically induced high frequency light-dark cycles regime and copes with the requirement imposed on the reaction model, i.e., it correctly describes both the quasi steady-state (flashing light enhancement) and dynamic phenomena (time course of the PSF model activated state).

We can conclude that our modeling framework provides an adequate and physically accurate description of microalgae growth in a MCS and thus it is suited for the optimal control problem formulation as well. To the extent of our knowledge, this is the first time when such advanced solution was successfully used in biotechnology of microalgae.

Our future plans are to further focus on the optimization problem and its formulation. Finally, we plan to analyze the shear stress impact on the microalgae growth.

Abbreviations

CARPT: 

Computer-automated radioactive particle tracking

CFD: 

Computational fluid dynamics

GCI: 

Grid convergence index

MCS: 

Microalgae culture systems

PBR: 

Photobioreactors

PDEs: 

Partial differential equations

PI: 

Photosynthesis-irradiance

PSF: 

Photosynthetic factory

UDF: 

User-defined function

Declarations

Funding

This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic - projects ”CENAKVA” (No. CZ.1.05/2.1.00/01.0024), ”CENAKVA II” (No. LO1205 under the NPU I program) and ’The CENAKVA Centre Development’ (No. CZ.1.05/2.1.00/19.0380). Publication costs were funded by University of South Bohemia.

About this supplement

This article has been published as part of BMC Systems Biology Volume 12 Supplement 5, 2018: Selected articles from the 5th International Work-Conference on Bioinformatics and Biomedical Engineering: systems biology. The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume-12-supplement-5.

Authors’ contributions

SP designed the model and drafted the manuscript, KP performed the simulations and collaborated in drafting the manuscript. JJ validated the results and supervised the manuscript preparation. All authors interpreted part of the data and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
Institute of Complex Systems, South Bohemian Research Center of Aquaculture and Biodiversity of Hydrocenoses, Faculty of Fisheries and Protection of Waters, University of South Bohemia in České Budějovice, Zámek 136, 373 33 Nové Hrady, Czech Republic
(2)
Czech Technical University in Prague, Faculty of Mechanical Engineering, Technická 4, Prague, 160 00, Czech Republic

References

  1. Su Y, Song K, Zhang P, Su Y, Cheng J, Chen X. Progress of microalgae biofuel’s commercialization. Renew Sust Energ Rev. 2017; 74:402–11. https://doi.org/10.1016/j.rser.2016.12.078.View ArticleGoogle Scholar
  2. Burlew J. Algal Culture From Laboratory to Pilot Plant. Washington: Carnegie Institution of Washington; 1953.Google Scholar
  3. Grobbelaar JU. Microalgal biomass production: challenges and realities. Photosynth Res. 2010; 106(1):135–44. https://doi.org/10.1007/s11120-010-9573-5.View ArticleGoogle Scholar
  4. Ooms MD, Dinh CT, Sargent EH, Sinton D. Photon management for augmented photosynthesis. Nat Commun. 2016; 7:12699. https://doi.org/10.1038/ncomms12699.View ArticleGoogle Scholar
  5. Richmond A. Handbook of Microalgal Culture: Biotechnology and Applied Phycology. Oxford: Blackwell Publishing; 2008.Google Scholar
  6. Chisti Y. Biodiesel from microalgae, biotechnology advancesbiodiesel from microalgae. Biotechnol Adv. 2007; 25:294–3062.View ArticleGoogle Scholar
  7. Bernardi A, Perin G, Sforza E, Galvanin F, Morosinotto T, Bezzo F. An identifiable state model to describe light intensity influence on microalgae growth. Ind Eng Chem Res. 2014; 53(16):6738–49. https://doi.org/10.1021/ie500523z. PMID: 25678739. https://doi.org/10.1021/ie500523z.View ArticleGoogle Scholar
  8. Park S, Li Y. Integration of biological kinetics and computational fluid dynamics to model the growth of nannochloropsis salina in an open channel raceway. Biotech Bioeng. 2015; 112(5):923–33. https://doi.org/10.1002/bit.25509.View ArticleGoogle Scholar
  9. Sato T, Yamada D, Hirabayashi S. Development of virtual photobioreactor for microalgae culture considering turbulent flow and flashing light effect. Energy Convers Manag. 2010; 51:31196–201.View ArticleGoogle Scholar
  10. Muller-Feuga A, Guédes RL, Pruvost J. Benefits and limitations of modeling for optimization of Porphyridium cruentum cultures in an annular photobioreactor. J Biotechnol. 2003; 103:153–63.View ArticleGoogle Scholar
  11. Bitog JP, Lee I-B, Lee C-G, Kim K-S, Hwang H-S, Hong S-W, Seo I-H, Kwon K-S, Mostafa E. Application of computational fluid dynamics for modeling and designing photobioreactors for microalgae production: A review. Comput Electron Agric. 2011; 76(2):131–47. https://doi.org/10.1016/j.compag.2011.01.015.View ArticleGoogle Scholar
  12. Bernard O, Mairet F, Chachuat B In: Posten C, Feng Chen S, editors. Modelling of Microalgae Culture Systems with Applications to Control and Optimization. Cham: Springer: 2016. p. 59–87. https://doi.org/10.1007/10_2014_287.View ArticleGoogle Scholar
  13. Kok B. Experiments on photosynthesis by chlorella in flashing light In: Burlew J, editor. Algal Culture From Laboratory to Pilot Plant, vol. 600. Washington: Carnegie Institution of Washington: 1953. p. 63–75.Google Scholar
  14. Nedbal L, Tichý V, Xiong F, Grobbelaar JU. Microscopic green algae and cyanobacteria in high-frequency intermittent light. J Appl Phycol. 1996; 8:325–33.View ArticleGoogle Scholar
  15. Terry KL. Photosynthesis in modulated light: Quantitative dependence of photosynthetic enhancement on flashing rate. Biotech Bioeng. 1986; 28:988–95.View ArticleGoogle Scholar
  16. Cornet J-F, Dussap CG, Gros J-B, Binois C, Lasseur C. A simplified monodimensional approach for modeling coupling between radiant light transfer and growth kinetics in photobioreactors. Chem Eng Sci. 1995; 50:1489–500.View ArticleGoogle Scholar
  17. Schugerl K, Bellgardt KH. Bioreaction Engineering, Modeling and Control. Berlin Heidelberg: Springer; 2000.View ArticleGoogle Scholar
  18. Beek WJ, Muttzal KMK, Heuven JWV. Transport Phenomena, 2nd edn. Chichester: Wiley; 1999.Google Scholar
  19. Alvarez-V ázquez LJ, Fernández FJ. Optimal control of bioreactor. Appl Math Comput. 2010; 216:2559–75.Google Scholar
  20. Fluent ANSYS. Theory guide. Ansys Inc. 2015.Google Scholar
  21. Pruvost J, Legrand J, Legentilhomme P, Muller-Feuga A. Simulation of microalgae growth in limiting light conditions: Flow effect. AIChE J. 2002; 48(5):1109–20. https://doi.org/10.1002/aic.690480520.View ArticleGoogle Scholar
  22. Han B-P, Virtanen M, Koponen J, Straškraba M. Effect of photoinhibition on algal photosynthesis: a dynamic model. J Plankton Res. 2000; 22(5):865. https://doi.org/10.1093/plankt/22.5.865.View ArticleGoogle Scholar
  23. Kmeť T, Straškraba M, Mauersberger P. A mechanistic model of the adaptation of phytoplankton photosynthesis. Bull Math Biol. 1993; 55:259–75.View ArticleGoogle Scholar
  24. Eilers PHC, Peeters JCH. A model for the relationship between light intensity and the rate of photosynthesis in phytoplankton. Ecol Model. 1988; 42(3):199–215. https://doi.org/10.1016/0304-3800(88)90057-9.View ArticleGoogle Scholar
  25. Eilers P, Peeters J. Dynamic behaviour of a model for photosynthesis and photoinhibition. Ecol Model. 1993; 69(1-2):113–33.View ArticleGoogle Scholar
  26. Papacek S, Celikovsky S, Stys D, Ruiz-León J. Bilinear system as modelling framework for analysis of microalgal growth. Kybernetika. 2007; 43:1–20.Google Scholar
  27. Papacek S, Matonoha C, Stumbauer V, Stys D. Modelling and simulation of photosynthetic microorganism growth: random walk vs. finite difference method. Math Comput Simul. 2012; 82(10):2022–32. https://doi.org/10.1016/j.matcom.2011.07.006. The Fourth IMACS Conference : Mathematical Modelling and Computational Methods in Applied Sciences and Engineering, Devoted to Owe Axelsson in ocassion of his 75th birthday.View ArticleGoogle Scholar
  28. Wu X, Merchuk JC. A model integrating fluid dynamics in photosynthesis and photoinhibition processes. Chem Eng Sci. 2001; 56:3527–38.View ArticleGoogle Scholar
  29. Wu X, Merchuk JC. Simulation of algae growth in a bench-scale bubble column reactor. Biotech Bioeng. 2002; 80(2):156–68. https://doi.org/10.1002/bit.10350.View ArticleGoogle Scholar
  30. Papacek S, Stumbauer V, Stys D, Petera K, Matonoha C. Growth impact of hydrodynamic dispersion in a couette–taylor bioreactor. Math Comput Model. 2011; 54(7):1791–5. https://doi.org/10.1016/j.mcm.2010.12.022. Mathematical models of addictive behaviour, medicine & engineering.View ArticleGoogle Scholar
  31. Rehak B, Celikovsky S, Papacek S. Model for photosynthesis and photoinhibition: Parameter identification based on the harmonic irradiation O 2 response measurement. TAC IEEE and TCAS IEEE. 2008:101–108.View ArticleGoogle Scholar
  32. Papacek S, Celikovsky S, Rehak B, Stys D. Experimental design for parameter estimation of two time-scale model of photosynthesis and photoinhibition in microalgae. Math Comput Simul. 2010; 80(6):1302–9. https://doi.org/10.1016/j.matcom.2009.06.033. Fifth IMACS Seminar on Monte Carlo Methods Applications of Computer Algebra 2007 (ACA 2007) special session on Nonstandard Applications of Computer Algebra Computational Biomechanics and Biology, a collection of papers presented at the 1st IMACS International Conference on the Computational Biomechanics and Biology ICCBB 2007.View ArticleGoogle Scholar
  33. Celikovsky S, Papacek S, Cervantes-Herrera A, Ruiz-Leon J. Singular perturbation based solution to optimal microalgal growth problem and its infinite time horizon analysis. IEEE Trans Autom Control. 2010; 55(3):767–72. https://doi.org/10.1109/TAC.2010.2040498.View ArticleGoogle Scholar
  34. Succi S. The Lattice Boltzmann Equation. New York: Oxford University Press; 2001.Google Scholar
  35. Bezzo F, Macchietto S, Pantelides CC. Computational issues in hybrid multizonal/computational fluid dynamics models. AIChE J. 2005; 51(4):1169–77. https://doi.org/10.1002/aic.10383.View ArticleGoogle Scholar
  36. Bezzo F, Macchietto S, Pantelides CC. General hybrid multizonal/cfd approach for bioreactor modeling. AIChE J. 2003; 49(8):2133–48. https://doi.org/10.1002/aic.690490821.View ArticleGoogle Scholar
  37. Štumbauer V, Petera K, Štys D. The lattice Boltzmann method in bioreactor design and simulation. Math Comput Model. 2013; 57:1913–8.View ArticleGoogle Scholar
  38. Čelikovský S. Procedure for estimation and reporting of uncertainty due to discretization in cfd applications. J Fluids Eng. 2008; 130(7). https://doi.org/10.1115/1.2960953.
  39. Taylor GI. Stability of a viscous liquid containing between two rotating cylinders. Phil Trans R Soc. 1923; A223:289–343.View ArticleGoogle Scholar
  40. Goodrich JW. An unsteady time-asymptotic flow in the square driven cavity. Technical report tech. mem. 103141. NASA. 1990.Google Scholar
  41. Bruneau C-H, Saad M. The 2d lid-driven cavity problem revisited. Comput Fluids. 2006; 35:326–48.View ArticleGoogle Scholar
  42. Davis EA. Turbulence In: Burlew J, editor. Algal Culture From Laboratory to Pilot Plant. Washington: Carnegie Institution of Washington: 1953. p. 135–138.Google Scholar
  43. Gudin C, Chaumont D. Cell fragility — the key problem of microalgae mass production in closed photobioreactors. Bioresour Technol. 1991; 38(2):145–51. https://doi.org/10.1016/0960-8524(91)90146-B.View ArticleGoogle Scholar
  44. Haut B, Amor HB, Coulon L, Jacquet A, Halloin V. Hydrodynamics and mass transfer in a couette–taylor bioreactor for the culture of animal cells. Chem Eng Sci. 2003; 58(3):777–84. https://doi.org/10.1016/S0009-2509(02)00607-3. 17th International Symposium of Chemical Reaction Engineering (IS CRE 17).View ArticleGoogle Scholar
  45. Luo H-P, Kemoun A, Al-Dahhan MH, Sevilla JMF, Sánchez JLG, Camacho FG, Grima EM. Analysis of photobioreactors for culturing high-value microalgae and cyanobacteria via an advanced diagnostic technique: Carpt. Chem Eng Sci. 2003; 58(12):2519–27. https://doi.org/10.1016/S0009-2509(03)00098-8.View ArticleGoogle Scholar

Copyright

© The Author(s) 2018

Advertisement