ARTICLE Biofilm Growth: A Multi-Scale and Coupled Fluid-Structure Interaction and Mass Transport Approach Mirella Coroneo, Lena Yoshihara, Wolfgang A. Wall Institute for Computational Mechanics, Technische Universita¨t Mu¨nchen, Boltzmannstr. 15 D-85747, Garching, Germany; telephone: þ49-89-289-15236; fax: þ49-89-289-15301; e-mail: [email protected]
Introduction ABSTRACT: In this paper, we propose a novel approach for modelling biofilm growth. It is based on a finite element method and includes both fluid–structure interaction (FSI) as well as scalar transport effects. Due to the different timescales of the involved phenomena, the growth of the biofilm structure is coupled with the FSI and mass transport through a multi-scale approach in time. In each hydrodynamic time step, first the non-linear FSI problem is solved followed by the scalar transport equations, using the information on velocities and deformations obtained in the FSI step. After a steady state solution is reached, information on mass fluxes and stresses are passed to the growth model. At this point, the growth is calculated for a biological time step larger than the hydrodynamic one and based on the mass flux through the interface and on normal and shear stresses on it. This type of approach can significantly contribute to the understanding of biofilm development in fluid flows, since the influence of hydrodynamic conditions and availability of nutrients is well known to have effects on biofilm development. Therefore, for the purpose of understanding biofilm macro-scale dynamics, it is essential to adopt a modeling approach, which takes into account all the relevant aspects, like fluid flow, structure deformation, mass transport and their effect on biofilm growth and erosion. First numerical examples demonstrate the suitability of the proposed model to catch the main features of a growing biofilm structure. Biotechnol. Bioeng. 2014;111: 1385–1395. ß 2014 Wiley Periodicals, Inc. KEYWORDS: biofilm growth; fluid–structure interaction; mass transfer; finite element method
Correspondence to: M. Coroneo Received 25 September 2013; Revision received 3 January 2014; Accepted 6 January 2014 Article first published online 4 February 2014 in Wiley Online Library (http://onlinelibrary.wiley.com/doi/10.1002/bit.25191/abstract). DOI 10.1002/bit.25191
ß 2014 Wiley Periodicals, Inc.
A better understanding of biofilm growth and behavior is essential for the development of successful strategies for biofilm control (Böl et al., 2013), that delay the formation of harmful biofilms or promote the formation of beneficial ones (Teodósio et al., 2012). The formation and development of biofilms, however, is influenced by many factors, including nutrient availability and hydrodynamic conditions (Battin et al., 2003; Garny et al., 2008; Pereira et al., 2002; Stoodley et al., 1999; Wijeyekoon et al., 2004). As a matter of fact, hydrodynamic conditions determine not only the rate at which nutrients are transported to the biofilm surface, but also the magnitude of forces acting on a growing biofilm (Simões et al., 2007), which will also have an influence on the biofilm structure (Hall-Stoodley and Stoodley, 2002). Understanding the mutual interaction of all the involved phenomena, however, has proven to be difficult, because of the many difficulties in carrying out experiments capable of isolating each single mechanism. For this reason, biofilm modeling is arising as a means of producing quantitative tools (Duddu et al., 2009). Biofilms, however, exhibit heterogeneous morphologies (de Beer Schramm, 1999) and the complexity and variability of their matrices impact their physical properties, also significantly increasing the difficulty of biofilm modeling (Pereira et al., 2002). The way they respond to environmental solicitations, however, is defined by their macro-scale characteristics. Indeed biofilms interact with the surrounding fluid as macro-scale materials and can be studied as flexible structures located in liquid flows. For this reason, continuum models, assuming that local properties can be considered smeared on all the spatial extension of the biofilm, can be used and are, moreover, the most suitable for studying big biofilm aggregates. The complex interaction between fluid flow, structure deformation, mass transfer and biofilm growth is fundamental for biofilm development. Most continuum models applied to biofilm modeling, however, neglect fluid-induced deformation of the biofilm structure (Alpkvist and Klapper, 2007a; Dockery and Klapper, 2001; Duddu et al., 2008, 2009;
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
1385
Eberl and Sudarsan, 2008; Graf von der Schulenburg et al., 2009; Radu et al., 2012; Smith et al., 2007), while only few account for it (Alpkvist and Klapper, 2007b; Böl et al., 2009; Taherzadeh et al., 2010; Towler et al., 2007), or for both the substrate transport and the biofilm transient interaction with the surrounding fluid (Cogan, 2008, 2011; Dillon et al., 1996; Klapper, 2004; Taherzadeh et al., 2012). Most of these studies are, moreover, restricted to twodimensions or to specific domain geometries (Cogan, 2011), disregard the advective effect (Alpkvist and Klapper, 2007a; Klapper, 2004), do not take into account the growth (Dillon et al., 1996; Taherzadeh et al., 2012) or calculate it only based on the concentration field (Dockery and Klapper, 2001). For biofilm streamers, for example, calculating the fluid-induced deformations is essential for the purpose of correctly calculating the mass transfer (Taherzadeh et al., 2012) and consequently the growth. For these structures it is also important to apply a multi-scale approach for the purpose of obtaining appropriate and detailed enough models for both the fluid flow and the mass transfer, as well as the biofilm growth. To date a general multi-scale modeling approach taking into account the effect of both local structure deformation and fluid flow on mass transfer and growth is still missing. In this work, we propose a novel multi-scale and coupled fluid-structure interaction (FSI) and mass transport approach for modeling biofilm growth. It is developed in our in-house research code BACI and is aimed at understanding the influence of operating conditions on biofilm growth and erosion and at developing biofilm control strategies.
Fluid–Structure Interaction The interaction between the biofilm structure and the surrounding fluid is modeled through an arbitrary Lagrangian–Eulerian (ALE) approach (Donea et al., 2004). As a result, the FSI problem domain consists of three fields. Apart from the structure VS and the fluid VF domains, sharing a common interface G, a third non-physical mesh field VG,F is present due to the ALE approach (Fig. 1a). Structure Field Biofilms are viscoelastic materials. They behave however elastically, when subjected to external forces over short periods of time, and in a viscous way, when external forces are applied over longer time period (Stoodley et al., 2002). In this work the interaction between fluid and structure is considered only at very small time scales and therefore the biofilm motion can be safely described by an elastic material model. In the structural domain the following non-linear structural elastodynamics equation is solved for the displacement field dS rS
d2 dS ¼ r ðF SÞ þ rS bS dt 2
in
VS ð0; TÞ;
ð1Þ
where rS represents the structural density, bS are the external body forces, VS denotes the undeformed domain, while T is
1386
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
Figure 1.
Definition of domains and their boundaries for the FSI (a) and biofilm growth (b) subproblems. The boundaries of solid and fluid domains consist of Dirichlet, Neumann, and FSI boundary conditions, that is @VI ¼ G ID [ G IN [ G with I 2 (F,S). The boundaries of VG,I are given by @VG;F ¼ G G;F E [ G for the FSI subproblem and by @VG;I ¼ G G;I E [ G g for the biofilm growth subproblem.
the hydrodynamic time period. The equation determines the structural displacements dS by prescribing an equilibrium between forces of inertia, body forces bS, and the internal forces determined from the second Piola–Kirchhoff stress tensor S and the deformation gradient F. The problem is then completed by the imposition of appropriate initial and boundary conditions. For details on nonlinear solid mechanics including stress measures the reader is referred to the monograph of Holzapfel (2000). Fluid Field In the fluid domain the following incompressible timedependent ALE version of Navier–Stokes equations are solved with the purpose of obtaining the pressure pF and the velocity uF fields rF
@uF þ rF ðcF rÞuF 2mr eðuF Þ @t þ rpF ¼ rF bF in VF ð0; TÞ; r uF ¼ 0
in
VF ð0; TÞ:
ð2Þ
ð3Þ
In the momentum Equation (2), e(uF) represents the strain rate tensor of the fluid, m denotes its dynamic viscosity, bF a prescribed body force, and cF the fluid ALE convective velocity, which represents the fluid velocity uF relative to the arbitrarily moving fluid domain cF ¼ uF uG;F :
ð4Þ
at the fluid and structure surfaces
The fluid grid velocity uG,F is defined by uG;F ¼
@w @t
VF ð0; TÞ;
in
ð5Þ
where w is an arbitrary and unique mapping for the deformation of the fluid domain dG,F d
G;F
ðx; tÞ ¼ w dG;F for G ; x; t
ðx; tÞ 2 V ð0; TÞ; ð6Þ F
while dG;F G represents the mesh interface displacement, which will be later related to the structure interface displacement dSG . Also in this case, appropriate initial and boundary conditions have to be set. It is moreover important to highlight that the fluid flow in the deforming fluid domain is considered to be laminar, since biofilm structures usually grow attached to surfaces and are considered to be in the hydrodynamic boundary layer. Hence, Navier–Stokes equations can be solved without any further (turbulence) modelling and with a reasonable computational effort. ALE Field For the purpose of calculating the fluid ALE convective velocity cF appearing in Equation (2), a definition of the mapping w introduced in Equation (6) has to be determined. In the FSI problems studied in the present work, the boundary of the fluid ALE mesh is coupled to the Lagrangian mesh of the structure and to an Eulerian mesh at the in- and outflow portions denoted by GE, that is, dG;F ¼ dS dG;F ¼ 0
on on
G ð0; TÞ;
ð7Þ
G G;F E ð0; TÞ:
ð8Þ
On the other hand, the fluid ALE mesh is allowed to deform arbitrarily within the domain and, in the present study, it was decided to treat the fluid ALE field as a quasielastostatic pseudo-structure (Wall and Ramm, 1998). The ALE equation of motion results then r s G;F ¼ 0
VG;F ð0; TÞ;
ð9Þ
s G;F ¼ lG;F tr eG;F I þ 2mG;F eG;F ; 1 G;F G;F T rd þ rd eG;F ¼ ; 2
ð10Þ
in
with
where tr(·) denotes the trace operator while lG,F and mG,F represent the Lamé constants of the pseudo-structure. FSI Coupling The deformation of the FSI surface G is controlled by kinematic and dynamic constraints. First, the fulfillment of the equilibrium of forces results in equal surfaces tractions h
hSG ¼ hFG
on
G ð0; TÞ:
ð11Þ
Second, at the interface the fluid grid velocity uG;F G and the fluid velocity uFG also have to match uFG ¼ uG;F G
on
G ð0; TÞ:
ð12Þ
Finally, a mass flow across G as well as a relative tangential movement of fluid and structure at G are prohibited, resulting in a fluid velocity uFG equal to the structure deformation rate, that is @dSG ¼ uFG @t
on
G ð0; TÞ:
ð13Þ
Scalar Transport On both the fluid and solid domains, the scalar transport was also solved. For this purpose, the water-substrate solution is considered diluted and in the fluid domain VF a convectiondiffusion equation is solved for the concentration FF @FF þ cF rFF r ðDF rFF Þ ¼ 0 @t in VF ð0; TÞ;
ð14Þ
where cF is the ALE convective velocity, which has already been introduced in Equation (4), and DF is the fluid diffusion coefficient. On the other hand, in the solid domain VS a diffusion-reaction equation in the conservative form is solved dFS þ FS ðr uS Þ r ðDS rFS Þ þ RS ¼ 0 dt in VS ð0; TÞ:
ð15Þ
In the previous equation, uS is the solid velocity, DS represents the solid diffusion coefficient, while RS denotes the reaction term. The microbial growth rate in aqueous environments is a function of the concentration of a limiting nutrient and is usually expressed by the Monod kinetic. If this kinetic is assumed to hold then the substrate consumption rate RS reads RS;M ¼ k
FS ; K þ FS
ð16Þ
where k and K represent the reaction rate constant and the half-saturation concentration of substrate, respectively. Furthermore, it is important to emphasize that, in contrast to the scalar equation for the fluid (14), in Equation (15) the convective velocity is not present, since the structure velocity S uS ¼ @d @t and the grid velocity are identical. Moreover, a conservative formulation is chosen for the mass transport in the solid since in general 5 · uS 6¼ 0.
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1387
The convection-diffusion transport equation in the liquid domain is coupled to the diffusion-reaction transport equation in the structural domain, with the constraints of equal fluxes h and concentrations F at the fluid-structure interface FS ¼ FF hS ¼ hF
on on
G ð0; TÞ; G ð0; TÞ:
ð17Þ ð18Þ
The initial boundary value problem is then completed by appropriate boundary conditions and initial concentration fields. Through the solution of the reported equations, the transport of the substrate from the bulk of the fluid till and through the fluid-structure interface and then its transport and reaction inside the structure is properly modeled. Biofilm Growth The formation and development of different biofilm structures are influenced by flow conditions and nutrients availability (Stoodley et al., 1999). As a result, a reliable model of biofilm growth needs to take into account variables not only connected to mass transfer but also to the flow conditions. The proposed biofilm model calculates the local amount of growth or erosion based on stresses and mass fluxes at the interface resulting from the FSI and mass transfer subproblems. The underlying hypothesis is that growth is the result of the combination of nutrients fluxes and stresses. Nutrients supply allows the biofilm to grow, while the presence of stresses could inhibit the growth and/or erode the structure. In particular, the local growth is calculated in term of displacements perpendicular to the biofilm surface in the following way ~ ~ d F ¼ d S ¼ K 1 J S K 2 s S K 3 t S on Gg ð0; T bio Þ;
ð19Þ
Here K1, K2, and K3 are constants, while JS is the mass flux through the interface, and sS and tS are the normal and shear stresses on it, respectively. Gg represents the part of the interface G allowed to grow, while Tbio denotes the biological time period. The values of the variables JS, sS, and tS are obtained from the coupled FSI and mass transfer model, either from the last step, when the problem reaches a steadystate solution, or as an average over a fixed period of time, if it reaches a periodic steady state solution, as it happens for flapping streamers (Coroneo et al., 2012). It is important to highlight that in the present model different and more complicated growth equations can be easily introduced and that Equation (19) constitutes only a first simple phenomenological example. In any case, the calculation of the biofilm growth based on the mass flux through the structure surface and on the stresses on it permits to catch the local information. In this way no further hypothesis on the homogeneity of the biofilm structure is needed. Clearly, this
1388
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
is only a first step toward a complete model for biofilm growth. For the purpose of permitting the biofilm to grow and of properly deforming the structure and the fluid mesh, in the growth step an ALE approach is additionally applied to the structure domain (Fig. 1b). For this purpose, the interface displacement due to growth calculated from Equation (19) is applied on the grid structure interface as a Dirichlet boundary condition ~ dG;S ¼ d S
on
G g ð0; T bio Þ;
ð20Þ
while considering fixed all the other boundaries of the structure ALE mesh. Within the domain, the ALE field is treated again as a quasi-elastostatic pseudo-structure (Wall and Ramm, 1998) and an equation of motion similar to Equation (9) is solved on the structure ALE mesh. Afterwards an ALE step is solved again on the fluid ALE mesh, but this time with different boundary conditions. In this step the boundary of the fluid ALE mesh is coupled with the Lagrangian mesh of the growing structure and to an Eulerian mesh at the in- and outflow portions denoted by GE, that is ~ dG;F ¼ d F dG;F ¼ 0
on
on
G g ð0; T bio Þ;
G G;F E ð0; T bio Þ:
ð21Þ ð22Þ
Multi-Scale Coupled Algorithm Monolithic schemes proved to be the most stable and efficient approaches for modeling complex biological problems involving the coupling of incompressible flows to soft structures (Klöppel and Wall, 2010; Küttler et al., 2010). For this reason, the fully coupled non-linear FSI problem is solved monolithically via an advanced multigrid approach (Gee et al., 2010). Moreover, since the mass transfer process does not influence the fluid flow and the structural deformation, the FSI and the transport process are coupled through a one-way coupling. Concerning the growth model, since biological phenomena happen at different time-scales with respect to FSI and scalar transport (Picioreanu et al., 1999), a multi-scale algorithm was developed. It consists of an inner timeloop solving FSI and scalar transport at fluid-dynamic time-scale and an outer timeloop solving the biofilm growth at biological time-scale. A similar idea has also been used previously for cardiovascular problems by Figueroa et al. (2009). As a result, at each fluid dynamic time step, the non-linear equations governing the fluid flow and the structure displacement are solved monolithically till when residuals rFSI reach a problem-specific tolerance eFSI. Subsequently, the local velocities and deformations, obtained from the FSI calculation, are passed to the mass transport subproblem and the dynamic convection-diffusion-reaction equations are solved on both the fluid and solid domains until residuals
rmass meet a problem-specific tolerance emass. These steps are repeated for a total time tmax, until a stationary or a periodic steady state condition is reached. Then, the mass flux normal to the interface and the normal and shear stresses on it are transferred to the growth and erosion subproblem for the calculation of the correct growth magnitude for a longer time step. Afterwards, nodes reference position are changed for all fields, and all the state vectors are reset. The new configuration is then ready for a new inner timeloop without any residual stress or changes in the mechanical properties. The proposed procedure is summarized in Algorithm 1 and Figure 2. It is, however, important to highlight that, when the changes in the structure shape or the structure growth are significant, remeshing could be necessary for the purpose of maintaining a high level of resolution or eliminating excessively distorted elements. Algorithm 1 Multi-scale coupling of FSI and mass transport subproblems to the growth model while tbio < tbio,max do while t < tmax do while k rFSI k> eFSI do Set up and solve the FSI problem Update dS, uF, and dG Calculate sS and tS Check convergence end while Transfer dS, dG,F, uS, uF, and cF to mass transport subproblem while k rmass k> emass do Set up and solve the mass transport subproblem Update FF and FS Calculate JS Check convergence end while Update t end while Transfer sS, tS, and JS to biofilm growth subproblem ~S Calculate d Set up and solve the biofilm growth subproblem Update node reference position Reset state vectors Update tbio end while
Figure 2.
Numerical Examples The validity of the proposed model is demonstrated by selected numerical examples. All the simulations proposed in the present study are based on the implementation of the algorithm discussed above in our in-house research code BACI (Wall and Gee, 2013). The field equations are discretized in space through finite elements and in time through implicit time integration schemes. The numerical solution of the resulting set of non-linear algebraic equations is then obtained through a Newton-type method. In the examples presented in the present work, a one-step-u time integration scheme with u ¼ 0.66 is used for the time discretization of fluid, structure and transport equations. For space discretization in most cases trilinear, hexahedral finite elements are used. Hence, when not differently highlighted, presented examples are based on three-dimensional models and discretizations, although a pseudo two-dimensional deformation and flow state are in some cases enforced by specific boundary conditions. The inner timeloop of the presented approach for coupling FSI and mass transport has already been successfully applied for simulating flapping biofilm streamers (Coroneo et al., 2012) for flow conditions and geometrical details similar to those experimentally investigated by Stoodley et al. (1998). The comparison of mass flux between a flapping and a fixed streamer has demonstrated the importance of coupling FSI and scalar transport for modeling the mass transfer on deformable biofilm structures under real operating and material conditions (Taherzadeh et al., 2012). For these reasons, the numerical examples proposed in the present study will not focus on the reliability and importance of coupling FSI and mass transfer. On the other hand, the following examples will highlight the capability of the model to reproduce the Monod non-linear kinetic and the ability of the multi-scale growth model to tackle simple growing biofilm structures. When quantitative results are of minor importance, spatial and temporal development of concentrations are computed with respect to an arbitrarily chosen prescribed value and the scalar concentration field F
Schematic of the coupling of FSI and mass transport to the growth model.
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1389
is dimensionless. In contrast, since the first and last examples reproduce a real application, all details are reported and the scalar concentration has the real values and dimensions.
through a steady-state convection-reaction equation in one-dimension
Kinetic Models
Simulation results for the two sets of simulations are reported in Figure 3 together with experimental data from Bekins et al. (1998). Results indicate the suitability of our stabilized implementation of the Monod kinetic substrate consumption to correctly predict the experimental behavior of the two investigated systems.
The rate of substrate utilization can be related to its concentration through a zero-order, a first-order, or a nonlinear Monod kinetic. In the three mentioned cases, the reaction rate R reads, respectively R0 ¼ k;
R1 ¼
k F; K
RM ¼ k
F : K þF
ð23Þ
For the purpose of highlighting the accuracy of the approach in catching concentration fields, simulations reproducing experimental data available in literature (Bekins et al., 1998) on biodegradation in contaminated aquifers are run with the three different kinetics. The first set of simulations concerns the methanogenic degradation of phenol in a batch laboratory microcosm, with phenol being the rate-limiting substrate. The substrate utilization rate k and the half saturation constant K are assumed equal to the experimental ones, being respectively, 1.39 mg/L day and 1.7 mg/L, and the system is studied through the solution of a time-dependent reaction equation @F þ R ¼ 0: @t
ð24Þ
The second set of simulations refers to a field site contaminated with creosote compounds in steady-state conditions and with non-transverse dispersion. In this case the experimental substrate utilization rate is k ¼ 0.337 mg/L day, the half saturation constant K ¼1.33 mg/L and the fluid velocity ux ¼ 1 m/day. This system is studied
Figure 3.
ux
@F þ R ¼ 0: @x
ð25Þ
Mass Conservation In this numerical example, the computational domain represents the section of a fluid channel, where on the wall a uniform biofilm layer is present. The computational domain consists of a cube (dimensions 10 mm 10mm 10 mm) in which the lower part (height 2 mm) consists of a flat biofilm structure VS and the upper one is fluid VF (height 8 mm), as reported in Figure 4. The interface between the biofilm and the fluid is denoted by G and represents the interface for FSI, mass transfer and biofilm growth. Dirichlet conditions are applied to the upper and lower boundaries for mass transport, while on all the other boundaries zero-flux conditions are prescribed. This in turn means that no inlet and outlet boundaries are present. As initial condition, a linear distribution of substrate concentration along the height of the domain and a zero velocity distribution are prescribed. With regard to the multi-scale approach, a biological time step 1,000 times larger than the fluid dynamic one was chosen. The substrate concentration distribution at the initial state and at the final step are reported in Figure 4, where the change of the interface G and, as a result, the biofilm growth is visible. In Figure 5 the substrate profiles along a vertical line in the whole domain at different time steps are reported. Because of
(a) Phenol concentration data (Bekins et al., 1998) from methanogenic degradation of phenol in a batch laboratory microcosm (symbols) and simulation results (lines). (b) Field data (Bekins et al., 1998) from a creosote-contaminated site (symbols) and simulation results (lines).
1390
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
Figure 4.
Distribution of substrate concentration at the initial (left) and final step (right).
the presence of the reaction term, the substrate concentration inside the domain diminishes in time, and inside the structure domain reaches values closer to the boundary one. As a consequence, the substrate concentration inside the overall domain is not in equilibrium and produces a mass flux from the fluid to the structure. The presence of the flux through the interface, in turn, produces biofilm growth. In Figure 5 the biofilm growth is also visible, since the dashed lines, referring to the structure domain, are longer at increasing time. Because of the presence of substrate reaction and biofilm growth, a check of the conservation of substrate mass inside the domain is necessary. For this purpose, the change of mass inside the domain, the mass flux through the boundaries, the change of mass due to reaction as well as the overall sum are reported in the following Figure 6. In particular, Figure 6b
Figure 5. Simulated concentration profiles in time along the vertical axis in the fluid (solid lines) and in the structure domains (dashed lines) at different time steps.
refers to the course during the first inner timeloop, and it shows a change in time of the flux through the boundaries and of the mass inside the domain for the first steps. Afterwards the solution reaches a steady-state condition and no more change of mass inside the domain is present anymore. Clearly, the overall sum is always zero and the mass is conserved in the inner timeloop. On the other hand, Figure 6a presents the changes at the end of each inner timeloop before the following growth step. In this case, no change of mass inside the domain is present, since each point represents a steady state solution. As a result, the change in mass due to reaction equals the mass flux out of the domain and the conservation of mass holds true for the multi-scale growth model, too. Influence of Phenomenological Law on Biofilm Growth The final biofilm shape is obviously determined by the values of the constants introduced in Equation (19). As a matter of fact, the presence of all three different terms and their values deeply influences the biofilm growth, while being determined by the particular biofilm under consideration and by its history. For the purpose of highlighting the effect of each contribution to the final biofilm shape, different simulations were run, only taking one term into account at a time. For these simulations a flow channel (dimensions 6 mm 9 mm 0.2 mm) was considered, inside which a circular biofilm structure VS (radius 0.5 mm and thickness 0.2 mm) is positioned. A constant velocity and concentration is imposed at the left boundary, while no-slip conditions are imposed at the top and bottom boundaries, in order to not develop any additional boundary layer. The scalar distribution in the fluid domain at the end of the first inner timeloop is reported in Figure 7. In this condition, the mass flux in the front part of the structure is higher than in its back part. For this reason if the growth is calculated only taking in account the mass flux, i.e. only the constant K1, the biofilm grows on all the surface but more on the front of it, as reported in Figure 7b. Also
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1391
Figure 6. P R ID
Time-dependent progress at the end of each inner timeloop in overall mass within the fluid and structure domain (a) and in the first timeloop (b) of the change R P FI dVI =Dt, of the mass flux out of the domain I FI c DI rFI nI ; of the change of mass due to reaction RS dVS and of the sum of all the previous terms.
normal stresses are higher in the front of the biofilm, producing an erosion of the structure, as is clearly visible when only the term relative to erosion due to normal stresses is taken into account (Fig. 7c). With regard to shear stresses, it has to be taken in account that their distribution is not uniform on all the front part of the structure and in particular that they are zero exactly on the mid point of the front part. A growth law, which only allows for the contribution of shear stresses, will obviously calculate a zero growth in that position, while producing erosion on the sides. This effect will produce a non-smooth structure shape, as reported in Figure 7d. The equation of growth (19), however, will never take in account only the contribution due to shear stresses, since in the reality the effect due to mass flux will always be present. For this reason, the highlighted problem will usually not be present. Effect of Operating Conditions on Biofilm Growth The model is able to predict the expected behavior not only when the structure is exposed to uniform fluxes but also when different concentration gradients are present. Actually, the operating conditions are known to have a big influence on biofilm growth and this effect is highlighted in this example. For this purpose a fluid domain VF (dimensions 10 mm 8 mm 0.2 mm) bound by a biofilm structure VS (dimensions 10 mm 2 mm 0.2 mm) at the bottom is taken into account. Different boundary conditions were applied to the same domain, while using constant material parameters. The interface G between the biofilm and the fluid represents the interface for FSI and mass transfer, as well as for biofilm growth. Only the interface G is allowed to move in any direction, while all the other boundaries are kept fixed. When not differently specified a zero velocity field is also applied at the fluid boundaries and a zero-flux conditions at the solid and fluid ones. As an initial condition, a zero velocity field is also always prescribed.
1392
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
Three different case were considered. (i) A constant unitary concentration is applied at the top of the fluid domain. (ii) A cosine scalar profile with a unitary maximum value at its center and zero-values at the sides is applied at the top boundary of the fluid domain. (iii) Conditions similar to the case (ii) were used, except for the fact that also a non-zero velocity x-component was applied at the left fluid boundary, with a parabolic profile in the vertical direction, a null value on G and a maximum value of 0.1 mm/ms on the top. In this case also flow through the right fluid boundary is allowed and a slip boundary condition is prescribed at the top boundary. Also for this example a biological time step 1000 times larger than the fluid dynamic one was used. In Figure 8 the scalar distribution in the whole domain and the fluidstructure interfaces are reported for the three different cases and at different biological time steps. The influence of the different boundary conditions on the final biofilm shape is evident. Moreover, when the scalar distribution is uniform, then a uniform growth is present (i), for not-uniform scalar fields, a larger growth is present in the direction of the higher scalar concentration (ii), while the presence of the velocity field distorts the scalar distribution producing a non-symmetric final biofilm shape (iii). FSI and Growth of a Finger-Like Biofilm Structure In this example, the proposed approach is applied to the simulation of a biofilm finger under realistic operating conditions and material properties. For this purpose the geometry reported in Figure 9 is taken into account. The computational domain represents a part of a flow channel (dimensions 0.6 mm 0.3 mm 0.01 mm) inside which a
Figure 7. Scalar distribution around a circular biofilm structure (top) and initial (a) and final biofilm shape superimposed to the initial one for K1 6¼ 0 and K2 ¼ K3 ¼ 0 (b), for K2 6¼ 0 and K1 ¼ K3 ¼ 0 (c), for K3 6¼ 0 and K1 ¼ K2 ¼ 0 (d).
Figure 8.
Concentration distribution in the fluid and structure domains for the three different cases after 0, 20, and 40 biological time steps.
finger-like structure is present (0.1 mm high and 0.04 mm wide). The simulation is run under the realistic operating conditions and material parameters summarized in Table I for the cases of rigid and deformable biofilm. At the inlet boundary parabolic velocity and concentration profiles are imposed, while at the top boundary and at the bottom one a slip and a no-slip condition for the velocity are prescribed respectively. As initial conditions a zero velocity field and a parabolic concentration profile in the vertical direction are imposed. The simulation was run for an hydrodynamic timestep of 5 ms and a biological timestep of 1day. The first FSI and mass transport timeloop was run for 500 steps until a first steady state solution was obtained. In the following inner timeloops, the velocity and scalar distribution calculated in the last inner timeloop, before the growth step, were imposed as initial conditions. In this way, the number of time steps in the following inner timeloops to reach a steady state solution was reduced to 20. In Figure 9 the initial biofilm shape and the one after 8 biological timesteps, that is, after 8 days, is reported for the two different cases. The different final biofilm shape for a fixed structure (Fig. 1b) and a deformable one (Fig. 1c) is apparent, demonstrating that the growth can be affected by the deformation and highlighting the importance of the proposed approach. The boundary and the operating conditions chosen in the present example, together with the
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1393
prescribed growth constants, produce a larger growth of the finger in the width direction instead of in height and the final finger results thicker but not much higher than the initial one. For different boundary and operating conditions, material parameters and biofilm growth constants, different biofilm shapes can be obtained. For these reasons, the results reported in Figure 9 only show a possible final biofilm shape. The present example demonstrates, however, the applicability of the present approach to the simulation of realistic biofilm structure under realistic material and operating conditions and the importance of the proposed coupled approach.
Conclusions
Figure 9. Definition of the computational domain (top), initial (a) and final shape of the finger-like biofilm, in the case of fixed structure (b), deformable structure in undeformed (c), and deformed configuration (d).
Table I. Material parameters, operating conditions and growth constants. Parameter Liquid Dynamic viscositya Densitya Inlet velocity Biofilm Densityb Young’s modulusc Poisson ratiob Substrate Diffusion coefficientd Uptake rate coefficiente Saturation coefficientf Inlet concentrationg Growth Normal flux constant Normal stress constant Shear stress constant
Symbol
Value
Unit
mF rF uFIN
103 103 1.5 102
kg m1 s1 kg m3 m s2
rS ES nS
3 102 102 4 101
kg m4 kg m1 s2
D k K FFIN
2.5 109 3 102 3 103 2.5 102
m2 s1 mol m3 s1 mol m3 mol m3
K1 K2 K3
4 105 102 1
s m3/mol s m2/kg s m2/kg
Water at 20 C. TanyolaSc and Beyenal (1997). c Böl et al. (2013). d Dissolved oxygen in water at 20 C. e Assumed for oxygen uptake by heterotrophic microorganisms (Henze, 2000). f Affinity for oxygen of heterotrophic microorganisms (Henze, 2000). g Relatively low dissolved substrate concentration. a
b
1394
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
To enable a better understanding of biofilm growth and development, we developed a novel and advanced computational model for the numerical simulation of FSI and mass transfer of growing biofilm structures. For the purpose of describing correctly and in enough detail the different phenomena involved, the fully coupled nonlinear FSI problem and the nonlinear multi-field mass transport model are encapsulated in the biofilm growth model through a multi-scale approach in time. The described new methodology enables the mutual coupling of (i) flow velocity and solid deformation (ii), mass transport on deforming fluid and solid domains, (iii) as well as growth, based on local mass transport and stresses. The proposed model takes in account the effect of nutrients flux and stresses on the biofilm surface through a simple phenomenological equation. Although the effect of mass transfer on biofilm growth as well as the influence of structure deformation on mass transport have already been considered before, none of the previous approaches were able to take into account all these important phenomena at the same time and without any restrictions. Furthermore, considering FSI, mass transport and growth at the same time is essential for understanding the formation and development of biofilm structures like streamers. The general validity of the proposed model was illustrated by selected numerical examples. The accurate implementation of the non-linear Monod kinetic was highlighted through comparison with experimental data. The mass conservation was shown. The effect of phenomenological law constants on biofilm growth was highlighted. The influence of the different operating conditions was discussed. The applicability of the novel approach to a realistic structure and at realistic operating conditions and material parameters was reported. As a result, selected examples successfully demonstrated the general suitability of the proposed model for modeling the FSI and the mass transfer of growing biofilm structures. After an appropriate calibration of the phenomenological law driven by experimental data and evidence, the application of this coupled and multi-scale approach to different real biofilm structures will be explored. Future work will also concern the introduction of a detachment model and its
application to new realistic scenarios also from different fields of biomechanics. References Alpkvist E, Klapper I. 2007a. A multidimensional multispecies continuum model for heterogeneous biofilm development. Bull Math Biol 69: 765–789. Alpkvist E, Klapper I. 2007b. Description of mechanical response including detachment using a novel particle model of biofilm/flow interaction. Water Sci Technol 55:265–273. Battin TJ, Kaplan LA, Newbold JD, Cheng X, Hansen C. 2003. Effects of current velocity on the nascent architecture of stream microbial biofilms. Appl Environ Microbiol 69:5443–5452. de Beer D, Schramm A. 1999. Micro-environments and mass transfer phenomena in biofilms studied with microsensors. Water Sci Technol 39:173–178. Bekins BA, Warren E, Godsy EM. 1998. A comparison of zero-order, firstorder, and Monod biotransformation models. Ground Water 36: 261–268. Böl M, Möhle RB, Haesner M, Neu TR, Horn H, Krull R. 2009. 3D finite element model of biofilm detachment using real biofilm structures from CLSM data. Biotechnol Bioeng 103:177–186. Böl M, Ehret AE, Bolea Albero A, Hellriegel J, Krull R. 2013. Recent advances in mechanical characterisation of biofilm and their significance for material modelling. Crit Rev Biotechnol 33:145–171. Cogan NG. 2008. Two-fluid model of biofilm disinfection. Bull Math Biol 70:800–819. Cogan NG. 2011. Computational exploration of disinfection of bacterial biofilms in partially blocked channels. Int J Numer Method Biomed Eng 27:1982–1995. Coroneo M, Yoshihara L, Bauer G, Wall WA. 2012. A finite element approach for the coupled numerical simulation of fluid-structure interaction and mass transfer of moving biofilm structures In: Eberhardsteiner J, et al. editor European Congress on Comp Meth Appl Sci Eng, ECCOMAS. Dillon R, Fauci L, Fogelson A. 1996. Modeling biofilm processes using the immersed boundary method. J Comput Phys 129:57–73. Dockery J, Klapper I. 2001. Finger formation in biolm layers. SIAM J Appl Math 62:853–869. Donea J, Huerta A, Ponthot JP, Rodriguez-Ferran A. 2004. Arbitrary Lagrangian–Eulerian methods. In: Stein E, de Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics, Vol. 1: Fundamentals. Chichester, West Sussex: John Wiley Sons. Duddu R, Bordas S, Chopp D, Moran B. 2008. A combined extended finite element and level set method for biofilm growth. Int J Numer Method Eng 74:848–870. Duddu R, Chopp DL, Moran B. 2009. A two-dimensional continuum model of biofilm growth incorporating fluid flow and shear stress based detachment. Biotechnol Bioeng 103:92–104. Eberl HJ, Sudarsan R. 2008. Exposure of biofilms to slow flow fields: The convective contribution to growth and disinfection. J Theor Biol 253:788–807. Figueroa CA, Baek S, Taylor CA, Humphrey JD. 2009. A computational framework for fluid-solid-growth modeling in cardiovascular simulations. Comput Methods Appl Mech Eng 198:3583–3602. Garny K, Horn H, Neu T. 2008. Interaction between biofilm development, structure and detachment in rotating annular reactors. Bioprocess Biosyst Eng 31:619–629. Gee MW, Küttler U, Wall WA. 2010. Truly monolithic algebraic multigrid for fluid-structure interaction. Int J Numer Method Eng 85:987–1016.
Graf von der Schulenburg DA, Pintelon TR, Picioreanu C, Van Loosdrecht MCM, Johns ML. 2009. Three-dimensional simulations of biofilm growth in porous media. AIChE J 55:494–504. Hall-Stoodley L, Stoodley P. 2002. Developmental regulation of microbial biofilms. Curr Opin Biotechnol 13:228–233. Henze M. 2000. Activated sludge models ASM1, ASM2, ASM2d, and ASM3. London, UK: IWA Publishing. p 128. Holzapfel G. 2000. Nonlinear solid mechanics. A continuum approach for engineering. Chichester, UK: Wiley. p 455. Klapper I. 2004. Effect of heterogeneous structure in mechanically unstressed biofilms on overall growth. Bull Math Biol 66:809–824. Klöppel T, Wall WA. 2010. A novel two-layer, coupled finite element approach for the nonlinear elastic and viscoelastic behavior of human erythrocytes. Biomech Model Mechanobiol 10:445–459. Küttler U, Gee M, Förster Ch, Comerford A, Wall WA. 2010. Coupling strategies for biomedical fluid-structure interaction problems. Int J Numer Method Biomed Eng 26:305–321. Pereira MO, Kuehn M, Wuertz S, Neu T, Melo LF. 2002. Effect of flow regime on the architecture of a Pseudomonas fluorescens biofilm. Biotechnol Bioeng 78:164–171. Picioreanu C, van Loosdrecht MCM, Heijnen J. 1999. A discrete-differential modelling of biofilm structure. Water Sci Technol 39:115–122. Radu AI, Vrouwenvelder JS, van Loosdrecht MCM, Picioreanu C. 2012. Effect of flow velocity, substrate concentration and hydraulic cleaning on biofouling of reverse osmosis feed channels. Chem Eng J 188:30–39. Simões M, Pereira M, Sillankorva S, Azeredo J, Vieira MJ. 2007. The effect of hydrodynamic conditions on the phenotype of Pseudomonas fluorescens biofilms. Biofouling 23:249–258. Smith BG, Vaughan JR BL, Chopp DL. 2007. The extended finite element method for boundary layer problems in biofilm growth. Commun Appl Math Comput Sci 2:35–56. Stoodley P, Lewandowski Z, Boyle J, Lappin-Scott HM. 1998. Oscillation characteristics of biofilm streamers in turbulent flowing water as related to drag and pressure drop. Biotechnol Bioeng 57:536–544. Stoodley P, Dodds I, Boyle J, Lappin-Scott H. 1999. Influence of hydrodynamics and nutrients on biofilm structure. J Appl Microbiol 85:19S–28S. Stoodley P, Sauer K, Davies DG, Costerton1 JW. 2002. Biofilms as complex differentiated communities. Annu Rev Microbiol 56:187–209. Taherzadeh D, Picioreanu C, Küttler U, Simone A, Wall WA, Horn H. 2010. Computational study of the drag and oscillatory movement of biofilm streamers in fast flows. Biotechnol Bioeng 105:600–610. Taherzadeh D, Picioreanu C, Horn H. 2012. Mass transfer enhancement in moving biofilm structures. Biophys J 102:1483–1492. TanyolaSc A, Beyenal H. 1997. Prediction of average biofilm density and performance of a spherical bioparticle under substrate inhibition. Biotechnol Bioeng 56:319–329. Teodósio JS, Simões M, Alves MA, Melo LF, Mergulh~ao FJ. 2012. Setup and validation of flow cell systems for biofouling simulation in industrial settings. Sci World J 12. Towler BW, Cunningham A, Stoodley P, McKittrick L. 2007. A model of fluidbiofilm interaction using a Burger material law. Biotechnol Bioeng 96:259–271. Wall WA, Ramm E. 1998. Fluid-structure interaction based upon a stabilized (ALE) finite element method In: Idelsohn SR, Onate E, Dvorkin EN, editors. Computational mechanics—New trends and applications. Barcelona, Spain: CIMNE. pp 1–20. Wall WA, Gee MW. 2013. BACI—A multiphysics simulation environment. Munich, Germany: Technische Universität München. Wijeyekoon S, Mino T, Satoh H, Matsuo T. 2004. Effects of substrate loading rate on biofilm structure. Water Res 38:2479–2488.
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1395
Introduction ABSTRACT: In this paper, we propose a novel approach for modelling biofilm growth. It is based on a finite element method and includes both fluid–structure interaction (FSI) as well as scalar transport effects. Due to the different timescales of the involved phenomena, the growth of the biofilm structure is coupled with the FSI and mass transport through a multi-scale approach in time. In each hydrodynamic time step, first the non-linear FSI problem is solved followed by the scalar transport equations, using the information on velocities and deformations obtained in the FSI step. After a steady state solution is reached, information on mass fluxes and stresses are passed to the growth model. At this point, the growth is calculated for a biological time step larger than the hydrodynamic one and based on the mass flux through the interface and on normal and shear stresses on it. This type of approach can significantly contribute to the understanding of biofilm development in fluid flows, since the influence of hydrodynamic conditions and availability of nutrients is well known to have effects on biofilm development. Therefore, for the purpose of understanding biofilm macro-scale dynamics, it is essential to adopt a modeling approach, which takes into account all the relevant aspects, like fluid flow, structure deformation, mass transport and their effect on biofilm growth and erosion. First numerical examples demonstrate the suitability of the proposed model to catch the main features of a growing biofilm structure. Biotechnol. Bioeng. 2014;111: 1385–1395. ß 2014 Wiley Periodicals, Inc. KEYWORDS: biofilm growth; fluid–structure interaction; mass transfer; finite element method
Correspondence to: M. Coroneo Received 25 September 2013; Revision received 3 January 2014; Accepted 6 January 2014 Article first published online 4 February 2014 in Wiley Online Library (http://onlinelibrary.wiley.com/doi/10.1002/bit.25191/abstract). DOI 10.1002/bit.25191
ß 2014 Wiley Periodicals, Inc.
A better understanding of biofilm growth and behavior is essential for the development of successful strategies for biofilm control (Böl et al., 2013), that delay the formation of harmful biofilms or promote the formation of beneficial ones (Teodósio et al., 2012). The formation and development of biofilms, however, is influenced by many factors, including nutrient availability and hydrodynamic conditions (Battin et al., 2003; Garny et al., 2008; Pereira et al., 2002; Stoodley et al., 1999; Wijeyekoon et al., 2004). As a matter of fact, hydrodynamic conditions determine not only the rate at which nutrients are transported to the biofilm surface, but also the magnitude of forces acting on a growing biofilm (Simões et al., 2007), which will also have an influence on the biofilm structure (Hall-Stoodley and Stoodley, 2002). Understanding the mutual interaction of all the involved phenomena, however, has proven to be difficult, because of the many difficulties in carrying out experiments capable of isolating each single mechanism. For this reason, biofilm modeling is arising as a means of producing quantitative tools (Duddu et al., 2009). Biofilms, however, exhibit heterogeneous morphologies (de Beer Schramm, 1999) and the complexity and variability of their matrices impact their physical properties, also significantly increasing the difficulty of biofilm modeling (Pereira et al., 2002). The way they respond to environmental solicitations, however, is defined by their macro-scale characteristics. Indeed biofilms interact with the surrounding fluid as macro-scale materials and can be studied as flexible structures located in liquid flows. For this reason, continuum models, assuming that local properties can be considered smeared on all the spatial extension of the biofilm, can be used and are, moreover, the most suitable for studying big biofilm aggregates. The complex interaction between fluid flow, structure deformation, mass transfer and biofilm growth is fundamental for biofilm development. Most continuum models applied to biofilm modeling, however, neglect fluid-induced deformation of the biofilm structure (Alpkvist and Klapper, 2007a; Dockery and Klapper, 2001; Duddu et al., 2008, 2009;
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
1385
Eberl and Sudarsan, 2008; Graf von der Schulenburg et al., 2009; Radu et al., 2012; Smith et al., 2007), while only few account for it (Alpkvist and Klapper, 2007b; Böl et al., 2009; Taherzadeh et al., 2010; Towler et al., 2007), or for both the substrate transport and the biofilm transient interaction with the surrounding fluid (Cogan, 2008, 2011; Dillon et al., 1996; Klapper, 2004; Taherzadeh et al., 2012). Most of these studies are, moreover, restricted to twodimensions or to specific domain geometries (Cogan, 2011), disregard the advective effect (Alpkvist and Klapper, 2007a; Klapper, 2004), do not take into account the growth (Dillon et al., 1996; Taherzadeh et al., 2012) or calculate it only based on the concentration field (Dockery and Klapper, 2001). For biofilm streamers, for example, calculating the fluid-induced deformations is essential for the purpose of correctly calculating the mass transfer (Taherzadeh et al., 2012) and consequently the growth. For these structures it is also important to apply a multi-scale approach for the purpose of obtaining appropriate and detailed enough models for both the fluid flow and the mass transfer, as well as the biofilm growth. To date a general multi-scale modeling approach taking into account the effect of both local structure deformation and fluid flow on mass transfer and growth is still missing. In this work, we propose a novel multi-scale and coupled fluid-structure interaction (FSI) and mass transport approach for modeling biofilm growth. It is developed in our in-house research code BACI and is aimed at understanding the influence of operating conditions on biofilm growth and erosion and at developing biofilm control strategies.
Fluid–Structure Interaction The interaction between the biofilm structure and the surrounding fluid is modeled through an arbitrary Lagrangian–Eulerian (ALE) approach (Donea et al., 2004). As a result, the FSI problem domain consists of three fields. Apart from the structure VS and the fluid VF domains, sharing a common interface G, a third non-physical mesh field VG,F is present due to the ALE approach (Fig. 1a). Structure Field Biofilms are viscoelastic materials. They behave however elastically, when subjected to external forces over short periods of time, and in a viscous way, when external forces are applied over longer time period (Stoodley et al., 2002). In this work the interaction between fluid and structure is considered only at very small time scales and therefore the biofilm motion can be safely described by an elastic material model. In the structural domain the following non-linear structural elastodynamics equation is solved for the displacement field dS rS
d2 dS ¼ r ðF SÞ þ rS bS dt 2
in
VS ð0; TÞ;
ð1Þ
where rS represents the structural density, bS are the external body forces, VS denotes the undeformed domain, while T is
1386
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
Figure 1.
Definition of domains and their boundaries for the FSI (a) and biofilm growth (b) subproblems. The boundaries of solid and fluid domains consist of Dirichlet, Neumann, and FSI boundary conditions, that is @VI ¼ G ID [ G IN [ G with I 2 (F,S). The boundaries of VG,I are given by @VG;F ¼ G G;F E [ G for the FSI subproblem and by @VG;I ¼ G G;I E [ G g for the biofilm growth subproblem.
the hydrodynamic time period. The equation determines the structural displacements dS by prescribing an equilibrium between forces of inertia, body forces bS, and the internal forces determined from the second Piola–Kirchhoff stress tensor S and the deformation gradient F. The problem is then completed by the imposition of appropriate initial and boundary conditions. For details on nonlinear solid mechanics including stress measures the reader is referred to the monograph of Holzapfel (2000). Fluid Field In the fluid domain the following incompressible timedependent ALE version of Navier–Stokes equations are solved with the purpose of obtaining the pressure pF and the velocity uF fields rF
@uF þ rF ðcF rÞuF 2mr eðuF Þ @t þ rpF ¼ rF bF in VF ð0; TÞ; r uF ¼ 0
in
VF ð0; TÞ:
ð2Þ
ð3Þ
In the momentum Equation (2), e(uF) represents the strain rate tensor of the fluid, m denotes its dynamic viscosity, bF a prescribed body force, and cF the fluid ALE convective velocity, which represents the fluid velocity uF relative to the arbitrarily moving fluid domain cF ¼ uF uG;F :
ð4Þ
at the fluid and structure surfaces
The fluid grid velocity uG,F is defined by uG;F ¼
@w @t
VF ð0; TÞ;
in
ð5Þ
where w is an arbitrary and unique mapping for the deformation of the fluid domain dG,F d
G;F
ðx; tÞ ¼ w dG;F for G ; x; t
ðx; tÞ 2 V ð0; TÞ; ð6Þ F
while dG;F G represents the mesh interface displacement, which will be later related to the structure interface displacement dSG . Also in this case, appropriate initial and boundary conditions have to be set. It is moreover important to highlight that the fluid flow in the deforming fluid domain is considered to be laminar, since biofilm structures usually grow attached to surfaces and are considered to be in the hydrodynamic boundary layer. Hence, Navier–Stokes equations can be solved without any further (turbulence) modelling and with a reasonable computational effort. ALE Field For the purpose of calculating the fluid ALE convective velocity cF appearing in Equation (2), a definition of the mapping w introduced in Equation (6) has to be determined. In the FSI problems studied in the present work, the boundary of the fluid ALE mesh is coupled to the Lagrangian mesh of the structure and to an Eulerian mesh at the in- and outflow portions denoted by GE, that is, dG;F ¼ dS dG;F ¼ 0
on on
G ð0; TÞ;
ð7Þ
G G;F E ð0; TÞ:
ð8Þ
On the other hand, the fluid ALE mesh is allowed to deform arbitrarily within the domain and, in the present study, it was decided to treat the fluid ALE field as a quasielastostatic pseudo-structure (Wall and Ramm, 1998). The ALE equation of motion results then r s G;F ¼ 0
VG;F ð0; TÞ;
ð9Þ
s G;F ¼ lG;F tr eG;F I þ 2mG;F eG;F ; 1 G;F G;F T rd þ rd eG;F ¼ ; 2
ð10Þ
in
with
where tr(·) denotes the trace operator while lG,F and mG,F represent the Lamé constants of the pseudo-structure. FSI Coupling The deformation of the FSI surface G is controlled by kinematic and dynamic constraints. First, the fulfillment of the equilibrium of forces results in equal surfaces tractions h
hSG ¼ hFG
on
G ð0; TÞ:
ð11Þ
Second, at the interface the fluid grid velocity uG;F G and the fluid velocity uFG also have to match uFG ¼ uG;F G
on
G ð0; TÞ:
ð12Þ
Finally, a mass flow across G as well as a relative tangential movement of fluid and structure at G are prohibited, resulting in a fluid velocity uFG equal to the structure deformation rate, that is @dSG ¼ uFG @t
on
G ð0; TÞ:
ð13Þ
Scalar Transport On both the fluid and solid domains, the scalar transport was also solved. For this purpose, the water-substrate solution is considered diluted and in the fluid domain VF a convectiondiffusion equation is solved for the concentration FF @FF þ cF rFF r ðDF rFF Þ ¼ 0 @t in VF ð0; TÞ;
ð14Þ
where cF is the ALE convective velocity, which has already been introduced in Equation (4), and DF is the fluid diffusion coefficient. On the other hand, in the solid domain VS a diffusion-reaction equation in the conservative form is solved dFS þ FS ðr uS Þ r ðDS rFS Þ þ RS ¼ 0 dt in VS ð0; TÞ:
ð15Þ
In the previous equation, uS is the solid velocity, DS represents the solid diffusion coefficient, while RS denotes the reaction term. The microbial growth rate in aqueous environments is a function of the concentration of a limiting nutrient and is usually expressed by the Monod kinetic. If this kinetic is assumed to hold then the substrate consumption rate RS reads RS;M ¼ k
FS ; K þ FS
ð16Þ
where k and K represent the reaction rate constant and the half-saturation concentration of substrate, respectively. Furthermore, it is important to emphasize that, in contrast to the scalar equation for the fluid (14), in Equation (15) the convective velocity is not present, since the structure velocity S uS ¼ @d @t and the grid velocity are identical. Moreover, a conservative formulation is chosen for the mass transport in the solid since in general 5 · uS 6¼ 0.
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1387
The convection-diffusion transport equation in the liquid domain is coupled to the diffusion-reaction transport equation in the structural domain, with the constraints of equal fluxes h and concentrations F at the fluid-structure interface FS ¼ FF hS ¼ hF
on on
G ð0; TÞ; G ð0; TÞ:
ð17Þ ð18Þ
The initial boundary value problem is then completed by appropriate boundary conditions and initial concentration fields. Through the solution of the reported equations, the transport of the substrate from the bulk of the fluid till and through the fluid-structure interface and then its transport and reaction inside the structure is properly modeled. Biofilm Growth The formation and development of different biofilm structures are influenced by flow conditions and nutrients availability (Stoodley et al., 1999). As a result, a reliable model of biofilm growth needs to take into account variables not only connected to mass transfer but also to the flow conditions. The proposed biofilm model calculates the local amount of growth or erosion based on stresses and mass fluxes at the interface resulting from the FSI and mass transfer subproblems. The underlying hypothesis is that growth is the result of the combination of nutrients fluxes and stresses. Nutrients supply allows the biofilm to grow, while the presence of stresses could inhibit the growth and/or erode the structure. In particular, the local growth is calculated in term of displacements perpendicular to the biofilm surface in the following way ~ ~ d F ¼ d S ¼ K 1 J S K 2 s S K 3 t S on Gg ð0; T bio Þ;
ð19Þ
Here K1, K2, and K3 are constants, while JS is the mass flux through the interface, and sS and tS are the normal and shear stresses on it, respectively. Gg represents the part of the interface G allowed to grow, while Tbio denotes the biological time period. The values of the variables JS, sS, and tS are obtained from the coupled FSI and mass transfer model, either from the last step, when the problem reaches a steadystate solution, or as an average over a fixed period of time, if it reaches a periodic steady state solution, as it happens for flapping streamers (Coroneo et al., 2012). It is important to highlight that in the present model different and more complicated growth equations can be easily introduced and that Equation (19) constitutes only a first simple phenomenological example. In any case, the calculation of the biofilm growth based on the mass flux through the structure surface and on the stresses on it permits to catch the local information. In this way no further hypothesis on the homogeneity of the biofilm structure is needed. Clearly, this
1388
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
is only a first step toward a complete model for biofilm growth. For the purpose of permitting the biofilm to grow and of properly deforming the structure and the fluid mesh, in the growth step an ALE approach is additionally applied to the structure domain (Fig. 1b). For this purpose, the interface displacement due to growth calculated from Equation (19) is applied on the grid structure interface as a Dirichlet boundary condition ~ dG;S ¼ d S
on
G g ð0; T bio Þ;
ð20Þ
while considering fixed all the other boundaries of the structure ALE mesh. Within the domain, the ALE field is treated again as a quasi-elastostatic pseudo-structure (Wall and Ramm, 1998) and an equation of motion similar to Equation (9) is solved on the structure ALE mesh. Afterwards an ALE step is solved again on the fluid ALE mesh, but this time with different boundary conditions. In this step the boundary of the fluid ALE mesh is coupled with the Lagrangian mesh of the growing structure and to an Eulerian mesh at the in- and outflow portions denoted by GE, that is ~ dG;F ¼ d F dG;F ¼ 0
on
on
G g ð0; T bio Þ;
G G;F E ð0; T bio Þ:
ð21Þ ð22Þ
Multi-Scale Coupled Algorithm Monolithic schemes proved to be the most stable and efficient approaches for modeling complex biological problems involving the coupling of incompressible flows to soft structures (Klöppel and Wall, 2010; Küttler et al., 2010). For this reason, the fully coupled non-linear FSI problem is solved monolithically via an advanced multigrid approach (Gee et al., 2010). Moreover, since the mass transfer process does not influence the fluid flow and the structural deformation, the FSI and the transport process are coupled through a one-way coupling. Concerning the growth model, since biological phenomena happen at different time-scales with respect to FSI and scalar transport (Picioreanu et al., 1999), a multi-scale algorithm was developed. It consists of an inner timeloop solving FSI and scalar transport at fluid-dynamic time-scale and an outer timeloop solving the biofilm growth at biological time-scale. A similar idea has also been used previously for cardiovascular problems by Figueroa et al. (2009). As a result, at each fluid dynamic time step, the non-linear equations governing the fluid flow and the structure displacement are solved monolithically till when residuals rFSI reach a problem-specific tolerance eFSI. Subsequently, the local velocities and deformations, obtained from the FSI calculation, are passed to the mass transport subproblem and the dynamic convection-diffusion-reaction equations are solved on both the fluid and solid domains until residuals
rmass meet a problem-specific tolerance emass. These steps are repeated for a total time tmax, until a stationary or a periodic steady state condition is reached. Then, the mass flux normal to the interface and the normal and shear stresses on it are transferred to the growth and erosion subproblem for the calculation of the correct growth magnitude for a longer time step. Afterwards, nodes reference position are changed for all fields, and all the state vectors are reset. The new configuration is then ready for a new inner timeloop without any residual stress or changes in the mechanical properties. The proposed procedure is summarized in Algorithm 1 and Figure 2. It is, however, important to highlight that, when the changes in the structure shape or the structure growth are significant, remeshing could be necessary for the purpose of maintaining a high level of resolution or eliminating excessively distorted elements. Algorithm 1 Multi-scale coupling of FSI and mass transport subproblems to the growth model while tbio < tbio,max do while t < tmax do while k rFSI k> eFSI do Set up and solve the FSI problem Update dS, uF, and dG Calculate sS and tS Check convergence end while Transfer dS, dG,F, uS, uF, and cF to mass transport subproblem while k rmass k> emass do Set up and solve the mass transport subproblem Update FF and FS Calculate JS Check convergence end while Update t end while Transfer sS, tS, and JS to biofilm growth subproblem ~S Calculate d Set up and solve the biofilm growth subproblem Update node reference position Reset state vectors Update tbio end while
Figure 2.
Numerical Examples The validity of the proposed model is demonstrated by selected numerical examples. All the simulations proposed in the present study are based on the implementation of the algorithm discussed above in our in-house research code BACI (Wall and Gee, 2013). The field equations are discretized in space through finite elements and in time through implicit time integration schemes. The numerical solution of the resulting set of non-linear algebraic equations is then obtained through a Newton-type method. In the examples presented in the present work, a one-step-u time integration scheme with u ¼ 0.66 is used for the time discretization of fluid, structure and transport equations. For space discretization in most cases trilinear, hexahedral finite elements are used. Hence, when not differently highlighted, presented examples are based on three-dimensional models and discretizations, although a pseudo two-dimensional deformation and flow state are in some cases enforced by specific boundary conditions. The inner timeloop of the presented approach for coupling FSI and mass transport has already been successfully applied for simulating flapping biofilm streamers (Coroneo et al., 2012) for flow conditions and geometrical details similar to those experimentally investigated by Stoodley et al. (1998). The comparison of mass flux between a flapping and a fixed streamer has demonstrated the importance of coupling FSI and scalar transport for modeling the mass transfer on deformable biofilm structures under real operating and material conditions (Taherzadeh et al., 2012). For these reasons, the numerical examples proposed in the present study will not focus on the reliability and importance of coupling FSI and mass transfer. On the other hand, the following examples will highlight the capability of the model to reproduce the Monod non-linear kinetic and the ability of the multi-scale growth model to tackle simple growing biofilm structures. When quantitative results are of minor importance, spatial and temporal development of concentrations are computed with respect to an arbitrarily chosen prescribed value and the scalar concentration field F
Schematic of the coupling of FSI and mass transport to the growth model.
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1389
is dimensionless. In contrast, since the first and last examples reproduce a real application, all details are reported and the scalar concentration has the real values and dimensions.
through a steady-state convection-reaction equation in one-dimension
Kinetic Models
Simulation results for the two sets of simulations are reported in Figure 3 together with experimental data from Bekins et al. (1998). Results indicate the suitability of our stabilized implementation of the Monod kinetic substrate consumption to correctly predict the experimental behavior of the two investigated systems.
The rate of substrate utilization can be related to its concentration through a zero-order, a first-order, or a nonlinear Monod kinetic. In the three mentioned cases, the reaction rate R reads, respectively R0 ¼ k;
R1 ¼
k F; K
RM ¼ k
F : K þF
ð23Þ
For the purpose of highlighting the accuracy of the approach in catching concentration fields, simulations reproducing experimental data available in literature (Bekins et al., 1998) on biodegradation in contaminated aquifers are run with the three different kinetics. The first set of simulations concerns the methanogenic degradation of phenol in a batch laboratory microcosm, with phenol being the rate-limiting substrate. The substrate utilization rate k and the half saturation constant K are assumed equal to the experimental ones, being respectively, 1.39 mg/L day and 1.7 mg/L, and the system is studied through the solution of a time-dependent reaction equation @F þ R ¼ 0: @t
ð24Þ
The second set of simulations refers to a field site contaminated with creosote compounds in steady-state conditions and with non-transverse dispersion. In this case the experimental substrate utilization rate is k ¼ 0.337 mg/L day, the half saturation constant K ¼1.33 mg/L and the fluid velocity ux ¼ 1 m/day. This system is studied
Figure 3.
ux
@F þ R ¼ 0: @x
ð25Þ
Mass Conservation In this numerical example, the computational domain represents the section of a fluid channel, where on the wall a uniform biofilm layer is present. The computational domain consists of a cube (dimensions 10 mm 10mm 10 mm) in which the lower part (height 2 mm) consists of a flat biofilm structure VS and the upper one is fluid VF (height 8 mm), as reported in Figure 4. The interface between the biofilm and the fluid is denoted by G and represents the interface for FSI, mass transfer and biofilm growth. Dirichlet conditions are applied to the upper and lower boundaries for mass transport, while on all the other boundaries zero-flux conditions are prescribed. This in turn means that no inlet and outlet boundaries are present. As initial condition, a linear distribution of substrate concentration along the height of the domain and a zero velocity distribution are prescribed. With regard to the multi-scale approach, a biological time step 1,000 times larger than the fluid dynamic one was chosen. The substrate concentration distribution at the initial state and at the final step are reported in Figure 4, where the change of the interface G and, as a result, the biofilm growth is visible. In Figure 5 the substrate profiles along a vertical line in the whole domain at different time steps are reported. Because of
(a) Phenol concentration data (Bekins et al., 1998) from methanogenic degradation of phenol in a batch laboratory microcosm (symbols) and simulation results (lines). (b) Field data (Bekins et al., 1998) from a creosote-contaminated site (symbols) and simulation results (lines).
1390
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
Figure 4.
Distribution of substrate concentration at the initial (left) and final step (right).
the presence of the reaction term, the substrate concentration inside the domain diminishes in time, and inside the structure domain reaches values closer to the boundary one. As a consequence, the substrate concentration inside the overall domain is not in equilibrium and produces a mass flux from the fluid to the structure. The presence of the flux through the interface, in turn, produces biofilm growth. In Figure 5 the biofilm growth is also visible, since the dashed lines, referring to the structure domain, are longer at increasing time. Because of the presence of substrate reaction and biofilm growth, a check of the conservation of substrate mass inside the domain is necessary. For this purpose, the change of mass inside the domain, the mass flux through the boundaries, the change of mass due to reaction as well as the overall sum are reported in the following Figure 6. In particular, Figure 6b
Figure 5. Simulated concentration profiles in time along the vertical axis in the fluid (solid lines) and in the structure domains (dashed lines) at different time steps.
refers to the course during the first inner timeloop, and it shows a change in time of the flux through the boundaries and of the mass inside the domain for the first steps. Afterwards the solution reaches a steady-state condition and no more change of mass inside the domain is present anymore. Clearly, the overall sum is always zero and the mass is conserved in the inner timeloop. On the other hand, Figure 6a presents the changes at the end of each inner timeloop before the following growth step. In this case, no change of mass inside the domain is present, since each point represents a steady state solution. As a result, the change in mass due to reaction equals the mass flux out of the domain and the conservation of mass holds true for the multi-scale growth model, too. Influence of Phenomenological Law on Biofilm Growth The final biofilm shape is obviously determined by the values of the constants introduced in Equation (19). As a matter of fact, the presence of all three different terms and their values deeply influences the biofilm growth, while being determined by the particular biofilm under consideration and by its history. For the purpose of highlighting the effect of each contribution to the final biofilm shape, different simulations were run, only taking one term into account at a time. For these simulations a flow channel (dimensions 6 mm 9 mm 0.2 mm) was considered, inside which a circular biofilm structure VS (radius 0.5 mm and thickness 0.2 mm) is positioned. A constant velocity and concentration is imposed at the left boundary, while no-slip conditions are imposed at the top and bottom boundaries, in order to not develop any additional boundary layer. The scalar distribution in the fluid domain at the end of the first inner timeloop is reported in Figure 7. In this condition, the mass flux in the front part of the structure is higher than in its back part. For this reason if the growth is calculated only taking in account the mass flux, i.e. only the constant K1, the biofilm grows on all the surface but more on the front of it, as reported in Figure 7b. Also
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1391
Figure 6. P R ID
Time-dependent progress at the end of each inner timeloop in overall mass within the fluid and structure domain (a) and in the first timeloop (b) of the change R P FI dVI =Dt, of the mass flux out of the domain I FI c DI rFI nI ; of the change of mass due to reaction RS dVS and of the sum of all the previous terms.
normal stresses are higher in the front of the biofilm, producing an erosion of the structure, as is clearly visible when only the term relative to erosion due to normal stresses is taken into account (Fig. 7c). With regard to shear stresses, it has to be taken in account that their distribution is not uniform on all the front part of the structure and in particular that they are zero exactly on the mid point of the front part. A growth law, which only allows for the contribution of shear stresses, will obviously calculate a zero growth in that position, while producing erosion on the sides. This effect will produce a non-smooth structure shape, as reported in Figure 7d. The equation of growth (19), however, will never take in account only the contribution due to shear stresses, since in the reality the effect due to mass flux will always be present. For this reason, the highlighted problem will usually not be present. Effect of Operating Conditions on Biofilm Growth The model is able to predict the expected behavior not only when the structure is exposed to uniform fluxes but also when different concentration gradients are present. Actually, the operating conditions are known to have a big influence on biofilm growth and this effect is highlighted in this example. For this purpose a fluid domain VF (dimensions 10 mm 8 mm 0.2 mm) bound by a biofilm structure VS (dimensions 10 mm 2 mm 0.2 mm) at the bottom is taken into account. Different boundary conditions were applied to the same domain, while using constant material parameters. The interface G between the biofilm and the fluid represents the interface for FSI and mass transfer, as well as for biofilm growth. Only the interface G is allowed to move in any direction, while all the other boundaries are kept fixed. When not differently specified a zero velocity field is also applied at the fluid boundaries and a zero-flux conditions at the solid and fluid ones. As an initial condition, a zero velocity field is also always prescribed.
1392
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
Three different case were considered. (i) A constant unitary concentration is applied at the top of the fluid domain. (ii) A cosine scalar profile with a unitary maximum value at its center and zero-values at the sides is applied at the top boundary of the fluid domain. (iii) Conditions similar to the case (ii) were used, except for the fact that also a non-zero velocity x-component was applied at the left fluid boundary, with a parabolic profile in the vertical direction, a null value on G and a maximum value of 0.1 mm/ms on the top. In this case also flow through the right fluid boundary is allowed and a slip boundary condition is prescribed at the top boundary. Also for this example a biological time step 1000 times larger than the fluid dynamic one was used. In Figure 8 the scalar distribution in the whole domain and the fluidstructure interfaces are reported for the three different cases and at different biological time steps. The influence of the different boundary conditions on the final biofilm shape is evident. Moreover, when the scalar distribution is uniform, then a uniform growth is present (i), for not-uniform scalar fields, a larger growth is present in the direction of the higher scalar concentration (ii), while the presence of the velocity field distorts the scalar distribution producing a non-symmetric final biofilm shape (iii). FSI and Growth of a Finger-Like Biofilm Structure In this example, the proposed approach is applied to the simulation of a biofilm finger under realistic operating conditions and material properties. For this purpose the geometry reported in Figure 9 is taken into account. The computational domain represents a part of a flow channel (dimensions 0.6 mm 0.3 mm 0.01 mm) inside which a
Figure 7. Scalar distribution around a circular biofilm structure (top) and initial (a) and final biofilm shape superimposed to the initial one for K1 6¼ 0 and K2 ¼ K3 ¼ 0 (b), for K2 6¼ 0 and K1 ¼ K3 ¼ 0 (c), for K3 6¼ 0 and K1 ¼ K2 ¼ 0 (d).
Figure 8.
Concentration distribution in the fluid and structure domains for the three different cases after 0, 20, and 40 biological time steps.
finger-like structure is present (0.1 mm high and 0.04 mm wide). The simulation is run under the realistic operating conditions and material parameters summarized in Table I for the cases of rigid and deformable biofilm. At the inlet boundary parabolic velocity and concentration profiles are imposed, while at the top boundary and at the bottom one a slip and a no-slip condition for the velocity are prescribed respectively. As initial conditions a zero velocity field and a parabolic concentration profile in the vertical direction are imposed. The simulation was run for an hydrodynamic timestep of 5 ms and a biological timestep of 1day. The first FSI and mass transport timeloop was run for 500 steps until a first steady state solution was obtained. In the following inner timeloops, the velocity and scalar distribution calculated in the last inner timeloop, before the growth step, were imposed as initial conditions. In this way, the number of time steps in the following inner timeloops to reach a steady state solution was reduced to 20. In Figure 9 the initial biofilm shape and the one after 8 biological timesteps, that is, after 8 days, is reported for the two different cases. The different final biofilm shape for a fixed structure (Fig. 1b) and a deformable one (Fig. 1c) is apparent, demonstrating that the growth can be affected by the deformation and highlighting the importance of the proposed approach. The boundary and the operating conditions chosen in the present example, together with the
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1393
prescribed growth constants, produce a larger growth of the finger in the width direction instead of in height and the final finger results thicker but not much higher than the initial one. For different boundary and operating conditions, material parameters and biofilm growth constants, different biofilm shapes can be obtained. For these reasons, the results reported in Figure 9 only show a possible final biofilm shape. The present example demonstrates, however, the applicability of the present approach to the simulation of realistic biofilm structure under realistic material and operating conditions and the importance of the proposed coupled approach.
Conclusions
Figure 9. Definition of the computational domain (top), initial (a) and final shape of the finger-like biofilm, in the case of fixed structure (b), deformable structure in undeformed (c), and deformed configuration (d).
Table I. Material parameters, operating conditions and growth constants. Parameter Liquid Dynamic viscositya Densitya Inlet velocity Biofilm Densityb Young’s modulusc Poisson ratiob Substrate Diffusion coefficientd Uptake rate coefficiente Saturation coefficientf Inlet concentrationg Growth Normal flux constant Normal stress constant Shear stress constant
Symbol
Value
Unit
mF rF uFIN
103 103 1.5 102
kg m1 s1 kg m3 m s2
rS ES nS
3 102 102 4 101
kg m4 kg m1 s2
D k K FFIN
2.5 109 3 102 3 103 2.5 102
m2 s1 mol m3 s1 mol m3 mol m3
K1 K2 K3
4 105 102 1
s m3/mol s m2/kg s m2/kg
Water at 20 C. TanyolaSc and Beyenal (1997). c Böl et al. (2013). d Dissolved oxygen in water at 20 C. e Assumed for oxygen uptake by heterotrophic microorganisms (Henze, 2000). f Affinity for oxygen of heterotrophic microorganisms (Henze, 2000). g Relatively low dissolved substrate concentration. a
b
1394
Biotechnology and Bioengineering, Vol. 111, No. 7, July, 2014
To enable a better understanding of biofilm growth and development, we developed a novel and advanced computational model for the numerical simulation of FSI and mass transfer of growing biofilm structures. For the purpose of describing correctly and in enough detail the different phenomena involved, the fully coupled nonlinear FSI problem and the nonlinear multi-field mass transport model are encapsulated in the biofilm growth model through a multi-scale approach in time. The described new methodology enables the mutual coupling of (i) flow velocity and solid deformation (ii), mass transport on deforming fluid and solid domains, (iii) as well as growth, based on local mass transport and stresses. The proposed model takes in account the effect of nutrients flux and stresses on the biofilm surface through a simple phenomenological equation. Although the effect of mass transfer on biofilm growth as well as the influence of structure deformation on mass transport have already been considered before, none of the previous approaches were able to take into account all these important phenomena at the same time and without any restrictions. Furthermore, considering FSI, mass transport and growth at the same time is essential for understanding the formation and development of biofilm structures like streamers. The general validity of the proposed model was illustrated by selected numerical examples. The accurate implementation of the non-linear Monod kinetic was highlighted through comparison with experimental data. The mass conservation was shown. The effect of phenomenological law constants on biofilm growth was highlighted. The influence of the different operating conditions was discussed. The applicability of the novel approach to a realistic structure and at realistic operating conditions and material parameters was reported. As a result, selected examples successfully demonstrated the general suitability of the proposed model for modeling the FSI and the mass transfer of growing biofilm structures. After an appropriate calibration of the phenomenological law driven by experimental data and evidence, the application of this coupled and multi-scale approach to different real biofilm structures will be explored. Future work will also concern the introduction of a detachment model and its
application to new realistic scenarios also from different fields of biomechanics. References Alpkvist E, Klapper I. 2007a. A multidimensional multispecies continuum model for heterogeneous biofilm development. Bull Math Biol 69: 765–789. Alpkvist E, Klapper I. 2007b. Description of mechanical response including detachment using a novel particle model of biofilm/flow interaction. Water Sci Technol 55:265–273. Battin TJ, Kaplan LA, Newbold JD, Cheng X, Hansen C. 2003. Effects of current velocity on the nascent architecture of stream microbial biofilms. Appl Environ Microbiol 69:5443–5452. de Beer D, Schramm A. 1999. Micro-environments and mass transfer phenomena in biofilms studied with microsensors. Water Sci Technol 39:173–178. Bekins BA, Warren E, Godsy EM. 1998. A comparison of zero-order, firstorder, and Monod biotransformation models. Ground Water 36: 261–268. Böl M, Möhle RB, Haesner M, Neu TR, Horn H, Krull R. 2009. 3D finite element model of biofilm detachment using real biofilm structures from CLSM data. Biotechnol Bioeng 103:177–186. Böl M, Ehret AE, Bolea Albero A, Hellriegel J, Krull R. 2013. Recent advances in mechanical characterisation of biofilm and their significance for material modelling. Crit Rev Biotechnol 33:145–171. Cogan NG. 2008. Two-fluid model of biofilm disinfection. Bull Math Biol 70:800–819. Cogan NG. 2011. Computational exploration of disinfection of bacterial biofilms in partially blocked channels. Int J Numer Method Biomed Eng 27:1982–1995. Coroneo M, Yoshihara L, Bauer G, Wall WA. 2012. A finite element approach for the coupled numerical simulation of fluid-structure interaction and mass transfer of moving biofilm structures In: Eberhardsteiner J, et al. editor European Congress on Comp Meth Appl Sci Eng, ECCOMAS. Dillon R, Fauci L, Fogelson A. 1996. Modeling biofilm processes using the immersed boundary method. J Comput Phys 129:57–73. Dockery J, Klapper I. 2001. Finger formation in biolm layers. SIAM J Appl Math 62:853–869. Donea J, Huerta A, Ponthot JP, Rodriguez-Ferran A. 2004. Arbitrary Lagrangian–Eulerian methods. In: Stein E, de Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics, Vol. 1: Fundamentals. Chichester, West Sussex: John Wiley Sons. Duddu R, Bordas S, Chopp D, Moran B. 2008. A combined extended finite element and level set method for biofilm growth. Int J Numer Method Eng 74:848–870. Duddu R, Chopp DL, Moran B. 2009. A two-dimensional continuum model of biofilm growth incorporating fluid flow and shear stress based detachment. Biotechnol Bioeng 103:92–104. Eberl HJ, Sudarsan R. 2008. Exposure of biofilms to slow flow fields: The convective contribution to growth and disinfection. J Theor Biol 253:788–807. Figueroa CA, Baek S, Taylor CA, Humphrey JD. 2009. A computational framework for fluid-solid-growth modeling in cardiovascular simulations. Comput Methods Appl Mech Eng 198:3583–3602. Garny K, Horn H, Neu T. 2008. Interaction between biofilm development, structure and detachment in rotating annular reactors. Bioprocess Biosyst Eng 31:619–629. Gee MW, Küttler U, Wall WA. 2010. Truly monolithic algebraic multigrid for fluid-structure interaction. Int J Numer Method Eng 85:987–1016.
Graf von der Schulenburg DA, Pintelon TR, Picioreanu C, Van Loosdrecht MCM, Johns ML. 2009. Three-dimensional simulations of biofilm growth in porous media. AIChE J 55:494–504. Hall-Stoodley L, Stoodley P. 2002. Developmental regulation of microbial biofilms. Curr Opin Biotechnol 13:228–233. Henze M. 2000. Activated sludge models ASM1, ASM2, ASM2d, and ASM3. London, UK: IWA Publishing. p 128. Holzapfel G. 2000. Nonlinear solid mechanics. A continuum approach for engineering. Chichester, UK: Wiley. p 455. Klapper I. 2004. Effect of heterogeneous structure in mechanically unstressed biofilms on overall growth. Bull Math Biol 66:809–824. Klöppel T, Wall WA. 2010. A novel two-layer, coupled finite element approach for the nonlinear elastic and viscoelastic behavior of human erythrocytes. Biomech Model Mechanobiol 10:445–459. Küttler U, Gee M, Förster Ch, Comerford A, Wall WA. 2010. Coupling strategies for biomedical fluid-structure interaction problems. Int J Numer Method Biomed Eng 26:305–321. Pereira MO, Kuehn M, Wuertz S, Neu T, Melo LF. 2002. Effect of flow regime on the architecture of a Pseudomonas fluorescens biofilm. Biotechnol Bioeng 78:164–171. Picioreanu C, van Loosdrecht MCM, Heijnen J. 1999. A discrete-differential modelling of biofilm structure. Water Sci Technol 39:115–122. Radu AI, Vrouwenvelder JS, van Loosdrecht MCM, Picioreanu C. 2012. Effect of flow velocity, substrate concentration and hydraulic cleaning on biofouling of reverse osmosis feed channels. Chem Eng J 188:30–39. Simões M, Pereira M, Sillankorva S, Azeredo J, Vieira MJ. 2007. The effect of hydrodynamic conditions on the phenotype of Pseudomonas fluorescens biofilms. Biofouling 23:249–258. Smith BG, Vaughan JR BL, Chopp DL. 2007. The extended finite element method for boundary layer problems in biofilm growth. Commun Appl Math Comput Sci 2:35–56. Stoodley P, Lewandowski Z, Boyle J, Lappin-Scott HM. 1998. Oscillation characteristics of biofilm streamers in turbulent flowing water as related to drag and pressure drop. Biotechnol Bioeng 57:536–544. Stoodley P, Dodds I, Boyle J, Lappin-Scott H. 1999. Influence of hydrodynamics and nutrients on biofilm structure. J Appl Microbiol 85:19S–28S. Stoodley P, Sauer K, Davies DG, Costerton1 JW. 2002. Biofilms as complex differentiated communities. Annu Rev Microbiol 56:187–209. Taherzadeh D, Picioreanu C, Küttler U, Simone A, Wall WA, Horn H. 2010. Computational study of the drag and oscillatory movement of biofilm streamers in fast flows. Biotechnol Bioeng 105:600–610. Taherzadeh D, Picioreanu C, Horn H. 2012. Mass transfer enhancement in moving biofilm structures. Biophys J 102:1483–1492. TanyolaSc A, Beyenal H. 1997. Prediction of average biofilm density and performance of a spherical bioparticle under substrate inhibition. Biotechnol Bioeng 56:319–329. Teodósio JS, Simões M, Alves MA, Melo LF, Mergulh~ao FJ. 2012. Setup and validation of flow cell systems for biofouling simulation in industrial settings. Sci World J 12. Towler BW, Cunningham A, Stoodley P, McKittrick L. 2007. A model of fluidbiofilm interaction using a Burger material law. Biotechnol Bioeng 96:259–271. Wall WA, Ramm E. 1998. Fluid-structure interaction based upon a stabilized (ALE) finite element method In: Idelsohn SR, Onate E, Dvorkin EN, editors. Computational mechanics—New trends and applications. Barcelona, Spain: CIMNE. pp 1–20. Wall WA, Gee MW. 2013. BACI—A multiphysics simulation environment. Munich, Germany: Technische Universität München. Wijeyekoon S, Mino T, Satoh H, Matsuo T. 2004. Effects of substrate loading rate on biofilm structure. Water Res 38:2479–2488.
Coroneo et al.: Biofilm Growth: A Coupled Multi-Scale Approach Biotechnology and Bioengineering
1395
