Biodigester Rapid Analysis and Design System (B-RADeS): A candidate attainable region-based simulator for the synthesis of biogas reactor structures

ContentslistsavailableatScienceDirect

Computers and Chemical Engineering

journalhomepage:www.elsevier.com/locate/compchemeng

Biodigester rapid analysis and design system (B-RADeS): A candidate attainable region-based simulator for the synthesis of biogas reactor structures

F. Abunde Neba

^a,d,e,f,∗

, Nana Y. Asiedu

^c

, Ahmad Addo

^b

, John Morken

^g

, Stein W. Østerhus

^d

, Razak Seidu

^e

aAbunde Sustainable Engineering Group (AbundeSEG), Kumasi, Ghana

bDepartment of Agricultural and Biosystems Engineering, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana

cDepartment of Chemical Engineering, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana

dDepartment of Civil and Environmental Engineering, Norwegian University of Science and Technology, Trondheim, Norway

eInstitute for Marine Operations and Civil Engineering, Norwegian University of Science and Technology, Alesund, ˚ Norway

fDepartment of Mathematics, Kwame Nkrumah University of Science and Technology, Kumasi, Ghana

gFaculty of Science and Technology, Norwegian University of Life Sciences, Drobakveien 31, 1432 Aas, ˚As, Norway

a rt i c l e i n f o

Article history:

Received 17 April 2019 Revised 10 September 2019 Accepted 16 October 2019 Available online 19 October 2019 Keywords:

Attainable regions Kinetic models

Macroeconomic parameters Process economics Biodigester design

a b s t r a c t

Anaerobicdigestersareseldomdesignedbasedonprocesskinetics,butrather onacombinationofhy- draulicandorganicloading,whichmaylimitoperationalperformance.Thisstudyfocusesontheincor- porationofprocesskineticsinthedesignofanaerobicdigesters,withintheattainableregionconceptual framework.Candidateattainableregionsforanaerobicdigestersareidentifiedusingthesoftwareenviron- mentBiodigesterRapidAnalysisandDesignSystem(B-RADeS),whichcouples,biodegradationkineticsas wellaseconomicparametersforthesynthesisofbiodigesterstructures. Byconsideringswine,palmoil andpharmaceuticalwastewaters,paybackperiodsof0.5,1&2years,andsubstrate,kineticmodeland/or economicparameters,apromisingdigesterstructure(andassociated hydraulicretentiontimes)issyn- thesized,consistingofaCSTRfollowedbyPFR(15days),CSTR(4.8hours)andaPFRwithbypassoffeed (3days).Theframeworkoffersgreatpromiseforwidespreadpracticalapplication.

© 2019TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Anaerobic digestionof solid waste and/or wastewater sludges has longbeen usedforstabilization ofthesewastes prior todis- posal. Among thebenefits involved inanaerobicwaste treatment compared to the aerobic counterpart are: improved dewaterabil- ity of the treated waste and generation of renewable bioenergy (Mes et al., 2003). The construction and operation of anaerobic treatment plants requires optimizing the techno-economic feasi- bility by defining optimal process configuration of anaerobic di- gesters. There exist several types of anaerobic digesters each of whichhavespecificcharacteristicsmakingthemmoreadequateto treatspecifictypesoforganicwastes (Maoetal.,2015).Fortreat- ment of solid waste andsludges, low-rate anaerobicsystems are more appropriate due totheir useof long butcoupled hydraulic andsludgeretentiontimestoensureastableoperationofthepro-

∗ Corresponding author.

E-mail address: [email protected] (F. Abunde Neba).

cess(Mesetal.,2003).Ontheotherhand,amajorbreakthroughin theanaerobictreatmenttechnology hasbeenthedevelopmentof high-ratesystems,whichusebiomassretentiontoemployshorter hydraulicretentiontimes,butthistechnologyismainlyadaptedto thetreatmentofwastewaters(Henzeetal.,2008).Designofhigh- ratesystems for wastewatertreatment has received considerable attentionover thepast years anda varietyof novelor improved digesterdesignsandhydrodynamicconfigurationshavebeenpro- posedinthe literature(Zhangetal., 2016; Mao etal., 2015). The motivationfordesigningnoveldigesterconfigurationshasbeento increase process stability, simplify construction and operation as wellasimproveprocesseconomics.However,fewresearchershave beenable to drawon anysystematicresearch intothe improved designandoperationstrategiesoflowrateanaerobicreactorsand theuseoflongprocesstimes(whichislinkedtoeconomicfeasibil- ityofasystem)remainsachallengetosuchsystems.Bothefficient andeconomicalperformancesof low-ratedigestersare extremely importanttopromotetheir widespreadadoption fortreatment of

https://doi.org/10.1016/j.compchemeng.2019.106607

0098-1354/© 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )

Nomenclature

Bt Annualsavingsfromelectricityconsumption($) C_Gen Costofbiogaselectricitygenerator($)

C_Inv Investmentcost($)

Ccon Costofdigesterconstruction($) Cm Annualcostofdigestermaintenance($) C_misc Costassociatedtomiscellaneousactivities($) C_pf Annualcostofbiogaspurification($) Ct Annualoperatingcost($)

K∗ KineticconstantofChenandHashimotobasedon specificmethaneproductionrate(−)

K_ha Maximum rate of substrate degradation by the acidogenicbacteria(day⁻¹)

K_Lr Monod-like half saturation constant for continu- ousmodeoperation(gVS/L/day)

Kam Maximumrateofsubstrateutilizationbytheace- togenic/methanogenicmicroorganisms

K_i Inhibitionconstantforbiogasyield(g/L)

Ks Monod-like half saturation constant for continu- ousmodeoperation(gVS/L)

Lr Organicloadingrate(gVS/L/day)

P_el Percentageofmethaneutilizedforelectricitygen- eration

Pr_pf Price for biogas purification per unit volume ($/m³)

Pt Netannualbenefit($) R² coefficientofdetermination

R²Adj Adjustedcoefficientofdetermination

R_f Recalcitrant fraction of initial volatile substrate thatisnon-biodegradable.

Rp Specific rate ofbiogas productionby acetoclastic methanogens(mL_biogas/gVS/day)

Rpm Maximum specific rate of biogas production by acetoclasticmethanogens(mL_b/gVS/day)

