Simulating biofilm deformation and detachment with the immersed boundary method

2.4 MB·English

by  Rangarajan Sudarsan

Simulating biofilm deformation and detachment with the
immersed boundary method
5
1 Rangarajan Sudarsan1,∗, Sudeshna Ghosh2, John M. Stockie2 and Her-
0
2 mannJ. Eberl1

n 1 Department of Mathematics and Statistics and Biophysics Interdepartmental Graduate Pro-
a
gram,UniversityofGuelph,ON,N1G2W1,Canada.
J
2DepartmentofMathematics,SimonFraserUniversity,Burnaby,BC,V5A1S6,Canada.
8
2

]
n
Abstract. Weapplytheimmersedboundary(orIB)methodtosimulatedeformationanddetach-
y
mentofaperiodicarrayofwall-boundedbiofilmcoloniesinresponsetoalinearshearflow. The
d
- biofilmmaterialisrepresentedasanetworkofHookeanspringsthatareplacedalongtheedges
u of a triangulation of the biofilm region. The interfacial shear stress, lift and drag forces acting
l
f onthebiofilm colonyarecomputedbyusing fluidstressjump method developedbyWilliams,
.
s FauciandGaver[Disc.Contin.Dyn.Sys.B11(2):519–540,2009],withamodifiedversionoftheir
c
exclusionfilter.Ourdetachmentcriterionisbasedonthenovelconceptofanaveragedequivalent
i
s continuumstresstensordefinedateachIBpointinthebiofilmwhichisthenusedtodetermine
y
acorrespondingvonMisesyieldstress;whereverthisyieldstressexceedsagivencriticalthresh-
h
p oldtheconnectionstothatnodearesevered,therebysignallingtheonsetofadetachmentevent.
[ In order to capture the deformation and detachment behaviour of a biofilm colony at different

stages of growth, we consider a family of four biofilm shapes with varying aspect ratio. Our
1
numericalsimulationsfocusonthebehaviourofweakbiofilms(withrelativelylowyieldstress
v
1 threshold) and investigate features of the fluid-structure interaction such as locations of maxi-
2 mum shear and increased drag. The most important conclusion of this work is that the com-
2 monly employed detachment strategy in biofilm models based only on interfacial shear stress
7
canleadtoincorrectorinaccurateresultswhenappliedtothestudyofshearinduceddetachment
0
of weakbiofilms. Our detachmentstrategybasedonequivalent continuum stressesprovidesa
.
1
unifiedandconsistentIBframeworkthathandlesbothsloughinganderosionmodesofbiofilm
0
detachment,andisconsistentwithstrategiesemployedinmanyothercontinuumbasedbiofilm
5
1 models.
: AMSsubjectclassifications:74D10,74F10,76D05,76M20.
v
i Key words: immersed boundary method, biofilms, detachment, von Mises yield stress, interfacial shear
X
stress,dragandliftforce.
r
a
1 Introduction
Thesubjectofthisworkistheflow-induceddeformationofabiofilmcolony,whichisamesoscale
collection of bacterial cells held together by an extracellular polymeric network (EPS) that is se-
cretedbythecells. Thedimensionsofabiofilmcolonycanbeanywherefromtenstohundredsof
microns,whereasthesizeofanindividualbacterialcellmakingupthecolonyisontheorderof1–5
∗Corresponding author. Email addresses: [email protected] (R. Sudarsan), [email protected] (S. Ghosh),
[email protected](J.M.Stockie),[email protected](H.J.Eberl)
http://www.global-sci.com/ GlobalSciencePreprint
2
microns;ourfocusisoncontinuummodelsthattreatthebiofilmasaviscoelasticsolidcontinuum
rather than incorporating the dynamics of individual bacteria. The flow-induced deformations
of the biofilm colony affect the fluid dynamic forces acting on it, and thereby also alter both the
extent and the mode of detachment (i.e., sloughing or erosion) that may be experienced by the
biofilm. Weareparticularlyinterestedinunderstandingwhetherbiofilmcoloniesgainanyprotec-
tionagainstdetachmentwhentheyareincloseproximitytoothercolonies. Tothisend,ouraimis
todeveloparobustnumericalmethodforsimulatingtheinteractionbetweenabiofilmcolonyand
thesurroundingfluidthatiscapableofcapturingthedifferentmodesofbiofilmdetachment.
Our approach is based on the immersed boundary (or IB) method, in which the biofilm con-
tinuum is replaced by a network of Hookean springs. Although the IB method has already been
usedbyseveralauthorsforstudyingbiofilmdeformationanddetachment[1,26],ourapproachof
computing an equivalent continuum stress at each IB node is markedly different and using it to
initiatedetachmentprovidesuswithawayofhandlingdetachmentinamannerthatisconsistent
withothercontinuummechanicsbasedmodelssuchas[13].
1.1 Bacterialbiofilms
Bacterialbiofilmsareaggregationsofmicrobesthatgrowonsurfacesinanaqueousenvironment.
They form when bacterial cells suspended in the fluid attach themselves to a surface and begin
producingan extracellular polymericsubstancein which thegrowing bacteria cells embed them-
selves. Biofilmsplaycontrastingrolesdependingonthescenariosinwhichtheyareencountered:
in waste-water treatment [39] and environmental engineering they play a helpful role; whereas
biofilms encounteredon medical devicesorfood-processingequipmentcan bedetrimentalin the
sensethattheycauseinfections[9]orcompromisefoodsafety[48].
In theseand other applications, biofilms experience vastly different physical conditions (tem-
perature and pH), hydrodynamics (ranging from creeping to turbulent flow) and chemical envi-
ronments(rich or sparsein nutrients,or saturatedwith antibiotics). Their adaptability to diverse
conditionsis believed to derive at least in part from the mechanical and chemical protection pro-
videdbythegel-likeEPSlayer. Itis importantfroman engineeringstandpointtounderstandthe
mechanical propertiesofbiofilms sothat wecan deviseeffective methodsfor notonlyenhancing
biofilm growth and survival but also removing them from surfaces. Consequently, over the past
10yearssignificantefforthasbeenexpendedtodevelopnovelrheologicalmethodsformeasuring
mechanical properties of biofilms [25]. Such experiments have typically arrived at different con-
clusionsonhowtocharacterizebiofilms, withsomeconcludingthattheybehaveaselastic[40]or
viscoelastic solids [28], while others liken biofilms to viscoelastic fluids [55]. Furthermore, even
within the same material class, measured material parameter values can vary over a fairly wide
range. The general consensus is that biofilms behave as (visco-)elastic solids when the applied
fluidshear stressis at ornear thestressesat which thebiofilm wasgrown,while at highervalues
of shear stress the biofilm can yield and behave as a viscoelastic fluid. In this paper, we restrict
ourselves to the low shear stress case and model the biofilm mechanically as a viscoelastic solid
embeddedwithinaviscousfluid.
1.2 Mathematicalmodelsof biofilm growth, deformation anddetachment
Mathematicalmodelsforsimulatingbiofilmgrowthmusttakeintoaccountawiderangeofmech-
anisticandotherdynamicalprocesses,includinggrowthanddeathofbacteria,attachmentofcells
to the substratum, transportof solutes(nutrients, metabolic products, antibiotics) within the sur-
3
rounding fluid and the biofilm itself, reaction kinetics, and removal of bacterial cells as clusters
(sloughing) or as individual cells (erosion) when the biofilm EPS matrix weakens in response to
fluid shear or chemical treatments. The earliest models developed in the 1980’s [29,61] assumed
that the biofilm is one-dimensional, while many subsequent modelling efforts have attemptedto
capturethe2Dor3Dmorphologythatdevelopsduringbiofilmgrowthprocesses.Thetreatmentof
thehydrodynamicsanditsinteractionwiththebiofilmhasreceivedvaryingdegreesoftreatment
in thesemodels, ranging from some modelsthat consider biofilm transportprocessesin isolation
anduseaspecifiednutrientconcentrationboundarylayer thicknessat thebiofilm-fluid interface,
whereasothermodelsperformfull fluidflow andnutrienttransportcalculations withor without
incorporating fluid-structure interaction effects. A variety of approaches have been developedto
modelthegrowthandspreadingofthebiomass, including individual-based models[32], cellular
automata [64], continuum models [2,13,30] and phase field models [67]. An exhaustive review
ofthedifferentmodelingapproachescanbefoundinthearticle[60]andperformancebenchmark
comparisonsareavailable in[14].
Amongthemorecompletemodelsarethosethat capturemulti-dimensionalgrowthandfluid
flow [13,15,16,44], although the role played by hydrodynamics in inducing biofilm deformation
anddetachmenthasmostoftenreceivedonlyadhocorapproximatetreatment. Theseapproxima-
tionsareaimedatcapturingthesolidmechanicsgoverningthedynamicsofthedeformingbiofilm
while simplifying as much as possiblethe complex fluid flow that surroundsit. Thefirst attempt
at solving the solid mechanics problem insidethe biofilm and usingit toinitiate detachment was
in [44] where the biofilm was treated as a linearly elastic material and a von Mises yield stress
criterion was used to initiate detachment; however, this work neglected the effects of biofilm de-
formation. In contrast, the particle-based biofilm model in [65] ignored the fluid and proposed
a method for initiating detachment using a detachment speed function that is based on the normal
velocityofthebiofilm-fluidinterface. Asimilarapproachhasbeenusedin[49],whichincludesthe
effectsofbothbiofilmgrowthandflowbyemployingadetachmentspeedthatdependsonthein-
terfacialshearstress. Thisapproachtoinitiating detachmentrestsontheassumptionthatregions
wheretheinterfacialshearstressishighestcorrespondtolocationswherethebiofilmstrainisalso
high.
Oneoftheaimsofthecurrentstudyistoverifythevalidityofthislastassumptionbyconsider-
ingafullfluid-structureinteractionsimulationthatiscapableofdeterminingstressesinbothfluid
and biofilm. We do not explicitly model biofilm growth but rather mimic the effects of growth
by considering a family of biofilm colony shapesof different aspect ratio, where each colony size
corresponds to a different instant of time during the growth of the same colony. In each case,
we investigate how the deformation affects both the mode and the extent of biofilm detachment.
Thisisincontrastwithotherapproaches[13,44]wherethesolidmechanicsaresimulatedbutany
deformationsofthebiofilmareneglected.
1.3 Fluid-structure interaction inbiofilms
During the last several years, significant progress has been made in the study of fluid-structure
interaction (FSI) in biofilms, driven by the increased availability of experimental data on biofilm
mechanical properties and the motivation to understand the role they play in biofilm survival.
MostFSIstudiesneglectbiofilmgrowthbytakingadvantageofanaturalseparationoftimescales,
in that growthprocessesare veryslow in relation tofluid motionand biofilm deformation. More
recently, a phase field method has been applied successfully in 2D [34] and 3D [47] to simulate
biofilm growth coupled with deformation arising from interaction with the surrounding flowing
4
fluidtreatingthebiofilmcontinuumasamultiphasepolymericgel. Withtheexceptionofthesetwo
works,themostcommonapproachinthebiofilmFSIliteraturecombinesaLagrangiandiscretiza-
tionofthebiofilmwith an Arbitrary Lagrangian Eulerian (orALE)formulation forthefluid. The
firststudyofthiskindappearedin[54]wheretheystudiedthedeformationofa2Dhemispherical
biofilmcolonyplacedinaturbulentflowusingtheANSYS® commercial softwarepackage. More
recently [53], a similar ALE approach was used to study flow-induced oscillations of 2D biofilm
streamers and their effect on mass transfer. A more realistic 3D biofilm model was studiedin [7]
that used sliced 3D confocal laser scanning microscopy data to construct the colony shapes, and
employed a nonlinear hyper-elastic constitutive model that accommodates detachment based on
avonMisesyieldstresscriterion.
Theaforementionedapproacheshavetheadvantageofbeingwell-establishedintheliterature
and capable of easily incorporating experimental parameters that measure biofilm rheology. The
primarydisadvantageistheirhighcomputationalcostandalgorithmiccomplexitythatresultfrom
needingtoconstantlyre-meshthefluiddomain as thebiofilm colonydeforms. Moreover,thisre-
meshingcostincreasesenormouslyifdetachmentisincorporatedinthemodel.
Motivatedbythedesiretodevelopasimplerandmoreefficientcomputationalapproach,Alp-
kvistandKlapper[1]proposedanalternateFSIstrategybasedontheimmersedboundary(orIB)
methodin which the biofilm is discretized at a set of moving Lagrangian points. TheseIB points
move relative toan underlyingfixedCartesian gridon which thefluid equationsare solved. The
elasticpropertiesofthebiofilmarecapturedbydistributingforcesontothefluidthatderivefrom
anetworkofHookeanspringsjoining theIBpoints. An incompressiblefluid pervadesbothfluid
and biofilm regions and the biofilm inherits the density and viscosity of the surrounding fluid,
so that the biofilm is actually a visco-elastic composite material consisting of elastic spring forces
and fluid viscous forces. This approach has the clear advantage that no re-meshing of the fluid
gridisrequired. Toillustratetheversatilityoftheirapproach,AlpkvistandKlappersimulatedthe
deformation ofboth 2D and 3D biofilm structures,while also incorporating a simple detachment
criterionbasedoncuttingindividualspringswhentheyarestretchedbeyondacriticallength. We
remark that other IB models for biofilms were developed prior to [1], namely the work of Dillon
andcollaborators[10,11];however,theseauthorswereconcernedwithslowflowandthedynam-
ics of individual bacterial cells aggregating and settling on the substratum, and hence the results
arerelevant todifferentphenomenaoccurringonmuch smaller spatialscalesandat muchearlier
stagesofbiofilmformation.
A number ofotherIB approaches have since appeared, such as [57] whoperformed3D simu-
lationsofbiofilmdeformation,comparingHookean(elastic,spring-only)andKelvin-Voigt(visco-
elastic, spring plus dashpot)models for the biofilm material. Using a parametric study of spring
stiffnessanddampingcoefficientandcomparisonswithexperiments,theyestablishedthatrealis-
ticbiofilmdeformationbehaviourcanbeobtainedusingtheIBmethod. Morerecently,Hammond
etal. [26]developedmoredetailed 2D and3D IBmodels forfragmentation ofabiofilm colonyin
which the location of actual bacterial cells from 3D images was used to determine coordinates of
IB points. They employed a similar spring network and detachment strategy as in [1], but they
allowedthedensityofthebiofilmtodifferfromthatofthefluid,andinsubsequentwork[27]also
extendedtheirapproachtohandlevariable viscosity.
DespitetheincreasingpopularityoftheIBmethodinbiofilmFSIstudies,twomainchallenges
remain to be addressed before the method is capable of simulating realistic biofilm deformation
and detachment. The first relates to connecting values of the IB spring parameters (elastic stiff-
nessanddampingcoefficients)toactualbiofilmmaterialproperties. Despitetheparametricstudy
in[57]thatshoweditwaspossibletodeterminesuitableparametersforgivenbiofilmcolonyshape
5
andIBspringnetworktopology,thereremainsasyetnoapriorimethodfordeterminingIBparam-
etersforagivenbiofilm.
The second challenge relates to initiating detachment in the IB framework in a way that is
consistentwithothermoreestablishedcontinuum-mechanics-basedbiofilm studiesmentionedin
Section 1.2. Whereas IB methods have so far used strain in any given spring as a measure for
initiating detachment, this is not a true measure of strain in the continuum mechanics context.
Indeed,Hammondetal.[26]demonstratedthatwhenusingspringstrainasadetachmentcriterion
the resulting biofilm shape is sensitive to the critical strain parameter, and so it is unclear how to
choosethisparametertomatchagivensetofbiofilmmechanicalproperties.
1.4 Objectivesand outline
In this paper, we aim to address the challenges identified at the end of the preceding section by
developing an immersed boundary approach for simulating biofilms that is capable of capturing
realisticdeformationanddetachmentbehaviours. Webeginwitha2DIBmodelinspiredbythatof
AlpkvistandKlapper[1],andextendthisworkguidedbytwomainobjectives. Ourfirstobjective
is to develop a novel approach for initiating biofilm detachment that is consistent with methods
employedin biofilm modelsusingcontinuummechanics basedcalculations toenact detachment.
Inthisway,wecanretainthesimplerspringnetworkrepresentationofthebiofilmcontinuumand
the advantages it offers, while also implementing a more physically realistic criterion for detach-
ment. Oursecondmajorobjectiveistoinvestigatetheeffectofflow-induceddeformationonboth
themechanical stability and mode ofdetachment experiencedby weak biofilm structureshaving
realistic shapes that resemble those grown under mass transfer limited conditions. This will al-
lowustobetterunderstandsuchfundamentalquestionsashowflow-induceddeformationaffects
theforces acting on biofilms, and how spatial clusteringof biofilm colonies can shield themfrom
detachmentbyreducingthehydrodynamicshearforces.
Our modellingapproach incorporatestheworkdonein severalprevious computationalstud-
ies of two-dimensional biofilms in [52,66] wherein we investigated the fluid shear-induced de-
tachment forces (drag and lift) acting on rigid biofilm colonies that are both uniformly and non-
uniformly spaced. In contrast to these earlier studies that were restricted to values of shear rate
exceeding 10 cm/s, we focus in this paper on more flexible weak biofilm colonies immersed in a
slowershearflowhavingshearratelessthan1 cm/s.
Theorganization oftheremainderofthispaperisasfollows. InSection2,wedefinetheprob-
lemgeometryand parameters,and describethegoverningequationsand correspondingnumeri-
calschemeforourbasicIBmodelframework. Section3containsthenovelbiofilm-relatedaspects
ofourIBmodelwherewederiveourapproachforcalculating fluidshearstressalongthebiofilm-
colonyfluidinterfaceandfromthatthedrag/liftforcesactingonthebiofilmcolony. InSection3.1
we describe how we implement a modified version of the exclusion filter devised by Williams et
al.in[63]foraccuratelyapproximatinginterfacial shearstressesoncurvedimmersedboundaries,
andthenin Section3.2wedescribehowweadopttheconceptofacontinuumstressaroundeach
IB node in the biofilm material. Section 3.3 discusses how these quantities are used to obtain a
more realistic detachment criterion for biofilm colonies in the IB context. Finally, in Section 4 we
performa number of numerical teststhat validate our numerical approach and also demonstrate
theadvantageofourIBmodelforsimulatingrealisticdeformationanddetachmenteventsinfluid
shear-inducedbiofilmdynamics.
6
2 Immersed boundary (IB) method
The IB method is both a mathematical formulation and a numerical method for simulating the
complex interaction betweena deformable solid material and a surroundingfluid. The approach
dates back to work of Peskin [42] who originally developedthe approach to simulate blood flow
interacting with heart muscle, and more recent theoretical and computational developments are
summarizednicely in thereview paper[43]. TheIBapproach hasbeenappliedextensivelytothe
study FSI in bio-fluid mechanics, including such problems as amoeboid locomotion [8], platelet
aggregation[19],andspermmotility[18].
The IB methodis a mixed Eulerian-Lagrangian approach wherein the fluid equations are dis-
cretizedonafixed,rectangular(Eulerian)meshwhereastheelasticstructureisdefinedbyasetof
(Lagrangian) IB points that move relative to the underlying fluid mesh as the structure deforms.
The effect of the immersed boundary on the fluid is represented using a singular source term in
the fluid momentum equations, which is distributed onto the underlying fluid grid by writing it
astheconvolutionofanelasticforcedensitywitharegularizeddeltafunction.
2.1 Problemgeometry
A diagram of the problem domain is given in Fig. 1, which depicts shape of a biofilm colony en-
countered during different stages of growth. The domain is delimited by two horizontal walls
separated by a distance H, where the bottom wall (serving as the biofilm substrate) is held sta-
tionary while the top wall moves horizontally with constant velocity U . The domain is filled
wall
with a viscous, incompressible fluid so that in the absence of biofilm, the top wall will generate
aflow that overtime approaches asteady-statecorrespondingtolinear (planar) shear with shear
rateG=U /H.
wall
Biofilm colonies having identical shape and size are placed along the bottomwall with a uni-
form spacing of D between adjacent colonies. In order to simulate biofilms at different growth
b
stagesoccuringundermasstransferlimitedconditions,weselectfourrepresentativebiofilmshapes
with increasing aspect ratio that have fixed width but increasing height. Motivated by shapes at
earlystagesofgrowthobservedinexperiments[12,31]and3Dsimulations[13,15,30,67],wechoose
anidealizedbiofilmshapedescribedbytheequation
x n1 y n2
+ =1 withn ,n >0andy>0, (2.1)
1 2
W /2 H
(cid:18) b (cid:19) (cid:18) b(cid:19)
whichishalfofasuper-ellipse withwidthW andheight H . Inthecasewhenn =n =2,Eq.(2.1)
b b 1 2
reduces to half of a regular ellipse, while further constraining W =2H yields a semi-circle with
b b
radius H .
b
We believe that this choice is a reasonable parameterization of the typical finger-like shapes
that dominate earlier stages of biofilm growth, in contrast with later stages that exhibit a classi-
cal mushroom-shaped colony having an ellipsoidal head on top of a thin cylindrical stem. For the
purpose of this study, we take values of n =n =2.5 and W =50µm, and select a sequence of
1 2 b
three initial shapes with H =25, 50 and 75 µm respectively. In addition, we consider two other
b
shapes: a semi-circle with diameter 40 µm (that is, n =n =2, W =40µm and H =20µm); and
1 2 b b
anirregularmushroom-likestructurewithW ≈60µmand H ≈280µm,whoseshapeisextracted
b b
fromfiguresin [1]. Thelattermushroomshapeis thetallest biofilm colonydepictedin Fig.1 and
ischosentoillustratetheeffectivenessofthedetachmentcriteriadevelopedinthisstudybydirect
comparisonswiththeresultsfrom[1,26].
7
Figure1: Theproblemdomainconsistingofastationarybottomwall(e-f)separatedbyadistance Hfrom atopwall(g-h)
that moves at velocity (U ,0). Biofilm colonies of width W and increasing height represent shapes at different stages
wall b
of growth. The walls are treated with IB points and the whole flow region is embedded in a slightly larger computational
domain (abcd) that is periodic in both x and y.
Rather than imposing wall boundary conditions directly on the fluid, we simplify the fluid
solver by treating the top and bottom walls using immersed boundaries and embed the problem
domain inside a slightly larger fluid domain denoted by dashed lines in Fig. 1. Because we are
interestedinstudyingrepeatingarraysofbiofilmscolonies,periodicboundaryconditionsareim-
posed in both the x and y directions so that the computational domain is Ω=[0,W]×[0,H] with
W=W +D (labeledabcdinthefigure),anditcontainsasinglebiofilmcolony.
b b
2.2 Governingequations
Thedomain Ω isfilledwithaNewtonianincompressiblefluidhaving constantdensityρ (g/cm3)
anddynamicviscosityµ(g/cms). Denotetheinfinitesimallythinimmersedboundariesrepresent-
ingthefixedbottomwallandmovingtopwall byΓ and Γ respectively,andlet Γ represent
bot top bio
thesolidelasticstructurecorrespondingtothebiofilmcolony. NotethatΓ isactuallyacomposite
bio
materialthatconsistsoftheelasticforce-generatingmaterialandfluidthatco-existwithinthesame
region. Wealsoassumeforthepurposesofthisstudythatthebiofilmisneutrallybuoyantandhas
thesamedensityasthefluid(althoughitisstraightforwardtoextendtheIBapproachtodealwith
variable density problems [22,26]). The fluid is therefore governed at all points x=(x,y)∈Ω by
theincompressibleNavier-Stokesequations
∂u
ρ +ρu·∇u=−∇p+µ∇2u+f, (2.2)
∂t
∇·u=0, (2.3)
8
whereu(x,t) (cm/s)isthefluidvelocityand p(x,t) ispressure(g/cms2).
Theeffectofthesolidboundariesonthefluidisencompassedinthefluidforcingtermf,which
we consider next. Assume that the configuration of the solid material making up both channel
walls and biofilm is described by a function X(q,t), where q is a generalized (dimensionless)pa-
rameterization that is either a scalar (q=s) in the case of the walls Γ and Γ , or else a vector
bot top
(q=(r,s)) for a solid region like Γ . The IB force in both cases is specified in terms of a discrete
bio
network of IB points connected by springs, and more detail on the precise form of these spring-
force connections for walls and biofilm is provided later in Sections 2.4 and 2.5. Assumethat the
force generated by any deformed configuration can be described by an IB force density function
F(X(q,t) dependingon the current stretchedconfiguration of the spring network. Then the fluid
force f may bedeterminedby spreadingtheforce densityatIBpointsontothefluidusingadelta
functionconvolution
f(x,t)= F(X,t) δ(x−X(q,t))dq, (2.4)
Γ
Z
whereδ(x)=δ(x)δ(y)istheCartesianproductoftwo1DDiracdeltafunctionsandΓ=Γ ∪Γ ∪
bio top
Γ representsthesetofallimmersedboundaries.
bot
The final equation required to close the system is an evolution equation for the immersed
boundaries,whichweassumemovewiththesamevelocityasthesurroundingfluid
∂X
= u(x,t)δ(x−X(q,t)dx. (2.5)
∂t Ω
Z
Thisissimplyanotherwayofstatingtheno-slipconditionforadeformable boundaryat location
X(q,t).
2.3 Numericalalgorithm
WenowdescribethebasicnumericalalgorithmforsolvingEqs.(2.2)–(2.5),whichisasemi-implicit
scheme very similar to the one employed in [50]. The fluid domain Ω is discretized on a regular
Cartesian mesh with coordinates x =(x ,y )=(ih ,jh ) for i=0,1,...,N and j=0,1,...,N , where
ij i j x y x y
h =W/N andh =H/N areconstantgridspacingsinthexandydirections(andweassumefor
x x y y
simplicitythath =h ). Thetimeintervalofinterest[0,T]islikewisedividedintoasequenceof N
x y t
equally-spacedpointsdenotedt =n∆tforn=0,1,...,N,where∆t=T/N isthetimestep. Letthe
n t t
discrete values of the velocity and pressure be denoted by un and pn respectively. Suppose that
ij ij
the immersed boundary is described by a set of N IB points whose locations at any time t are
b n
givenbyXn=(Xn,Yn)forℓ=1,2...,N . ThecorrespondingforcedensitiesaredenotedbyFn,with
ℓ ℓ ℓ b ℓ
theprecisespecificationoftheimmersedboundarydiscretizationandforcedensitycalculation Fn
ℓ
beinggiveninthethefollowingtwosections.
Wenowdescribeouralgorithmforupdatingfluidgridquantitiesun and pn andtheIBconfig-
ij ij
urationXn fromtimet totimet . Thealgorithmproceedsinfourmainsteps:
ℓ n n+1
Step1: Compute the force density Fn based on the current IB configuration Xn as described in
ℓ ℓ
Section2.5.
Step2: Spreadtheforcedensityontofluidgridpointsusingadiscreterepresentationofthedelta-
functionconvolutioninEq.(2.4)
N
b
fnij=∑Fnℓδh(xij−Xnℓ)A, (2.6)
ℓ=1
9
whereδ (x)isaregularizeddeltafunctiongivenby
h
1 x y
δ (x)= φ φ , (2.7)
h
h h h h
x y (cid:18) x(cid:19) (cid:18) y(cid:19)
with
1 1+cos(πr) , if|r|62,
φ(r)= 4 2 (2.8)
(0,(cid:0) (cid:1) otherwise.
The scaling factor A has units of length for forces generated by the 1D wall interfaces (for
which (2.6) approximates a line integral) whereas A has units of area for the 2D biofilm
region. To ensure that Eq. (2.6) is a consistentrepresentationof the correspondingintegrals
undergridrefinement,AisinverselyproportionaltothenumberofIBpoints. Detailsonthe
preciseexpressionusedforAineachcaseareprovidedinSection2.5.
Step3: IntegratetheNavier-StokesequationsusingChorin’ssplit-stepprojectionscheme:
(1)
a. Compute an intermediate velocity u by updating the velocity only for the contribu-
ij
tionfromtheIBelasticforce:
u(1)−un
ρ ij ij =fn. (2.9)
 ∆t  ij
 
(2) (3)
b. Compute intermediate velocities u and u by applying convection and diffusion
ij ij
termsusinganalternatingdirectionimplicit(ADI)approach:
u(2)−u(1)
ρ ij ij +unD0u(2) =µD+D−u(2), (2.10)
 ∆t ij x ij  x x ij
 (3) (2) 
u −u
ρ ij ij +unD0u(3) =µD+D−u(3). (2.11)
 ∆t ij y ij  x x ij
 
TheoperatorsD+andD−refertothestandardfirst-orderforwardandbackwarddiffer-
x x
enceapproximationsofthe x-derivative, and D0 isthestandardsecond-ordercentered
x
differenceapproximation. Analogousdefinitionsapplyforthey-derivativeapproxima-
tions D+, D− and D0. Equations (2.10) and (2.11) thus represent periodic tridiagonal
y y y
linearsystemsfortheintermediatevelocities.
(3)
c. Project theintermediatevelocity u ontothe space ofdivergence-freevector fields by
ij
firstsolvingaPoissonequationforthepressure pn+1
ij
ρ
∇ ·∇ pn+1= ∇ u(3), (2.12)
h h ij ∆t h ij
where ∇ =(D0,D0) is a centered approximation of the gradient operator and the dis-
h x y
creteLaplacian ∇ ·∇ yieldsawidefinitedifferencestencilthatspansfourgridpoints
h h
in each direction as described in [50]. Because of the periodic boundary conditions,
10
the pressure Poisson equation is solved most efficiently using a fast Fourier transform
(FFT).Finally,thevelocityprojectioniscompletedviathecorrection
∆t
un+1=u(3)− ∇ p . (2.13)
ij ij ρ h ij
Step4: Evolvetheimmersedboundarytotimet using
n+1
Xℓn+1=Xnℓ+∆t∑uinj+1δh(xij−Xnℓ)hxhy. (2.14)
i,j
TheIBalgorithmexplainedaboveisfirst-orderaccurateinbothtimeandspace. Despitetheuseof
second-order differences for spatial derivatives, the spatial accuracy reduces to first order owing
to the particular choice of interpolation used for spreading the fluid velocity onto the immersed
boundary[37].
This algorithm has the advantage that it is simple and easy to code, although it suffers from
afairly strict stability restrictionon the time stepowingto the explicit treatment oftheIB forcing
term in Step 1 of the algorithm. We implement the above algorithm using MATLAB®, and when
that is combined with the extra cost of implementing a realistic detachment criterion (see details
Section 3.3) our approach is restricted to fairly short-time simulations. The scope of this paper is
thereforelimited to the studyofbiofilm deformation and detachmentin the early stagesup until
aquasi-steadybiofilmconfigurationinreached. Anyfull-scaleimplementationofthedetachment
criteria described in this study would therefore benefit from a more efficient IB implementation
suchasthefullyimplicitapproachin[38]oroneoftheparallelapproachesdevelopedineither[24]
or[62].
2.4 Discrete representationof walls andbiofilm
In this section, we describe the discretization of the immersed boundaries (walls and biofilm) in
terms ofa networkof IB points connected by springs. Section 2.5 will then provide details of the
forcedensitycalculations basedonthisdiscreteIBconfiguration.
Webegin with thehorizontal walls at locations y=0 and H (labeled e-fand g-h in Fig. 1) that
are each replacedby a1Dperiodicarray ofIBpointsstretchingacross thedomain. Adjacentwall
pointsare not explicitly connectedto each other, but instead each IB point is connectedby a very
stiffspringto a correspondingtether or target point that initially occupies thesame positionas the
IB point. The tether forces generated by these springs constrain the IB points to remain close to
their target locations and hence mimic a solid wall. The initial IB point spacing h is chosen
wall
sothat h <1max(h ,h ), which helpstocontrolnumerical errorsthatwould otherwiselead to
wall 2 x y
significantleakageoffluidbetweenIBpoints[43].
Nextweconsiderthediscreterepresentationofthebiofilmcolony,whichisobtainedbytrian-
gulatingthebiofilmregion,placingIBpointsatthenodallocations,andconstructinganetworkof
Hookeanspringscorrespondingtotheedgesinthetriangulation. Thespringstiffnessischosento
approximatetheelasticpropertiesofanactualbiofilmmaterial. Thebiofilmisattachedtothelower
wall byconnecting pointsalong thebottomofthecolony toacorrespondingwall IBpointwith a
stiff spring, where each pair of points initially occupies the same location. To ensure that such a
spring network mimics a mechanically isotropic biofilm, we use an initial triangulation that has
approximatelyuniformshapeand size oftheelements. Forthis purpose,we usetheopen-source
MATLAB package DistMesh of Persson and Strang [41]. The average spacing between biofilm IB