S_i Initialconcentrationofsubstratetakenupbyace- togenic/methnogenicmicororgansims(g/L) So Initialsubstrateconcentration(gVS/L)

S_u Concentration of acidified substrate produced by acidogenicbacteria(g/L)

T_el Annualtimeperiodforuseofelectricity(d) Tpr Annualtimeperiodforbiogaspurification(d) V_D Volumeofdigester(mL)

X^TX Characteristicmatrix

Xam concentration of acetogenic/methanogenic mi- croorganisms

Y_PS Biogasyieldcoefficient(mL_biogas/gVS)/(gVS_utilized/L) Y_XS Cellyieldcoefficient

b_e Unitconversioncoefficient(kWh/m³CH4) dyt/dt Rateofchangeinbiogasyield

−dS/dt Rate of decrease in the concentration of hy- drolyzedsubstrate(gVS/L/day)

k_i Inhibitionconstantforcellgrowth(g/L)

k_n(s) Rateofsubstratedegradationbyacidogenicbacte- ria(g/g/day)

ks Monods’ half saturaion contant forsubstrate up- take

m_VS Mass of volatile solids added into the digester (gVS)

nt Projectlifespan(years)

rp Modifiedrateofbiogasproductionbyacetogenic/

methanogenicmicroorganisms(mL_b/gVS/day)

s_βˆ

i Approximate standard error of parameter esti- mates

t_υ,α/2 Studentt-distributionparameter(−) t_υ,α/2 student’st-distributionparameter yt Biogasyield(mL_biogas/gVS)

y_tm Maximumattainablebiogasyield(mL_biogas/gVS) (

β

−

β

ˆ) ^{Deviation between the real and the estimated}

modelparamters

β

ˆ ^{Vectorofestimatedmodelparameters}(−)

β

ˆ ^{Vectorofestimatedmodelparameters}

μ

^{max Maximum specific growth rate of aceto-}

genic/methanogenicbacteria(day⁻¹)

σ

^{2 Truevariance}

χ

^{2 Reducedchi-square}

AR AttainableRegions

B-RADeS BiodigesterRapidAnalysisandDesignSystem CSTR ContinuousStirredTankReactor

DSR DifferentialSidestreamreactor GUI GraphicalUserInterface

GUIDE Graphical User Interface Development Environ- ment

IDEAS InfiniteDimEnsionAlState-space NRT NetworkResidenceTime

NRT-C-AR Network Residence Time Constrained Attainable Region

PFR PlugFlowReactor

ν

^{Numberofdegreesoffreedom}

C Two-dimensionalstatevector

F F-distributionvarincecomparismparameter Pr_el Feed in tariff rate for biogas based electricity

($/kWh)

Q Volumetricflowrate(mL/day) RMSE Rootmeansquareerror S Substrateconcentration(g/L)

X Acidogenicbacteriaconcentration(g/L)

b Fraction ofinitial volatile solids remaining in ef- fluent.

cov(

β

^{) Covarianceofestimatedmodelparameters}(−) cov(

β

^{) Covariancematrixofestimatedmodelparameters}

gVS Gramvolatilesolids gVS Gramvolatilesolids

k KineticconstantofChenandHashimoto(−) m Measure of microbialadoption to stationarypro-

cessesbymutation(−)

n Couldprovideausefulmeasureofmicrobialcoop- erativity(−)

p Numberofmodelparameters r Discountrate(%)

r(C) Two-dimensionalreactionratevector

α

Significancelevel

β

^{Vectorofmodelparameters} (−)

β

^{Vectorofmodelparameters}

μ

^{Specific growth rate of} acetogenic/methanogenic bacteria(day⁻¹)

τ

^{Retentiontime(days)}

sludgesandsolidwastes.Usingempiricalmethodstooptimizede- signofanaerobicdigestersoftenrequiresconstructionofexpensive prototypesystemsandtimeconsumingstudies,whichhasbeena key motivationforrelianceonmodel-basedtechniques (Yu etal., 2013).Use ofthemodels isagainhighlydependent onthe avail- ability of kinetic coefficients and hence modelling requirements for design of anaerobic digesters are often simplified to a mini- mal number of inputs and experimental states (most commonly

biogas yield) required for model identification (Batstone, 2006).

The simplified models employed for digester design can be sin- gle stage, basedon theratelimitingstepapproach, ortwo-stage, based on lumping the process in to acid-forming and methane- producingmicroorganism.Thesinglestage modelspreviouslyem- ployed indigesterdesignincludefirstordermodels (Linke, 2006; MomohandNwaogazie,2011)andthemodelsbasedonmaximum bacteria growthrate(Fdez.-Güelfo etal., 2011; Fdez-Güelfoetal., 2012) while thetwo stage models includethe biogas yieldmod- elspresentedbyMomohetal.(2013).Althoughsimplifiedmodels have been used extensively in digester design,published articles are limitedtomainly todetermination ofdigestercapacitybased on parameters such a VS loading,temperature, etc. with firstor- der models beingmostly used(Momohet al., 2013; Wang etal., 2007).Fordesign,criticalissuesarehydrodynamics,aswellasthe behaviourofsolids,whichrequiresatleasttwo-stage,withhydrol- ysis and biological steps (Batstone, 2006). This studyfocuses on hydrodynamicsandanapproachtooptimizehydrodynamicconfig- uration ofanaerobicdigestersbasedon simplifiedmodels willbe a highlypracticableasit willrequirelessexperimental measure- ments toestimate kinetic constants. Althoughthestudydevelops two-stage models which can be identified based using only bio- gas yieldmeasurements, theemphasis ofthepaperis notneces- sarily on the models buton how the authors usethe models to develop new hydrodynamicconfigurations for operatinglow rate anaerobic digesters.The design objectiveis to minimize the pro- cesstimeaswellaspaybackperiodbyconsideringbiodegradation andmixingastheonlypermittedfundamentalprocessesoccurring inthedigester.Theapproachisbasedontheconceptofattainable regions (AR), which is a techniquefor process synthesis andop- timization that incorporates elementsof geometryto understand hownetworksofchemicalreactorscanbedesignedandimproved.

“Following thisinitialwork, manyother researchersadvanced AR research.Glasseretal.(1987)proposedageometricapproachthat identified candidate AR’s satisfying a number of necessary con- ditions that the AR must possess. Burri et al. (2002) demon- strated that, within the InfiniteDimEnsionAl State-space (IDEAS) conceptual framework, construction of the true AR, and increas- ingly accurate AR approximants, can be carried out through Infinite Linear Programming (ILP), and a sequence of approx- imating finite Linear Programs (LP) respectively. Subsequently, Manousiouthakis et al. (2004) developed, within IDEAS, neces- sary andsufficient conditionsthat thetrue AR mustsatisfy, pro- posedtheShrink–Wrap algorithmforARconstruction,andestab- lished,thisalgorithm’sequivalencetotheaforementionedLPbased AR construction methods. They also demonstrated that the true AR can be potentially larger than the candidate AR’s identified bygeometricmethods.ZhouandManousiouthakis(2006)demon- strated that the true AR for reactor networks involving only re- action and mixing may be smaller than the true AR for reactor networks also incorporating diffusion effects (e.g. by considering non-ideal dispersion reactor models). Zhou and Manousiouthakis carried out pollution prevention studies using the AR approach (Zhou and Manousiouthakis, 2007a), extended the AR approach to reactor networks involving variable density fluids (Zhou and Manousiouthakis, 2007b), discussed the dimensionality of the space in which AR construction can be pursued (Zhou and Manousiouthakis, 2008), and extended the AR approach to non- isothermal reactor networks (Zhou and Manousiouthakis, 2009).

Around thesame time,Posada andManousiouthakis(2008), pro- posed AR construction methods for reactor networks with mul- tiple feeds, while Davis et al. extended the AR approach to batch reactor networks (Davis et al., 2008). More recently, Con- ner and Manousiouthakis extended the AR approach to gen- eralprocess networks(ConnerandManousiouthakis,2014),while Ming etal.(2016)summarizedmanyoftheARliteratureresults.”

Geometrically,theattainableregionrepresentstheregionbounded by the convex hull for the set of points achievable by the fun- damental processes occurringin the system(Asieduet al., 2015; HildebrandtandGlasser, 1990; Hildebrandtetal., 1990).Oncethe ARhasbeendetermined,thelimitsofachievabilitybythesystem forthegivenkineticsandfeedpoint isknownandthe boundary oftheAR canthenbeusedtoanswerdifferentdesignand/orop- timizationquestionsrelatedtothesystem.Ourrecentpublication, AbundeNeba etal.(2019) hasbeen firstof itskind laying down theoretical framework foruse ofattainable regions to model op- timal configurations of multistage anaerobic digesters. The study employedafour-statedynamicmodelofanaerobictreatmentpro- cess and the attainable region analysis has been based on con- centration (state) space but not residence time. The lack of res- idence time makes it impossible to size the digester structure or perform economic feasibility studies on the optimal digester structure.Inaddition,thefour-state model(comparedtothesim- plified model in this study) poses requirement for more experi- mentalmeasurements hence limitingits application tosituations whereprocessmeasurementsarelimited.Thecurrentstudyisde- signed to illustrate how simplified models (requiring onlybiogas yield measurements) of the anaerobic treatment process can be usedforattainableregionanalysisinvolvingresidencetimespace.

“AlHusseini and Manousiouthakis (2013) were the first to incor- porateresidence-time considerationsintheAR conceptualframe- work, by first introducing a production normalized, capital cost measure fora reactor network, that they termed “NetworkResi- denceTime” (NRT), and definedas “the ratio of the sum of the volumesof all reactors participating in the reactor network over the total volumetricflowrate entering the network.” They subse- quently introduced the Network Residence Time Constrained At- tainableRegion(NRT-C-AR),whichtheythenproceededtoquantify usingaLinearProgrammingformulationwithintheIDEASconcep- tualframework.TheadvantageofconstructinganARoftheanaer- obictreatmentprocess that incorporatesresidencetimeconsider- ations,overtheAR presentedinourprevious study,isthatiten- ablesthecouplingofbiodegradation kinetics,economic feasibility objectivesandcountryspecificmacroeconomicparametersforthe synthesisofbiogasdigesterstructures.Byitsuseofattainablere- gions,knowledgeofallpossiblestates,forallpossibledigestercon- figurationscanbeobtainedconsideringbiodegradationandmixing astheonlyfundamentalprocessesoccurringinthedigester.Unlike previous studies where economic analysisis performed to deter- mine the feasibility parameters of a predefined digester configu- ration,thisstudyratherdetermines the biodigesternetwork con- figuration required to achieve a given economic objective based on the macroeconomic situation of a given country. Finally, the studyseeks to deploy the theoretical framework into a software in order to save time and effort for designers who are plan- ninganddesigningbiogasplantsfordifferentprocessoreconomic scenarios.

2. Attainableregiontheoryforprocesssynthesisand optimization

TheAttainableRegion(AR)theoryisa techniquethat incorpo- rateselementsofgeometryandmathematicaloptimization,tode- signandimproveoperationofchemicalreactors(Mingetal.,2016).

The power of the AR approach to process optimization is that the answer to all possible optimization problems, even the ones notconsidered are firstdetermine,andthen we lookforwaysof achievingthatanswer.Inreactoroperationknowledgeofallpossi- blereactorstatesforallpossiblereactorconfigurations,eventhose thathavenotyetbeendevised,isobtained.Foratwo-dimensional system, theconvexhullfortheset ofallpointsachievable byall possiblecombinationsofCSTR+PFRandmixingdefinestheattain-

ableregion.Forhigherdimensionalsystems,theattainableregion is the convex hull for the set of points generated by all possi- blecombinationsofCSTR,PFR,DSR andmixinglines. Theconvex hullis understood as the smallestsubset of a set of pointsthat can be used to generateall other pointsby reactionand mixing (Ming etal., 2016). Geometrically,aconvexhullisafiniteconvex polytopeenclosedbyafinitenumberofhyperplanes,whichisin- terpretedinatwo-dimensionalspaceasthesmallestpolygonen- closedby planar facetssuch thatall of theelementslie on orin theinteriorofthepolygon(Asieduetal., 2015).OncetheAR has beendetermined,thelimitsofachievabilitybythesystemforthe givenkinetics andfeedpoint is known, whichcan then be used toanswerdifferentdesignoroptimizationquestionsrelatedtothe system.

Given a setof reactions andassociated kinetics,thefollowing fivekeysteps needstobe performedinordertocompletean at- tainableregionanalysis(Mingetal.,2016):

➢ Definethereaction,dimensionandfeedset.

➢ Generate the AR using combinations ofthe fundamental pro- cesses.

➢ InterprettheARboundaryintermsofreactorequipment.

➢ Define the objectivefunction andoverlay this onto theAR to determinepointofintersectionwiththeARboundary.

➢ Determinethespecificreactorconfigurationrequiredtoachieve theintersectionpoint.

Some necessary conditions for AR derived from the work of Glasseretal.(1987)canbesummarizedasfollows:

➢ TheARincludesallfeedstothesystem.

➢ TheARisconvex.

➢ NoprocessvectorpointoutoftheARboundary.

➢ No rate vectorsinthe complementof the AR whenextended backwardintersectstheAR.

The objective ofthissection is toanalyze theaforementioned necessary requirements with respect to its application to the anaerobic treatment process. However, AR analysis requires that theprocess kinetics isknown andwe thereforebegin by model- ingthekineticsoftheanaerobictreatmentprocess.

2.1.Reactionkineticsoftheanaerobictreatmentprocess

In the presentpaper,the mathematical models describingthe kineticsofsubstrateutilizationandmethaneproductioninanaer- obictreatmentprocessaredevelopedbasedontheapproachpre- sented by Momoh et al. (2013). The approach assumes that the ADprocess takesplaceinthree stages.(i)hydrolysis/acidogenesis of the organic substrates in wastewater by acidogenic bacte- ria to produce acidified substrate; (ii) uptake of acidified sub- stratebyacetogenic/methanogenicbacteriaand(iii)acidifiedsub- strate assimilation, growth and biogas production by the aceto- genic/methanogenicbacteria.

2.1.1. Enunciationoftheprocessmodel

Fig. 2 presents the algorithm used to develop and validate thesimplified two-stagebasedmodestopredictbiogasyield.The modeldevelopmentinvolvesfivemainaspects,whichinclude:

➢ Developsubstratedegradationmodel.

➢ Formulatesubstrateuptakemodel.

➢ Choosemicrobialkineticmodel.

➢ Derivemodelsforsubstrateassimilationandbiogasproduction.

➢ Identificationofthedevelopedmodel.

Considering theseaspects led to a seriesofordinary differen- tialequations to predict biogas yield based on microbial growth kinetics

Stage1:Hydrolysisandacidogenesis

Many constituents of organic wastes behave ascomplex sub- strates(polysaccharides, proteins, fatsetc.). The Grau model pre- sentedinEq.(1),whichhaswidely beenused tomodelmultiple substrateremovalkinetics(Kimetal.,2006;Liu,2006)wasthere- foreadoptedforthisstudy.

−dS

dt =kn(^s)X

_S

So

ⁿ

(1)

where −dS/dtrepresentsthe rateof decrease inconcentration of substrate beinghydrolyzed, S is the concentration of initial sub- strate left atevery instant following onsetof hydrolysis,Soisthe concentration of initial substrate, k_n(s) is the rate of substrate degradationbyacidogenicbacteria,Xistheconcentrationofacido- genicbacteriaandndefinesthedegreeofadaptationbyacidogenic bacteriaforsubstratedegradation.

The multicomponentsubstratedegradation model isbased on theassumption thatthe differentcomponentsare simultaneously removedandtransportedintothecells(Grau etal.,1975).Assum- inghydrolysisandacidogenesisarecatalyzedbyacidogenicbacte- ria,whoseconcentrationisconstant,thenEq.(1)canbere-written asEq.(2).

−dS dt =Kha

_S

So

ⁿ

(2)

Where,K_ha isthe maximumrateofsubstratedegradationby aci- dogenic bacteria.Since anaerobicdigestionisa biological process andtheAR approachconsiders biodegradationandmixingasthe only fundamentalprocessesoccurring inthe digester, itbecomes importanttoconsiderthenon-biodegradablepartofthesubstrate.

ThemodelisthenmodifiedasshowninEq.(3)

−dS

dt =kn(^s)X

_S

So−R_f

ⁿ

(3)

Eq.(3)representsthekinetics ofsubstratedegradation,where R_fistherecalcitrantfractionofinitialvolatilesubstratethatisnon- biodegradable.

Stage2:Substrateuptakebyacetogenic/methanogenicmicroor- ganism

The hydrolytic model of Momoh et al. (2013), Eq. (4), which represents a modified version of the hydrolytic model presented bypreviousstudies(Barthakuretal.,1991;FaisalandUnno,2001; Zinatizadehetal.,2006)wasadopted.

Su=SoAf

b−Rf

n

(4) Thismodeltakesintoconsiderationtheacidifiedsubstratepro- ducedaftersubstratedegradationbyacidogenicbacteriaaswellas the uptakeofacidified substrateby acetogenic/methanogenicmi- croorganism.Surepresentstheactualamountofthesubstratethat wasacidified andutilized bythe acetogenic/methanogenicbacte- riawhilebisthefractionofinitialvolatilesolidsremaininginef- fluent. The coefficient A_f=K_ha/(^Kam(¹−

α

)+K_ha) ^{represents the} rate limiting coefficient for very slow (case of 0<

α

<1) or very fast(case of

α

=1) metabolism ofacidified substrateby the ace- togenic/methanogenicbacteria(Momoh etal., 2013).The constant Kam is the maximum rate of substrate utilization by the aceto- genic/methanogenicmicroorganisms.

Stage3:Kineticsofbacteriagrowthandbiogasproduction The attainable region is unique for a given kinetics and a change in organic substrate can cause a change the ki- netic model used to describe the growth of microorganisms.

Table 1presents a listof microbial growthmodels considered to modelsubstrateassimilation. Thetable hasbeen assembledfrom Kythreotou et al.(2014),who presented a comprehensive review ofsimpletoscientificmodelsforanaerobicdigestion.Asexpected, the different models have different characteristics often making

Table 1

Microbial growth models selected to model substrate assimilation.

Author Model equation Eq. No. Remark

Monod, 1949 μ⁼ μmax Su

ks+Su (5) Describe growth processes for low substrate concentration

Moser, 1958 μ⁼ μmax Sum

ks+Sum (6) Integrates effect of microbial adoption to stationary processes by mutation Tessier model μ⁼ μmax(1 −e ^−Skus) (7) An exponential function used to describe cell growth processes

Chen & Hashimoto, 1978 μ= _kS^μmaxSu

i+(1−k)Su (8) Considers cell concentration depending on the level of substrate degradation Haldane, 1930 μ^{= μmaxSu}

(ks+Su)(1+S u/k i) (9) For growth process affected by the allosteric effectors present in the acidified substrate

Andrews, 1968 μ= ^μmaxSu

ks+Su(1+S u/k i) (10) Based on Haldane for enzyme inhibition at high substrate concentrations Aiba et al., 1968 μ= ^μ_k^max_s+S^S_u^uexp(−S u/ k i) (11) An empirical correlation to describe substrate inhibition

Dagley & Hinshelwood, 1983 μ= μmax Su

Ks+Su(1 −k iS u) (12) An empirical correlation with critical inhibitor concentration of growth stop Ierusalimsky, 1967 μ⁼ μmax Su

ks+Su ki

ki+Su (13) Haldane model for product inhibition Moser, 1981 und Bergter, 1983 μ= μmax Sum

ks+Sum ki

ki+Su (14) Production inhibition model with effect of microbial adoption

them more adequate to describe microbial growthof specific ef- fluentand/ordigesterconditionsratherthanothers.

k_s isthe Monods’ half saturaion contant forsubstrate uptake,

μ

^{max is themaximumspecific growthratefor} methatnogenicar- chae, m is the coefficient of acetogenic/methanogenic microbial adaptation forcooperativity,S_i isthe initialconcentration ofsub- strate taken up by acetogenic/methnogenic micororgansims, k is the kinetic constant of Chenand Hasshimoto, and k_i is the sub- strate concentration where bacteria growth isreduced to 50% of themaximumspecificgrowthrateduetosubstrateinhibition

Taking the case of the Monod model for growth of aceto- genic/methanogenic microorganisms, and using product and cell yield coefficients, the rateofbiogas productioncan be expressed byEq.(15)

dyt

dt =YPS

YXS

μ

maxSu

ks+Su

Xam (15)

Wheredyt/dtistherateofchangeinbiogasyield,yt isthebio- gasyield,Y_PS isthebiogasyieldcoefficient,Y_XSisthecellyieldco- efficientandXam istheconcentrationofacetogenic/methanogenic microorganisms. If we consider the growth rate of the aceto- genic/methanogenic bacteria is very slow or relatively constant while dyt/dt can be described as the specific biogas yield rate (Rp) at the end of biogas production (Momoh et al., 2013), then Eq. (15) can re-written as Eq. (16). The parameter R_pm= (XamY_PS/Y_XS)isthemaximumspecificrateofbiogasproduction.

Rp= RpmSu

ks+Su

(16)

Hence,bysubstitutingEq.(4)intoEq.(16)andrearranging,we obtainEq.(15).

Rp= RpmSo ks Af

(

^b−Rf

)

^n+So

(15)

The termks/A_f(^b−R_f)^{n representsthe Monodhalf saturation} constantintermsofthefractionofacidifiedsubstratetakenupby acetogenic/methanogenicbacteria(representedasK_S)andthefinal biogasyieldmodelconsideringtheMonodkineticsispresentedby Eq.(16).

Rp= RpmSo

KS+So

(16)

Eq.(16)describesthebiogasyieldinanaerobicdigesterconsid- eringabatchoperationmode.Incaseswherethesysteminoper- atedincontinuousmode,theinitialsubstrateconcentration(So)is convertedtoloadingratebymultiplyingthefactor(Q/V)asshown

Table 2

Two-stage based models to predict biogas yield rate.

Author Model equation Eq. No.

Monod, 1949 R p= R pm So

KS+So (19) Moser, 1958 R p= R pm Som

KS+Som (20) Tessier model R p= R pm(1 −k pe ^−SKoS) (21) Chen & Hashimoto, 1978 R p= _K ^RpmSo

∗Si+(1−k)So (22) Haldane, 1930 R p= ^RpmSo

(Ks+So)(1+S o/K i) (23) Andrews, 1968 R p= R pm So

Ks+So(1+S o/K i) (24) Aiba et al., 1968 R p= R pm So

Ks+Soexp(−S o/ K i) (25) Dagley & Hinshelwood, 1983 R p= R pm So

Ks+So(1 −K iS o) (26) Ierusalimsky, 1967 R p= R pm So

Ks+So Ki

Ki+So (27) Moser, 1981 und Bergter, 1983 R p= R pm Son

Ks+Son Ki

Ki+So (28)

inEq.(17).Thefactor(Q/V)istheratioofvolumetricflowrate(Q) tovolumeofthedigester(V).

Rp= Rpm

(

^QSo/V

)

(

^KsQ/V

)

⁺

(

^QSo/V

)

⁽¹⁷⁾

Theresultingcontinuousmodecounterpartofthebiogasyield modelconsideringMonodkineticsisshownbyEq.(18).

Rp= RpmLr

KLr+Lr

(18)

Where,thevariableLristheorganicloadingrateintothebiodi- gester andKLr is the Monod’shalf saturation constant definedin termsof theorganic loadingrate. Similar process wasapplied to develop theother biogas yield ratemodels by assuming that the growthprocessoftheacetogenic/methanogenicmicroorganismcan bedescribedusingtheothergrowthmodelspresentedinTable1. However, a parameter of kp wasintroduced to the Tessier based modelasacoefficienttotheexponentialterm,whichservesasan indexoftheprocessingspeedofR_pasitapproachesR_pmduetothe changeinSo orLr (dependingon themodeofoperation).Thede- rivedbiogas yieldmodels considering atwo-stagebiodegradation kineticsispresentedinTable2.

K is the kinetic constant of Chen and Hasshimoto defined in terms of specific biogas yield, K_i is the substrate concentration wherespecificbiogasyieldrateisreducedto50%ofthemaximum specificbiogasyieldrateduetosubstrateinhibitionandnprovides ausefulmeasureofmicrobialcooperativitytobiogasproduction.

2.1.2. Parameterestimationandstatisticalmethods

The kinetic constantsof the different models were estimated using the Matlabnonlinear regression solver ‘nlinfit’ (Mathworks NatickNA).Inassessingthevariabilityofthemodelidentification process,we usedthe kerneldensityestimatesandtheparameter confidenceregions.Itisalsointerestingtonotethatmarginalcon- fidenceintervalsoftenusedbyseveralresearchersdonotaccount forcorrelationsbetweentheparameterestimates.Therefore,their useinparameterestimationcansometimesbemisleadingifthere isstrongcorrelation betweenseveralparameter estimates.Inthis study,we rather illustrate the use ofjoint confidence regions in assessingreliabilityofparameterestimatesinleastsquare regres- sion.

The100(¹−

α

)^{%jointconfidenceregionandthemarginalcon-} fidence intervals of the parameter estimates is computed using Eqs.(29)and(30)respectively

β

−

β

ˆ

T

X^TX

β

−

β

ˆ

≤p

σ

^2F(¹−α),p,(ⁿ−p) (29)

β

ˆ^±tυ,α/2s_βˆ

i (30)

Where s_βˆ

i isthe approximatestandard erroroftheparameter estimatesgivenbyEq.(31)

s_βˆ

i=

diag

(

^co

v ₍ β ) )

⁽³¹⁾

X^TX=s_βˆ

i

2/co

v

₍

β

)^isthecharacteristicmatrix,cov(

β

^)istheco-

variancematrixofestimatedmodelparameters,(

β

−

β

ˆ)^isthede- viationbetween thereal(

β

^{) andthe estimatedmodelparamters} (

β

^ˆ),t_υ,α/2isthestudent’st-distributionparameter,

ν

^isthenum-

berofdegreesoffreedom (n−1),wherenisthenumberofdata pointsused tocompute thevariance andaverage,F_{(1−α),p,(n−p)} is theF-distribution varincecomparism parameter,p isthe number ofmodelparametersbeingestimated,(ⁿ−p)^{isthemodeldegrees} offreedom,and

σ

^{2isthetruevariance,and}

α

=0.05isthelevel ofsignificance.

2.2.ARanalysisoftheanaerobictreatmentprocess

2.2.1. Reactiondimensionandprocessvectors

Before it is possible to construct the AR, the engineer must firstdeterminethespacewhereintheARmustreside(bychoosing uniquespeciescomponents inthe systemthat willrepresentthe AR). Theseare variables requiredto characterize the state ofthe system,inthiscase, ananaerobictreatment processandmustbe sufficientto describe the dynamicsof thefundamentalprocesses chosentodescribethesystem. Thesevariableswouldincludebio- gasyield (y_t) and retention time (

τ

^{). A key criteriafor selecting}

variablesinARisthat theymustobeythelinearmixinglaw.The conceptof“NetworkResidenceTime” (NRT),asintroduced byAl- Husseini and Manousiouthakis (2013) defines the residence time of a reactor network asthe ratio of the sum of the volumes of all reactor units constituting the network to the total volumet- ricflowrateentering the network.Using thisdefinition,it canbe shownthat the residence time of a network comprising two re- actorunitsobeythe linearmixinglaw, Eq.(32).Thisimplies the overallresidencetimeofthesystemmustlieinastraightlinebe- tween the residence times of the individual reactors,

τ

1 and

τ

2

comprisingthesystem.

τ

mix=

ατ

1+

(

¹−

α ) τ

2 (32) Where

τ

mixistheoverallresidencetimeofthesystemcompris- ingtwoindividualreactors.Thedevelopedsimplifiedkineticmod- elspredictbiogasyield(yt),whichisgivenintermsofvolumeof biogasproduced(mL)pergram ofvolatilesolidsaddedtothedi- gester(gVS).yt=Vg/m_VS.Assumewehavetwodigestersofknown

biogasyield,we canobtaintheactualvolume ofbiogasproduced for digesters 1 and2 byV_g1=y_t1m_VS1 andV_g2=y_t2m_VS2 respec- tively.Conservationofmassmaybeusedtocalculatethetotalbio- gasyieldforbothdigesters.Conservationofmassensuresthatthe totalmassofvolatilesolidsinthemixtureisequalto thesumof theindividualmassescontainedinbeakers1and2,whichisgiven by m_VST=m_VS1+m_VS2. Computingthe biogas yield of theentire systemisequivalenttodeterminingthebiogasyieldforamixture ofdigesters 1and2 since the densityofthe liquid phase ofthe digestercanbeassumedconstant.Thebiogasyieldofthemixture (y_tM)isgivenbytheratioofthetotalvolume ofbiogasproduced tothetotalmassofvolatileacidsaddedasshownbyEq.(33). ytM=yt1mVS1+yt2mVS2

mVST

(33)

If we set

α

=m_VS1/m_VST then Eq. (33) can be written as Eq. (34), which is similar to the linear mixing law. What this meanspracticallyisthatifwemixthecontentsoftheliquidphase oftwodigesters,eachofwhichcontainsagivenquantityofvolatile solidsadded,then thetotalbiogasyieldofthemixturewillliein astraightlinejoiningthatofbothdigesters.

ytM=

α

^yt1+

(

¹⁻

α )

^{yt2 (34)}

Theprocessofcombiningthecontentsoftwoparalleldigesters (ordigesternetworks)ofdifferentvolatilesolidscontentsresultsin alinearmixinglawmeasuredintermofbiogasyield.Thisimplies biogasyieldmaybeusedintheconstructionofcandidateARsina similarmannertothatforconcentration.

Thebiogasyieldandtheretentiontimegroupedtogetherform avectorcalledthecharacteristicvector;C=[yt

τ

^],whosedimen- sion determines the dimension of the optimizationproblem. We therefore havea 2-D optimization problemwith the objectiveof minimizing[

τ

^{],timeparameter.}

Ifwe assume that assubstrateis consumedrateof changeof biogas yield is directly correlated withthe quantity of biogas to thebiogasyieldy_t,suchthatthedrivingforceforgasproductionis disappearingwhenthebiogasyieldgradually approachesitsmax- imum (ytm) then for themass ofvolatile solidsadded to the di- gester.ThisismodelledbyEq.(35)

rp=dyt

dt =Rp

1− yt

ytm

(35)

Where, rp is the modified rate of biogas production by ace- togenic/methanogenicmicroorganisms. The reactionratevector is thereforegivenby−−→

r(^C)=[rp r_τ]

2.2.2. GeneratetheARusingcombinationsofthefundamental processes

The attainable region (AR) represents the set of all possible states that can be achieved by a combination two fundamental processes,biodegradationandmixinginthecaseoftheanaerobic treatmentprocess.InARtheory,mixingisperformedbyacontin- uousstirredtankreactor(CSTR)whilebiodegradation(reaction)is performedby the plug flow reactor (PFR), sincethe operation of bothreactors respectivelymimicsthetwo fundamentalprocesses.

At steady state operation, the general mathematical model of a CSTRandPFRaregivenrespectivelybyEqs.(36)and(37) respec- tively.

C=C_f+

τ

^r

(

^C

)

⁽³⁶⁾

dC

dx =r

(

^C

)

⁽³⁷⁾

Cisthetwo-dimensionalstatevectormadeofbiogasyieldand theresidencetime, Eq.(38)whiler(C)isthe reactionratevector, whichcanbeillustratedtobegivenbyEq.(39).

C=[yt

τ

^{]T (38)}

r

(

^C

)

^{=[rp 1]T (39)}

During constructionofAR, thePFR trajectory istheset points generatedbynumericallysolvingthePFRequationwhiletheCSTR locusisthesetofpointsgeneratedbysolvingtheCSTR equation.

The convex hullfor the set ofall points achievable by all possi- blecombinations ofCSTR+PFR definestheattainable region.The convexhullisunderstoodasthesmallestsubsetofasetofpoints thatcanbeusedtogenerateallotherpointsbyreactionandmix- ing(Mingetal.,2016).Geometrically,aconvexhullisafinitecon- vexpolytopeenclosedbyafinitenumberofhyperplanes,whichis interpreted in a two- dimensional space as the smallestpolygon enclosedbyplanarfacetssuchthatalloftheelementslieonorin theinteriorofthepolygon(Asieduetal.,2015)

The candidate attainable region wasconstructed with Matlab usingthefollowingfive-steps

Step1:DeterminePFRtrajectoryfromfeed Step2:DeterminetheCSTRlocusfromfeed

Step3:DeterminePFRtrajectoryfromeachCSTRpoint Step4:Constructtheconvexhullofthesetofachievablepoints Step5:VerifytheobtainedARagainstthenecessaryconditions ofARandifanyconditionisnot metreturnextendtheAR byrunningaPFRfromthepointofdisagreement

The PFR equations are solved using theMatlab ode45 routine forsolvingnon-stiff differentialequationswhilethesystemofnon- linearCSTR equationswere solved using‘fsolve’routine The con- vexhalloftheentiresetofgeometricpointsisobtainedbyusing theMatlab‘convhull’routine,whichimplementstheQuickhullal- gorithm(Mathworks,NatickNA).

It is importantto mentionthat even though thisconstruction approachhasbeenappliedinacoupleofstudies(Mingetal.,2013; Ming etal., 2016),the NRT-C-AR obtainedisa candidateandnot thetrue NRT-C-AR.ForatrueAR, theInfiniteDimEnsionAlState- space(IDEAS)conceptualframeworkisappliedtoobtainageneral linear programming formulation for the construction of the true NRT-C-AR, as shown in (AlHusseini and Manousiouthakis, 2013).

However,theinterestofthestudyisnotnecessarilyonthemethod usedforARconstruction,butonhowtheconceptofattainablere- gionscanbeappliedtooptimizedoperationoftheanaerobictreat- mentprocess. Also,evenifjustacandidateAR isobtained,itcan still be usedforprocess synthesisandoptimizationonly thatthe totalityofoutputsisnotobtained.

2.2.3. InterprettheARboundaryintermsofreactorequipment TheARboundaryiscomposedoftwotypesofgeometries:mix- inglinesreferredtoaslineationsandmanifoldsofPFRtrajectories referredtoasprotrusions.TheroleofPFRsontheARboundaryis togeneratetheouter extremitieswhereasCSTRsandDSRs(inthe caseof higherdimensional constructions)are used asconnectors to thesePFRs (Ming et al., 2016). This implies the AR boundary isdefinedintermsofreactorstructures, andfortwo-dimensional constructions, the boundary is composed of combinations of the twofundamentalreactortypesandmixinglines.ThePFRandCSTR eachexhibit uniquegeometricinterpretationandhencedetermin- ing the reactor configurations that form the AR relates to inter- pretingthesurfacesthatformtheARboundaryusingitsgeometric properties.

2.2.4. Definetheobjectivefunctionandtheoptimalreactorstructure TheAR,whichdefinesthelimitsofachievabilitybythesystem forthegivenkineticsandfeedpointcanbeusedtoansweroneor

moredesignoroptimizationquestionsrelatedtothesystem.Two- dimensional ARsinvolving residencetime are particularly impor- tantforunderstanding the minimumreactorvolume requiredfor agivenoutput.Sincethe constructionandoperationofanaerobic digesters generallyrequirescapitalinvestment, itwouldbe inter- estingtousetheARconceptindeterminingtheprofitabilityofthe plant. However, we needto develop a suitable objectivefunction thatincorporatesbiogasyieldandresidencetime(ordigestervol- ume).

Theeconomic evaluationconsiders thatbiogas generatedfrom the anaerobic digester is utilized for electricity generation. The total annual income, benefit (Bt) from installing the biomethane plantis determined by Eq.(40), which is benefit dueto savings fromelectricityconsumption.

Bt=0.9P_el×T_el×be×Pr_el×gVS×yt (40) WhereP_elis thepercentageofmethaneutilized forelectricity generation(whichis100%inthecaseofthecurrentstudy),T_el is theannualtimeperiodforuseofelectricity,beistheunit conver- sioncoefficient,Pr_elisthefeed-intariff rateforbiogasbasedelec- tricity,gVSisthegrammassofvolatilesolidsfedinthedigester.

Thetotalannualexpensesoroperatingcost(Ct)iscomputedby Eq.(41).Theoperatingcosts areassumedtobeafunction oftwo factors:therepairandmaintenancecosts,Eq.(41a),whichistaken toby1%ofthecostofconstruction(0.01Ccon)andthecostofbio- gasupgrading,whichisafunctionofthebiogasvolume,Eq.(41b).

Ct=Cm+C_{p f} (41)

Cm=0.01Ccon (41a)

Cp f=gVS×yt×Tpr×Prp f (41b) Where,Cmistheannualcostofdigestermaintenance,C_pfisthe annualcostofbiogaspurification,Ccon isthecostofdigestercon- struction,T_pristheannualtimeperiodforbiogaspurificationand Pr_pfisthepriceforbiogaspurificationperunitvolume.

Thecostofinvestment/constructioniscomputedusingtherates ofacommercialbiogascompany inGhana,statingthecostofdi- gesterconstructiontobe$300percubicmeter(Mohammedetal., 2017). This includes administrative, transport costs, consultancy feesandotherlogisticaspects.Thefinalexpressionofthetotalan- nualcostofdigesteroperationisgivenbyEq.(42).

Ct=3V_D+gVS×yt×Tpr×Pr_{p f} (42) WhereV_D is the volume ofdigester. The annual profit(Pt), is definedasthedifferencebetweenannualbenefit,B_t,duetosavings fromelectricity consumptioin andthe annualoperating costs, Ct. ThisisexpressedbyEq.(43).

Pt=Bt−Ct (43)

SubstitutingtheexpressionsforB_t andC_t intoEq.(43)theex- pressionfortheannualprofitcanbewrittenasinEq.(44)

Pt=gVS×yt

0.9P_el×T_el×be×Pr_el−Tpr×Pr_{p f}

−3VD (44) SincetheARisconstructedinresidencetimespace,itisneces- sarytoexpressthe volumeofdigester(V_D) intermsofresidence time

τ

^{andvolumetricflowrateQ.TheexpressionforPtasafunc-}

tionofresidencetimebecomes;

Pt=gVS×yt

0.9P_el×T_el×be×Pr_el−Tpr×Pr_{p f}

−3

τ

^{Q (45)}

Theeconomicevaluationofthedigesterinvestmentisbasedon thepaybackperiod(PBP)(Gittinger,1986)andthedecisionruleis thatonegenerallyacceptsprojectsthat requireshorternumberof yearstorecovertheinvestment.Thepaybackperiodisgivenbythe annualprofits,generatedovernyears,neededtorecoverthetotal

Fig. 1. B-RADeS software development procedure.

Table 3

Summary of input parameters used for the economic evaluation.

S.N Parameter Unit Value

1 Discount rate % 10

2 Average cost of digester and infrastructure $/ m ³( Base ) 300 3 Biogas-based electricity generator (500 kW) $/4 PCS 600

4 Biodigester lifespan years 20

5 Upper calorific value of methane gas MJ / m ³ 39.8

6 Density of methane kg / m ³ 0.75

7 Electricity equivalent of methane kWh / m ³ 11.06 8 Feed-in tariff rate for biogas based electricity $/ kWh 17.5

futurecostofthebiogasplant,determinedusingthecompounded interestformula.ThepaybackperiodisevaluatedusingEq.(46).

ntPt=CInv

(

^1+r

)

^{nt (46)}

CInv=Ccon+CGen+Cmisc (46a)

WhereC_Invisthecostofinvestment,risthediscountrate,C_Gen isthecostofbiogasgenerator,C_miscisthemiscellaneouscostand n_t is the project lifespan. Eq. (46) can be rearranged to express yt asafunction of

τ

^{,Eq.(47) whichmaybe plottedovertheAR}

boundaryascontours(fordifferentspecifiedvaluesofn)todeter- minethepoint(s)ofintersectionwiththeARboundary.Thesein- tersectionpointsrepresenttheoptimaloperatingpoint(whichcan interpretedintoanoptimalreactorstructure)requiredtoachievea specifiedpaybackperiod.

yt

( τ )

^{= CInv}

(

^1+r

)

^nt+3

τ

^Q

ntgVS×

0.9P_el×T_el×be×Pr_el−Tpr×Pr_{p f}

⁽⁴⁷⁾

Table 3 presents of summary of the parameter sets that are usedtoperformtheeconomicevaluationofdesigningaconstruct- ingamethaneplant.

3. Developmentofcomputationalmodel

Thedesignofthegraphicaluserinterface(GUI)wasdoneusing theMatlabGUIDE(GraphicalUserInterfaceDevelopmentEnviron- ment).Thisisdoneusingicon-basedprogramingusingseveralob- jectssuch aspush buttons,statictexts, edittexts, pop-upmenus

andaxeshandles.GUIDEgeneratesaGUIandthem-filethatcon- tainsthecodetohandletheinitializationandlaunchingoftheGUI.

After creationofthe GUI,it wasprogrammedbyentering theal- gorithmsintothe variouscallback functionsintheMatlabm-file.

ThestepsofcreatingtheB-RADeSGUIinMatlabareshowninthe flowchartinFig.1.

4. Resultsanddiscussion 4.1. B-RADeSuserinterface

Fig. 2presents theGraphical User Interfaceof theBiodigester RapidAnalysisandDesignSystem(B-RADeS).Thismulti-levelpro- cessdesignandsimulationtoolcanbeusedtofindthemosteffi- cientdesignofmulti-stageanaerobicdigesternetworkstoachieve adefinedeconomicandprocessobjective.B-RADeShasseveralat- tributes that make it useful for a scientifically and economically objectiveprocessdesignandanalysisplatformforusebyengineers todotheircalculationsduringdesignandoperationofmulti-stage digesters.ThemainfeaturesofB-RADeSareasfollows:

➢ B-RADeS is based on peer-reviewed models that describe growth kinetics of anaerobic digestion microorganisms. It in- cludes ten simple biokinetic models derived based on biogas yieldanalogy.

➢ It does not rely on published kinetic coefficients, but it in- cludesasectionwheretheuserdetermineskineticcoefficients requiredfordigestersynthesisfromownexperiments.Uponin- putofexperimentaldata,B-RADeSautomaticallyscansthrough the10 modelsandranks them inorderof bestfit usingboth quantitativeandqualitativetechniques.

➢ Data requirements are simple: Only experimental measure- mentsofbiogasyield arerequiredtodeterminekinetic coeffi- cients,constructattainableregionsawellassynthesizedigester networks.

➢ It takes into account the country-specific macroeconomic pa- rameters(interestrate,electricityfeedintariff rateandannual workingdays) intothedesign process,whichisa keymotiva- tionforinvestors.

➢ It isbased ona systematicmethodological framework forthe design of multistage digester networksusing the global opti-