DepartmentofChemicalandBiomolecularEngineering

CN4118RFinalYearProject

AY2010/2011

Comparisonoffinitedifferenceandfinitevolume
methods
&the
Developmentofaneducationaltool
forthe
Fixedbedgasadsorptionproblem

FinalReport
by
YeoPuZhongOliver(U080118W)

CN4118RFinalReport

Abstract

Finite difference (FDM) and finite volume methods (FVM) are two important
categories of numerical methods used to solve partial differential equations. These
methodsareevaluatedandcomparedwithoneanother,withspecialemphasisputinto
thelatter.Thefixedbedgasadsorptionproblemisthemainfocusoftheprojectasitis
consideredasanimportanttechnologyinairpollutioncontrol.

Thisprojectinvestigatesexplicitschemesonfixedgrids,whicharebynomeans
advancetechniques.FVMwasfoundtobemoresuperiorinaccuracyandefficiencythan
FDM in most cases, especially those involving sharp discontinuities. It also conserves
mass rigorously, something the FDM cannot achieve. Various FVM schemes provide
advantages under specific conditions. For example, the Superbee scheme is especially
accurateandefficientwithsharpdiscontinuitiesbutfailwithsmoothsolutions.

The development of the educational tool in the form of a standalone simulator

usingMATLABisbrieflymentioned.
 CN4118RFinalReport

TableofContents

1. Introduction........................................................................................................................................1
PARTI:LITERATUREREVIEW
2. GasSeparationTechnologies......................................................................................................3
3. AdsorptioninGasSeparation.....................................................................................................5
4. NumericalMethodsforSolvingPartialDifferentialEquations....................................9
5. Theory.................................................................................................................................................10
PARTII:METHODOLOGY
6. NumericalExperiments...............................................................................................................33
PARTIII:RESULTS&DISCUSSION
7. Experiment1....................................................................................................................................37
8. Experiment2....................................................................................................................................44
9. Experiment3....................................................................................................................................49
10. Improvements.................................................................................................................................60
PARTIV:DEVELOPMENTOFGUI
11. GUI........................................................................................................................................................61
12. AlgorithmSpecifictoNonisothermalEquations.............................................................62
PARTV:CONCLUSION.................................................................................................................................63

NotationsandAbbreviations.................................................................................................................64
References........................................................................................................................................................65

CN4118RFinalReport

ListofFigures
Figure1Typesofadsorptionisotherms..............................................................................................12
Figure2TypeIisothermatdifferenttemperatures.......................................................................12
Figure3Shellandtubedepictionofvesselsection.......................................................................14
Figure4GriddefinitionforFDM.............................................................................................................21
Figure5GriddefinitionforFVM.............................................................................................................24
Figure6Godunovsschemeforfirstorderapproximation.........................................................25
Figure7Summaryoffirstorderschemes...........................................................................................26
Figure8FVMinterpretationwithlinearfunction...........................................................................26
Figure9Asummaryofdelimiterfunctions[35]..............................................................................29
Figure10FDMonplugflowmodelusingdifferentN....................................................................37
Figure11Dissipativeeffectoffirstorderupwindscheme.........................................................38
Figure12Effectsofgridsizeonfirstorderupwindscheme......................................................38
Figure13(a)SecondorderFVMschemes(b)HighresolutionFVMschemes...................40
Figure14ComparisonofLerrorsbetweendifferentschemes.................................................42
Figure15ComparisonofL2errorsbetweendifferentschemes...............................................43
Figure16Comparisonofmaxerrorsbetweendifferentschemes...........................................43
Figure17ComparisonofLerrorsbetweendifferentschemes.................................................44
Figure18Generaltrendobservedinloglogplotsforexperiment2......................................46
Figure19TrendobservedforSuperbeescheme.............................................................................47
Figure20ImprovementinFDMoverFVMasPedecreases........................................................48
Figure21Loglogplotsofmassbalanceerrorsagainstgridsizeforexperiment2.........48
Figure22Generaltrendobservedinloglogplotsforexperiment3......................................50
Figure23LoglogplotsofLerrorsagainstgridsizeforexperiment3................................50
Figure24Loglogplotsofmassbalanceerrorsagainstgridsizeforexperiment3.........51
Figure25Simulatedsolutionsforextremecases............................................................................53
Figure26ComparisonbetweensnapshotsofFDMandvanLeerschemes.........................54
Figure27Computationaltimeplotsfork=0.1................................................................................55
Figure28Computationaltimeplotsfork=1....................................................................................56
Figure29Computationaltimeplotsfork=10.................................................................................57
Figure30Computationtimebyequivalenterrorproduced.......................................................59

ListofTables
Table1Assumptionsforthederivationofgoverningequations..............................................13
Table2Componentsforgasandsolidphaseenergybalancewithrespectto(4)............19
Table3Componentsforcolumnwallenergybalancewithrespectto(4)...........................19
Table4Listofcoefficientsassociatedwith(27)and(28)...........................................................20
Table5NumericaldifferentiationforFDM........................................................................................22
Table6NumericaldifferentiationforFVM.........................................................................................29
Table7Functionsselectedforprofilingandtheirdescriptions................................................35
Table8SampledataforfirstorderupwindschemewithN=10.............................................39
Table9SampledataforFDMschemewithN=10..........................................................................42

 CN4118RFinalReport

1. Introduction

Thesimulationofaprocesshasbecomepartandparcelofaprocessengineers
workeversincethearrivalofmodernprocessorscapableofsolvingcomplexequations
inareasonableamountoftime.Theexpertiseinmodelingachemicalprocess,however,
began far earlier before that. Both endeavors arose from a need to predict a process
outcomeandimproveonexistingprocesses.

The importance of modeling and simulation became apparent in the

development of adsorption technology. Arguably, the underlying phenomenon in
adsorption has been understood fairly well; simple Langmuir models and more
complicated nonideal, multicomponent models have been developed to explain a
multitude of complex systems. However, as with all models, simplifying assumptions
havetobemade,whichbecameuntenableduetothedynamicnatureoftheadsorption
process. Unlike other major separation techniques, adsorption does not operate in a
steadystatemanner.Thispreventedthedevelopmentofadesignmethodaspredictive
as those used in distillation and absorption, impeding largescale adoption of the
technology[1].

Analytical solutions which work well with the aforementioned steadystate

technologieswouldnotbeadequateforadsorption.Numericalmethods,orsolvingthe
equationsbybruteforce,becameincreasinglypopularandfeasiblewithmorepowerful
computers. Simulators developed for specific systems such as the pressure swing
adsorption(PSA)requirecomplexboundaryconditions.Itallowsonetoquicklypredict
outcomesfromchangesinprocessconditions(e.g.feedcomposition,temperatureetc.)
andmakecorrespondingadjustments.Thecomputeralsoreducescapitalinvestmentin
thepreliminarydesignofanoveloperationcyclescheme.Therefore,acompetentmodel
put into an accurate simulator can allow a process engineer to optimize a process as
complexasthePSArapidlywhileinspiringindustrialconfidenceatthesametime.

This project encompassed two distinct phases. The firstphase involved solving
thenonlinearpartialdifferentialequations(PDEs)usingtwomethods:finitedifference
and finite volume methods. The performances of both codes were compared against
each other and the advantages of each method were highlighted. A simple fixed bed
adsorption model was selected for this phase: isothermal, isobaric, single component
andaxiallydispersedplugflowmodel.Theperformanceofeachsimulatorwasbasedon

1
CN4118RFinalReport

the errors produced relative to analytical solutions, their mass balance analysis and
computationaltime.Thefinitevolumecodeswerethenadaptedforthenextphase.

In the second phase of the project, a simulator was being developed as an

educational tool. A slightly more complex model was implemented, giving the user
optionsforadiabaticoperation.Thissimulator,intheformofastandaloneapplication,
incorporates a graphic user interface (GUI) developed using MATLAB. The purpose of
this application is to aid students in understanding the basic principles of adsorption.
One of its features helps the user visualize concentration and temperature waves
moving through the adsorber. Currently, commercially available modeling and
simulation software, such as COMSOL, can provide higher capabilities and robustness,
but their costs do not justify their use solely as an educational tool. An application
created using MATLAB allows one to distribute the software and the accompanying
proprietarycomponentsroyaltyfree.

2
 CN4118RFinalReport

PARTI:LITERATUREREVIEW

2. GasSeparationTechnologies

Separation processes are hugely important in a chemical plant and can be

responsiblefor40%to70%ofitscapitalandoperatingcosts[2].Apartfromtheirchief
functionofimprovingproductivitiesofunitoperations,separationequipmentisbeing
deployed as endofpipe solution to the problem of air pollution. Public power,
combustion plants and industrial combustions contribute to almost 70% of total CO2,
SO2andNOxemissionsinEurope,outstrippingtransportation(20.8%),thenextlargest
sourceofemissions[3].

Flue gas treatment as an endofpipe solution requires the removal of CO2, SO2
and NOx from postcombustion exhaust gas. The climate change conferences at
Copenhagen and Cancun had brought about growing awareness in combating climate
change and alluded to the possible global adoption of carbon tax and other economic
toolstoregulateemissions.Withthisinmind,researchintogasseparationmaybecome
vital,justastheairpollutioncontrolindustrybecomemoreprofitable.

Severalgasseparationtechnologiesexistandarealreadyindispensibleinmany
chemical processes. The sections below provide a brief review into three such
technologiescryogenicdistillation,absorptionandmembranetechnology.Adsorption,
thefocusofthispaper,willbepresentedinthefollowingchapter.

2.1 CryogenicDistillation

Based on the most common and widely used separation technique, cryogenic
distillationoperatesatlowtemperatures(below150C,typicalboilingpointsofgases).
It differentiates itself from other technologies by having no need of a physical
separating agent between different phases, thus eliminating the regeneration step. It
canproducehighpurityproductsbutisenergyintensiveandincurshighmaintenance
costs. Being a mature technology, the potential for cost cutting is limited. Thorough
understandingofchemicalthermodynamics[4],improveddesignofprocessequipment
and more complex process cycles [5] may reduce overall power consumption and
capital cost. Although it is unsuitable for flue gas treatment due to the enormous

3
CN4118RFinalReport

amountofenergyrequiredtocoolthehotfluegas,itremainsthepreferredtechnology
inproducingultrahighpurityoxygenandnitrogen.

2.2 Absorption

Gasseparationusingabsorptionisthephenomenonwheremasstransfertakes
place between the gas phase and a nonvolatile liquid phase. This equilibrium staged
process offers numerous advantages. The design of its unit operation is simple and
provides flexibility in continuous operation. Because of its cost effectiveness and
operability, gas absorption via chemical reaction, with monoethanolamine (MEA)
solution as the solvent, is widely employed in flue gas treatment. However, the high
capital cost of fresh solvent, substantial operating costs from its energy intensive
regeneration steps, as well as the large throughput of solvent involved, compares it
unfavorablytowhatcompetingtechnologiespromisetooffer[6].Severalinvestigators
are directing their efforts at finding more energy efficient ways of regenerating used
solvents[7],processmodifications[8,9]anddevelopmentofsuperiorsolvents[10,11].

2.3 MembraneTechnology

In 2000, the gas separation membrane module business was estimated at

US$125millionwithanannualgrowthrateof8%[12].Itisarelativelynewtechnology
withhugepotentialforgrowth.Itisslatedtobeinvolvedinwiderangingapplications,
fromhydrogenandhydrocarbonrecoverytoacidgastreatment.

Its advantages include simple modular design, ease of operation, long lifetime
andcontinuousoperationatnearambienttemperatures.Themaindifficultiesfacedin
its development are the conflict between membrane selectivity and permeability,
reliability issues caused by hydrocarbon condensation as well as high capital costs in
research and implementation. Currently, applications in large scale processes are still
uncommon[12].Mostmodulesareincorporatedwithamineabsorptionoradsorbents
toimproveexistingprocesses[13].

Despiteitsversatility,preliminarystudiessuggestthatitisapoorcandidatefor
fluegastreatmentwhichrequirehighpurityproduct[14].Thisisduetolimitationsin
producinghighpurityproductsineitherpermeateornonpermeatewithasinglestage;

4
 CN4118RFinalReport

a multistage design would render the operation uneconomical due to increased

operatingcostfromextracompressionsteps.

3. AdsorptioninGasSeparation

Principles,ApplicationsandFutureDevelopment

Adsorption,togetherwithchromatographyandionexchange,formauniqueset
of separation techniques that are characterized as heterogeneous (fluidsolid),
unsteadystate operations. Although the underlying phenomenon was first recorded
scientifically in the late 19th century, adsorption only had limited applications in air
purification up until the early 1960s [15]. Developments over the next 40 years have
establishedadsorptionastheleadingprocessintheremovaloftoxinsfromwastewater
andgas,withglobalsalesofadsorbentsvaluedatUS$1.15billionin1997[2].Thefuture
goalofthistechnologywouldbetodevelopcapabilityformoretypesofseparationsand
to compete economically with established techniques such as distillation and
absorption.

3.1 Principles

Adsorptionutilizesporoussolids asthemassseparatingagent,normallyinthe
form of a fixed packed bed; species are preferentially adsorbed onto the solid surface
while less strongly adsorbed species are allowed to pass through the column. The
species which are adsorbed are called adsorbates while the porous solids are called
adsorbents.Masstransferfromthefluidphasetothesolidphaseisusuallyequilibrium
controlled; steric and kinetic effects, two other adsorption mechanisms, are less
dominant.

The steric effect arises due to the uniformly shaped pores found in zeolite
molecular sieves which only allow certain shapes of molecules to diffuse into, and
subsequentlyadsorbedby,theadsorbent.Thekineticeffectreliesonthedifferencesin
diffusionratesofdifferentmolecules.Theequilibriumeffectdictatesthatmasstransfer
is driven by the difference between the fluid concentration and solid loading.
Equilibrium models are developed to describe the equilibrium relationships. The
Langmuir model and its isotherms (measurements are done at constant temperature)
arecommonandreadilyapplicabletoawiderangeofadsorbateadsorbentsystems.

5
CN4118RFinalReport

3.1.1 Adsorbents

Thestateofadsorptiontechnologyarguablydependsonthatofadsorbents.Gas
separation by adsorption has been advanced by the invention of synthetic zeolites in
1959 (by Milton) and the carbon molecular sieve (CMS) years later. These new
adsorbentscouldseparatespeciesbytheirmolecularshapeswhereasearlyadsorbents
suchasactivatedcarbonwereusedprimarilytopurifywastegasandwastewater[16].

Activatedcarbonishighlyporous,givingitalargesurfaceareatovolumeratio.
Its nonpolar, or slightly polar, surface allows it to adsorb organic compounds much
morestronglythanwatermolecules.Thispropertyeliminatestheneedtodrythegas
streampriortothefeedstep.Regenerationofactivatedcarbonrequireslesserenergy
ascomparedtootheradsorbentsduetothelowerheatofadsorptionthestrengthof
whichtheadsorbateisbondedtotheadsorbent.

CMSaremanufacturedfromthesamestartingmaterialasactivatedcarbon,but
additionaltreatmentwithhydrocarboncrackinggaveititsuniqueproperties.Itspore
structureistightlycontrolledandmoreuniformthanthatofactivatedcarbon.Kinetic
effect excludes larger molecules, while surface interactions remain similar to that of
activatedcarbon.CMSaresignificantlymoreexpensivethanactivatedcarbonbecause
oftheextramanufacturingstepandpatentprotection.

Zeolites are crystalline aluminosilicates whose pores have exact dimensions.

Thus zeolites can be designed for separation based on steric effect. However, the
equilibriumeffectismoredominantintheircurrentapplications.

Silicagelisanamorphoussolidmadeupofcolloidalsilica(SiO2).Itishydrophilic
innatureandhasalargesurfaceareatovolumeratio.Heatofadsorptionofwateron
thegelislow,whichhelpsduringregeneration(lowerregenerationtemperature)and
limitstemperaturechangesduringoperation.Itcanbedamagedbyliquidwater,alow
riskduringgasseparation.

Activated alumina, or dehydrated alumina hydrate, has excellent water

adsorption properties and is more durable than silica gel in the presence of liquid
water. However, prolonged use leads to coking and contamination, reducing its
adsorptioncapacity.

6
 CN4118RFinalReport

3.1.2 Regenerationmethods

Another mover of adsorption technology is the ability to regenerate the spent

adsorbentquicklyandusingaslittleenergy(ordesorbent)aspossible.Asadsorptionis
inherently an unsteadystate process, techniques are being developed for the sole
purposeofregeneration,allowingtheprocesstooperatecontinuously.Conceptssuchas
thermalswing,pressure/vacuumswingandsimulatedmovingbed(SMB)areusefulin
optimizing energy use and enhancing product purity. As SMB is currently limited to
liquid separations, only temperature swing adsorption (TSA), pressure swing and
vacuumswingadsorption(PSAandVSA)wouldbediscussedhere.

TSA is the oldest and most developed adsorption technique [15]. It works
because the adsorptive capacity of adsorbents decreases as temperature increases,
adsorbing more at low feed temperatures and desorbs readily at high regeneration
temperatures.However,thethermalcyclescancauseexcessivethermaldecomposition
at very high temperatures and the activated carbon bed can catch fire under certain
conditions (e.g. concentration above lower explosion limit, large heat of adsorption,
etc.)[17].

Pressureswingandvacuumswingadsorption(PSAandVSA)wasfirstdevelop
bySkarstromand,independently,byGuerindeMontgareuilandDomine[15].AllPSA
and VSA processes work on the same principle: high pressure adsorption and low
pressure regeneration. They are equilibrium controlled separation where different
amountofadsorbatecanbeadsorbedatdifferentpressuresthroughoutacycle.Most
PSA processes are operated adiabatically where pumps provide the compression or
suction energy required to achieve the separation [1]. PSA can achieve higher purity
separationthanmembranetechnologyandiseconomicalforlowtomoderatecapacity
units [1]. Activated carbon used with PSA is useful for separation of toxic gas species
fromfluegasesduetotheirhighcapacity,fastkineticsandrelativeeaseofregeneration
[18].APSAcyclecanbeveryfast;onetotwominutesiscommoninindustry.Thisleads
to high productivity and smaller adsorbers, reducing capital cost. A VSA cycle, on the
other hand, has long desorption step, thus reducing productivity. However, VSA can
producehigherpurityproductsincertainapplications[17].Thedesignandoperationof
PSAandVSAcyclescanbeveryversatile;thisflexibilityallowsforquickadaptationto
changesinfeedcompositionandmarketconditions[19].

7
CN4118RFinalReport

3.2 Applications

3.2.1 Byadsorbents

Activated carbon, being inexpensive to manufacture and a strong adsorbent of

organiccompounds,iswidelyusedinpollutioncontrolinvolvinghumidgases.Itisalso
thedefactoadsorbentusedinpersonalprotectionagainstatoxicatmosphere.Itisused
inH2purificationaswell.CMSareusedforproducingpurenitrogenfromair.

Zeolites are used for drying air and natural gas and the removal of carbon
dioxidefromairorfluegas.Theyarealsousedinsmallscaleairseparation,producing
oxygenandnitrogen.

Both silica gel and activated alumina are used in gas drying and are
complimentarytozeolites.Theyareusedtoreducewatercontentatlowertemperature
before the gas stream is passed through zeolites at higher temperature for complete
drying.

3.2.2 Byregenerationmethods

TSA is limited by slow heating and cooling rates. The long TSA cycles make it
unproductiveforbulkseparations.Itisemployedforpurificationprocessessuchasthe
removalofVOCs[20]anddryingofgases[21],wherethetimetakenforthefeedstepis
comparabletothatoftheregenerationstep.

OneofthebestestablishedprocessesinPSAisH2purification,withcapacitiesup
to 100 kN.m3/h and purities up to 99.9995%, with good energy efficiency and high
recovery (70 90%). O2 production with moderate purity (93 95%) and small
flowrates of up to 1 kN.m3/h using zeolites is common as well. Minor applications
include H2CO2, CO2CH4, CO from blast furnace gas, H2CH4, CO2N2 and ArO2
separations. VSA is commonly used for N2 production (up to 99.5% purity) at
moderateflowrate(~10kN.m3/h),usingCMSastheadsorbent.Pilarczyketal.reported
that new applications have been developed for recovery of CO2 from flue gases,
achievingCO2gasstreamsofhighpurity(>98%)andrecovery(5372%)[22].

8
 CN4118RFinalReport

3.3 FutureDevelopment

Theuseofsolidadsorbentsinfluegastreatmentisonlyintheearlydevelopment
stage but some studies have shown that they have superior thermal properties as
compared to liquid solvents used in absorption [23], important considerations with
regards to energy consumption. Nonetheless, many acknowledge that developing
adsorbents capable of meeting present demands in flue gas treatment remains an
ongoing process [19, 24, 25]. The cost of manufacturing adsorbent is now negligible
compared to the equipment cost, thereby providing incentive for the development of
moreexpensivebutefficientones[1].

In todays R&D in the technology, improving the energy demand of adsorption

process is one of its key goals. This has been achieved for O2 production via PSA; the
power consumption of these units has decreased from 0.6 kWh/N.m3 in 1983 to 0.35
kWh/N.m3in1990,bringingittothelevelofthebestcryogenicunits[1].Otherwaysto
bringdownoperationcostsaretodevelopadsorbentswithhigherselectivity,expansion
ofcapacitytoreapeconomyofscaleaswellasmorerefinedoperationschemes.Capital
costreductioncanbeachievedthroughmorecompactdesignofequipment.

4. NumericalMethodsforSolvingPartialDifferentialEquations

Understandingaprobleminadsorptionrequiressolvingitsgoverningequations,
which are usually partial differential equations (PDEs). These differential equations
differ from ordinary differential equations (ODEs) in that they have two or more
independent variables. Solving PDEs analytically, though not impossible, are tedious
and impractical. To that end, various numerical approaches are available to give
approximate solutions. The application of numerical methods to solve problems
involving fluid flow, such as adsorption, is known as computational fluid dynamics
(CFD).

Two methods currently dominate the CFD community: finite volume methods
(FVM)andfiniteelementmethods(FEM)[26].FVMareincreasinglybeingusedforCFD
inrecentyearsduetoparticularadvantagesinitsformulation[27].Anothermethod,
alsotheoldestandsimplest,isthefinitedifferencemethod(FDM).Theformulationof
these methods differs in grid definition and interpretation of numerical quantities.
Other methods such as orthogonal and spline collocation exist in literature and are

9
CN4118RFinalReport

knowntobemoreefficientandaccuratethantraditionalmethods[1,28,29].However,
theyarenotcommonlyappliedduetotheircomplexity.

TheapproachofFDM,FVMandFEMistoapproximatethePDEsbydiscretizing
the spatial dimension to convert the PDEs into ODEs in time. The ODEs can then be
solvedbypowerfulODEsmethodssuchastheRungeKuttamethod.

TraditionalFDMcodescannothandlediscontinuities,abaneinCFDwheresharp
fronts can occur [30, 31]. Another issue with FDM is that mass is not rigorously
conserved. Although accuracy can be improved as time step and grid size approach
zero and by using higher order of approximation, computation time increases as well
[30, 32]. FVM codes rely on integration to give a more fundamental interpretation of
mass(orenergy)fluxesbetweeninterfaces[33];therefore,conservationofmass(and
energy)isrigorous.FVMcanalsobeeasilyusedonnonuniformgrids[34,35].Several
papers comparing FDM and FVM were published. FVM codes were found to be more
accurate and as a result, coarser grids can be used, saving computation time [36, 37].
Botte al et. concluded in their analysis that although FVM conserve mass and is more
accurate in some cases, generally speaking, FDM codes are more accurate and
computation times for either codes are comparable [32]. The FDM can also be more
efficientthanFVMorFEMiftheproblemcanbesolvedwithauniformstructuredgrid
[38].

Inthispaper,asimpleadsorptionmodelwillbesolvedusingFVMandFDMusing
fixed grids (adaptive grids, or stencils, allow superior tracking of sharp fronts but is
beyond the scope of this project). FEM was not performed due to constraints in the
project timeline; its complex nature would require significant effort in formulating its
codes. Three numerical experiments were performed to evaluate the performances of
FVMandFDMcodes.

5. Theory

5.1 LangmuirModel

The Langmuir model proposed by Langmuir in 1918 [39] is used to describe

adsorbateadsorbentsystemsbasedonthekineticapproach.Therelationshipbetween
theamountofadsorbateadsorbed(q)atequilibriumwasexperimentallydeterminedto
be dependent on concentration and temperature. Equilibrium isotherms are plots

10
 CN4118RFinalReport

where temperature is fixed and q is dependent on concentration only. Six types of

isothermsareknown(Figure1);TypeI,alsocalledthefavorableisotherm,isofinterest
hereandcanbedescribedbytheLangmuirmodel.

5.1.1 Monolayeradsorptiononhomogeneoussurfaces

ThreeassumptionsaremadebasedonLangmuirsanalysis[15]:

1. Theadsorbedmoleculeisheldatdefinite,localizedsites.

2. Eachsitecanonlyadsorbonemoleculeoratom.

3. Theheatofadsorption,H,isconstantthroughoutthesurfaceandnointeraction
existsbetweenneighboringadsorbedmolecules.

The isotherm equation is derived based on equal rates of adsorption and

desorption at equilibrium. The rate of adsorption is dependent on the fluid phase
concentrationoftheadsorbate(C)andthefractionofremainingsites(1 ).Therate
ofdesorptionisonlydependentonthefractionofoccupiedsites( ).

Rateofadsorption ka C 1
Rateofdesorption kd

Atequilibrium,ka C 1 kd ,or

q bC ka
, b (1)
qs 1 bC kd

b is the Langmuir constant and is found to follow an Arrhenius type dependence on

temperature.

b b0 HRT
(2)

b0 canbedependentorindependentontemperaturefordifferentsystems.

LinearisothermsareusedincaseswhereHenryslawisvalid.Atlow
concentration(i.e.bC 1),(1)reducestoalinearform(3)whereKistheHenrys
constant.Thecorrespondingisothermiscalledalinearisotherm.Itssimplicityisalso
usefulforarrivingatananalyticalsolutionusedforvalidatingthesimulator.

q KC, K qs b (3)

11
CN4118RFinalReport

Figure1Typesofadsorption isotherms
Ordinateamountadsorbed,q;Abscissarelativepressure,P/P0.Source:[40]
5.1.2 Effectoftemperature

Theisothermsarecompresseddownwardsastemperatureincreases(Figure2).
This is because adsorption is an exothermic reaction. By Le Chateliers principle, the
equilibrium shifts towards higher rate of desorption as temperature increases. This is
thebasisforthermalregenerationcycles,wherethetemperaturewithintheadsorberis
manipulatedforfasterdesorptionrateorhigheradsorptioncapacity.

T1
Amountadsorbed

T2

T1<T2

AdsorbateConcentration

Figure2TypeIisothermatdifferenttemperatures

12
 CN4118RFinalReport

Table1Assumptionsforthederivationofgoverningequations.
Assumption Comments

1. Homogeneouspacking(nochanneling) ValidifcarefullypackedandDcol/dp>
30

2. Negligibleradialgradients

3. Nochemicalreactions Adsorbentscanactascatalyst

4. Neglectkineticandpotentialenergy

5. Noradiantheattransfer Absolutelytrueonlyforisothermal;
usuallylumpedwithconvectiveheat
transferfornonisothermalcase

6. Noelectricalormagneticfields

7. Nophasechangesexceptsorption Validforgassystem

8. Velocityconstantacrosscrosssection Reasonableassumptionif:no
channelling,noradialgradients,no
viscousfingering

9. Noirreversibleadsorption

10. Rigidpacking

11. Constantsorptionproperties Nodegradationorpoisoning

12. Nobreakageordissolutionofpacking

13. Negligibleconductionincolumnwalls Validforlargediametercolumns

14. Constantadsorbatemoleculardiffusivity Validforisobaricandisothermalsystem;

maybeapplicabletoadiabaticsystemif
temperaturechangeisinsignificant.
Adaptedfrom[17],Table176,pg673.

5.2 AxiallydispersedPlugFlowModelwithNegligiblePressureDrop

The axiallydispersed plug flow model was probably first described in 1949 by
Glueckauf et al. [41] and is widely used in the study of fluid mixing in chemical
processingvesselsbecauseofitssimplicity[42].Itischaracterizedbythespreadingofa
pulse as it progress through the vessel. This phenomenon is analogous to a diffusion
processofthetracerinbothdirectionsfromthecentreofthetracer,superimposedonto
theplugflowprocessofthesamecentremovinginonedirection[35].Forasufficiently
longvessel,radialdispersionisassumedtobenegligibleandthereisnovariationinthe

13
CN4118RFinalReport

radialdirection.Alongitudinaldispersioncoefficient,DL1,canbeusedtorepresentthis
spreadingprocess,reducingittoaonedimensionalproblem.

In deriving the governing equations, several assumptions have to be made for

theanalysistobevalid.TheseassumptionsaresummarizedinTable1.

5.2.1 Isothermal

InadditiontotheassumptionsconsideredinTable1,temperatureisassumedto
beconstantfortheentiresystem.Consideranarbitrarysection(z)ofapackedvessel
(withuniformcrosssectionalarea,S)asshowninFigure3,thegeneralmaterialbalance
forasinglecomponentAwithinthiscontrolvolumeatanarbitrarytimeinterval(t)is

Accumulation=InputOutput Source (4)

Figure3Shellandtubedepictionofvesselsection

LetCA andudenote the fluid phase concentration of A and the fluid velocity in the z
direction respectively; subscripts z andz+zrepresent spatial position of respective
quantities. denotesthevoidfraction.Theindividualcomponentsin(4)are

Accumulation CA Volumeoffluidinsection CA z S
Inputduetobulkconvection CA,z Volumeentering CA,z uz t S
Outputduetobulkconvection CA,z+z Volumeleaving CA,z+z uz+z t S

1Thelongitudinaldispersioncoefficientcanbeobtainedviasuitableempiricalcorrelations.

14
 CN4118RFinalReport

Fickslawonmoleculardiffusioncanbeappliedforthedispersioneffect,withDL
CA
assumedtobeconstantthroughoutthelengthofthevesseland isthelocal
z
concentrationgradientofcomponentA,
CA
Inputduetoaxialdispersion NA,z DL S t
z z

CA
Outputduetoaxialdispersion NA,z+z DL S t
z z+z

Leavingthefinalterm(Source)asitisforthewhileanddividingbyz S t,

CA CA
CA DL DL CA,z+z uz+z CA,z uz 1
z z+z z z
Source
t z z z S t

Takingthelimitsofz0andt0,

CA CA CA u 1 (5)
DL [Source]
t z z

(5) represents the general equation for the isothermal axiallydispersed plug flow
model.Thefinalterm,[Source],isrelatedtothecontrolvolumez Sandtimeinterval
t; hence these terms are omitted for later sections. It is positive if species A is
generated in the control volume (desorption) and negative if it is being consumed
(adsorption).

5.2.2 ExtensiontoadsorptionusingLangmuirmodel

In time interval t , solid loading q increases by qA while fluid phase

concentrationdecreasebythesameamount.qA isdefinedtohavethesameunitsasCA .
Asitdescribesquantitiesdependentonthesolidphasebutisaffectingthefluidphase,
theexpressionfor[Source]is

qA qA
Source Volumeofpackedsolid 1 z S
z S t z S t

Substitutingbackinto(5),

CA CA CA u 1 qA (6)
DL
t z z t

15
CN4118RFinalReport

Evidently,thefinaltermin(6)describestherateofmasstransferbetweenthefluidand
solid phases. A detail approach to quantifying the rate involves obtaining parameters
for film diffusion and diffusion within the pores of the particles. A lumped parameter
mass transfer coefficient, k,is often usedfor simplification and practical reasons [17].
The lumped parameter model assumed a linear driving force (LDF) caused by
concentration differences between the solid and fluid phases. One approach converts
the fluid phase concentrationCA into a hypothetical solid loadingqA , the amount
adsorbedthatwouldbeinequilibriumwithCA .

qA (7)
k qA qA
t

Applying(1)forqA ,

bCA
qA qs (8)
1 bCA

OverallmaterialbalanceisderivedforaninertcarrierbydenotingCI asthefluidphase
concentration of the inert carrier in (6) and summing the two material balances
together.

CI CI CI u
DL (6)
t z z
Note:Nomasstransferofinertcarrierintosolidphase

CA CI CT andCT isinvariantintimeandspace(isobaricandisothermal),

u 1 qA (9)
CT
z t

(6), (7) and (9) provide the governing equations for the isothermal model and can be
expressed in their dimensionless forms (10), (11) and (12) by defining the following
variables:

z CA u L u tu qA kL qs 1
Z ,y , ,u , ,x , ,
L A CT DL u L A qs u CT

whereListhelengthofthevessel,CT istheoverallfluidconcentrationandu isthefeed

interstitialvelocity.(13)and(14)providetheboundaryconditionsforadsorptionand
desorptionrespectively.

16
 CN4118RFinalReport

yA 1 yA yA u xA (10)

Z Z
xA qA byA CT
xA xA , xA (11)
qs 1 byA CT

u xA (12)

z
yA xA yA
0, 0, yA yA0 , 1 (13)
Z
Z Z

yA xA yA
0, 0, yA 0, 1 (14)
Z
Z Z

Byidealgaslaw,

P Q
CT ,whileu
V RT S

where P, T and Q are the operating pressure, temperature and feed flow rate
respectively. R is the molar gas constant. The quantitiesyA andxA can be normalized
through dividing by feed mole fraction (yA0 ) and the equilibrium fractional coverage
associated withyA0 (xA0 ) respectively. In subsequent sections,yA0 is defined as one of
theoperatingparameterswhilexA0 isobtainedfrom(15).

byA0 CT
xA (15)
1 byA0 CT

5.2.3 Nonisothermal

The nonisothermal model contains two more governing equations: one heat
balanceforthefluidandsolidphasesandoneheatbalanceforthecolumnwall.Before
goingintotheheatbalances,thematerialbalanceshavetobemadeintomorecomplete
forms.(6)remainsvalidandshouldbeexpandedbyidealgaslaw.CA yA PRTandR
andPareinvariantwithtimeandspace,(6)becomes

yA yA T yA 1 yA T u yA u T RT 1 qA
DL (16)
t T t z T z z z T z P t

17
CN4118RFinalReport

Overallmaterialbalanceisobtaineddifferently,withCT afunctionofspaceandtime.

CT CT u 1 qA (17)

t z t

Uponexpansionbyidealgaslaw,(17)becomes

1 T u u T RT 1 qA (18)

T t z T z P t

Substituting(18)into(16),

yA yA 1 yA T u u RT 1 qA
(19)
DL yA y 1
t z T z z z z P A t

Energybalanceforthefluidandsolidphasesisderivedusingthesameprocedureasin
section5.2.1,withthecomponentsin(4)listedinTable2.

1 T 1 T CT T
s Cps qA Cpa Cpg
t t t
(20)
Kz T CT T u 1 qA 2hi
Cpg H T Tw
z z t ri

Table3providesthereferenceforenergybalanceinthecolumnwall.
Tw Tw 2ri hi 2ro ho (21)
w Cpw Kw T Tw T Ta
t z ro ri ro ri w

(20)canberearrangedforfurthermanipulations

1 T Kz T CT CT
g Cpg s Cps qA Cpa Cpg T u
t z t z
Tu 1 qA 2hi (22)
Cpg g H T Tw
z t ri
Note:CT g

Byonlyexpandingthefirsttermontherighthandside,(17)canberewrittenas

CT CT u 1 qA (23)
u g
t z z t

18
 CN4118RFinalReport

Substituting(23)intothesecondtermontherighthandsideof(22)andrearranging,
thefinalformofthegasandsolidphaseenergybalanceisgivenas

1 T Kz T u Tu
g Cpg s Cps qA Cpa g Cpg T
t z z z
(24)
1 qA 2hi
Cpg T H T Tw
t ri

(18),(19),(21),(24)arethematerialandheatbalancesofthemodelandtogetherwith
(7),theLDFmodel,constitutethegoverningequations.

Table2Componentsforgasandsolidphaseenergybalancewithrespectto(4)
Terms Components

Accumulation Energygainedinthefluidphase
Energygainedinthesolidphase
Energygainedintheadsorbedphase

Input Energytransferredintocontrolvolumebyconvectionanddiffusion

Output Energytransferredoutofcontrolvolumebyconvectionanddiffusion
Energytransferredintocolumnwallbyconvection

Source Energygeneratedbyadsorption(positive)orconsumedby
desorption(negative);storedasadsorbateinternalenergy

Table3Componentsforcolumnwallenergybalancewithrespectto(4)
Terms Components

Accumulation Energygainedbywallsection

Input Energytransferredintowallsectionbyconduction
Energytransferredintowallsectionfromfluidbyconvection

Output Energytransferredoutofwallsectionbyconduction
Energytransferredoutofwallsectionintoambientairbyconvection

Source Notapplicable

As with section 5.2.2, dimensionless forms of (18), (19), (21) and (24), together with
(11),areusedinthesimulator.Definingadditionalvariables,

T Tw Ta 1
T , Tw , Ta , g Cpg s Cps qA Cpa
T T T

19
CN4118RFinalReport

u 1 T u T xA (25)
T
Z T T Z
yA 1 yA 1 yA T u u xA
yA T yA 1 (26)
Z T Z Z Z Z

Tw Tw (27)
1 2 T Tw 3 Tw Ta
Z
T T u Tu xA
4 5 T 6 7 T 8 T Tw (28)
Z Z Z

Theinflowandoutflowboundaryconditionsare(29)and(30).

AtZ 0,
1 yA

Z Z Z
(29)
T T
4 5
Z Z Z
Tw Ta , 1
AtZ 1,
yA (30)
T Tw
0
Z Z Z

Afulllistofcoefficientsusedin(27)and(28)isgiveninTable4.

Table4Listofcoefficientsassociatedwith(27)and(28)
Coefficient Coefficient
Kw g Cpg
1 5
w Cpw u L
2ri hi L 1
H qs
2 6
2
w Cpw u ro r2i T
2ro ho L 1
3 7 Cpg qs
w Cpw u r2o r2i
Kz 2hi L
4 8
u L ri u

5.3 FiniteDifferenceMethod(FDM)

Inthisandthefollowingsections,generalformulationofthenumericalmethods
will be presented. will be used to denote the variable of interest. Nodes are

20
 CN4118RFinalReport

discretized locations where are defined. The quantities obtained via numerical
methods are not continuous; that is to say, a set ofN 1 values form the entire
solutionofthePDEs.Inthispaper,onlyfixedgridswillbeused.

ExplicitschemesareusedfortheFDMandFVM.Thismeansthatthevaluesof
obtained for the next time step is only dependent on values at the current time step.
Explicit schemes are more efficient and can give reasonable accuracy for practical
purposes.(10)isreferredtoagainintermsof ,therefore,

1
u Source'
Z Z
u
Note:Source' Source
Z

TheconvectiontermwasexpandedandcombinedwiththeSourcetermtosimplifythe
expression. This can only apply for the FDM since the fluid velocityudoes not have
physicalsignificanceinFDM;thisisnottrueforFVM.

5.3.1 Griddefinition

Figure4depictsthegraphicalrepresentationofanadsorbercolumnoflengthL
dividedintoN cellswithN 1nodes,withthefirstnodeatZ 0andthelastnodeat
Z 1.EachnodeisseparatedfromtheotherbythegridsizeZ.Boundaryconditions
areimposedatthefirstandlastnodes.Forsubsequenttreatmentof ,subscriptsirefer
tothenodelocationandsuperscriptsnrefertothecurrenttimestep.

Figure4 GriddefinitionforFDM
5.3.2 Spatialderivatives

Numerical differentiation is used to approximate the spatial derivatives as

algebraic expressions. Forward differencing is used for the first grid while backward
differencing is used for the last grid. Centered differencing is used for the grids in
between.Secondorderaccurate( Z )formulationswereusedasthestandardhere.
TheformulaearepresentedinTable5.

21
CN4118RFinalReport

Table5NumericaldifferentiationforFDM
Firstderivative Secondderivative
Nodes

Z Z
3 4 2 5 4
1
2Z Z

2
i
2Z Z

4 N 3 N 1 N 1 2 N 1 5 N 4 N 1 N 2
N 1
2Z Z
5.3.3 Remainingtermsandequations

Thesourcetermfollowsthealgebraicexpressionin(11).Overallmassbalance
(12) is implemented with backward difference using the inflow boundary condition
1.

5.3.4 FDMalgorithm

The fluid phase concentration, solid phase loading and fluid phase velocity are
updated using (31), (32) and (33) respectively. Boundary conditions for adsorption
operation prescribed by (34) (yA 0for desorption) and the source term (35) are

performedatthestartofeachtimestep.

1 yA yA xA
yA yA u yA 1 , 2N (31)
Z Z

xA (32)
xA xA , 1 N 1

xA (33)
u u Z , 1, 2 N 1

ZyA yA (34)
yA , yA yA0 , yA yA
Z 1
xA byA CT
xA * xA , xA * , 1N 1 (35)
1 byA CT

22
 CN4118RFinalReport

5.4 FiniteVolumeMethod(FVM)

InFVM,thedifferenceslieinthegriddefinitionandtheinterpretationofhowthe
quantities should be evaluated. Numerous schemes are available in the literature,
categorizedintoexplicitimplicitandlinearquadraticschemes.Asmentionedinsection
5.3,thescopeofthispaperislimitedtoexplicitschemesandfixedgridsonly.

Anotherformof(10)isderived.

1 (36)
u Source
Z Z
The terms in the parenthesis has significance in determining the stability limits of
explicit FVM algorithms. In physical terms, the Peclet number determines whether
diffusion or convection is more dominant. (36) can be expressed in a more general
form:

Source, (37)
Z Z

Solutions derived based on , which has coefficient , can be extended towith

coefficient.

5.4.1 Griddefinition

Figure5depictsthegraphicalrepresentationofanadsorbercolumnoflengthL
dividedintoNcellswithNnodes.Eachnodeislocatedatthecenterofitscelland is
definedastheaveragevalueofthequantitywithinthecell.Inflowboundaryconditionis
imposed on the left interface of the first cell and the outflow boundary condition is
imposed on the right interface of the last cell. In addition to the notations introduced
earlierinsection5.3.1,subscriptslandrrepresentsintermediatevaluesof attheleft
andrightinterfacesrespectively.Theseintermediatevaluesarenotretainedinthefinal
solution, with the exception of the extreme left and right boundaries, denoted by
subscriptsLandRrespectively.

Another difference with FDM is the introduction of ghost cells at either

boundary.Twoghostcellsateachboundaryarerequiredfortheimplementationofhigh
resolutionschemes.

23
CN4118RFinalReport

Figure5 GriddefinitionforFVM
5.4.2 FVMalgorithm

The FVM algorithm can be described by the 3steps Godunovs scheme. The
algorithm first define the interface values, reconstructing discrete solution into a
piecewisecontinuousone.Theflowenteringandleavingthecontrolvolume(CV)ina
singletimestepisthenapproximatedbyintegratingthecurrentsolution atthecell
overthedistancecoveredbythefluid.Thenetflowisaddedtotheinitialintegraland
finally,theintegralisdividedbythecellvolumetoobtainthenewsolution .Figure
6illustratesthealgorithmforasingletimestepforafirstorderscheme.

Thisinterpretationcloselyresemblesthephysicalimplicationoftheequation.It
can be thought of as the deconstruction of the equation, tracing back the steps of
derivingitinthefirstplace(section5.2.1).OneadvantageoftheFVMcanbeobserved
here.Thefluxattheleftboundaryoftheithcellisequaltothefluxattherightboundary
ofthepreviouscell.Thiscontinuityensuresthatthequantityisconservedrigorously.

5.4.3 Firstorderschemes

Severalmethodsareavailabletoapproximatetheintermediatefluxterms , ,
and , , . In firstorder schemes, the cell centered average value is presumed to
apply across the entire cell (Figure 5(a)). As such, can be used directly to
approximate values at the interface. The simplest of all is the firstorder upwind
scheme.Inthisscheme,thedirectionoftheflowistakenintoaccountandinflowfluxes
are only dependent on the information propagating from upstream. Using the
convection term once again to illustrate, for a lefttoright flow, the upwind scheme
essentiallyequates to , and to , andsoforth.

24
 CN4118RFinalReport

Discretized Algebraic
Description
Expression Expression
Step1

Defineapiecewisecontinuoussolution

Z, Z
FlowintoCV Fluxatleftboundary
CrosssectionalArea , , S
Time
FlowoutCV Fluxatrightboundary
CrosssectionalArea , , S
Time
Step2

, ,
NetFlowoutofCV V
, , S Z

VolumeintegralofquantityinCVat
S Z

S Z
Volumeintegral NetFlowoutofCV , ,

, ,
VolumeintegralofquantityinCVat
S Z
1
Step3

1
Timederivativeofquantity
Z , , Z

, ,
Figure6Godunovsscheme forfirstorderapproximation

Another approach is to use central differencing to approximate values at the

interfacetobehalfthesumofcellcenteredvaluesonbothsides.Thismethodwillgive
rise to unstable flux and produce oscillations in the solution as a result. The Lax
Friedrichsmethodmodifiesthecentraldifferencingwithanadditionaltermwhichacts
as numerical diffusion. A summary of firstorder schemes making use of the result
obtainedinFigure6ispresentedinFigure7.

25
CN4118RFinalReport

Upwind
Z Z

CentralDifferencing

(Unstableflux) Z 2Z

Z 2Z
LaxFriedrichs
2

2
Figure7Summaryoffirstorderschemes

These firstorder schemes are similar to the FDM numerical differentiation; in

fact, the LaxFriedrichs scheme is similar to that presented Table 5 for the central
differencing of the first spatial derivative. However, firstorder schemes in FVM are
problematic. The upwind scheme is dissipative and will smooth out any sharp front
while the other two central differencing schemes can produce oscillations, if they are
stableinthefirstplace.Higherorder(orhighresolution)schemesderivedwithslightly
moreeffortcanbevastlysuperiortothesesimplefirstorderones.

5.4.4 Higherorderschemes

Figure8FVMinterpretationwithlinearfunction

Higherorderschemesapproximatethesolutionwithinacellasafunction.The
simplestway,whichissecondorderaccurate,istousealinearfunctionthatintersects
theaveragevalue atthenodeslocation(Figure5(b)).Eachcellisassignedavaluefor
theslope( )ofitslinearfunction.TheresultobtainedinFigure6nolongerappliesbut

26
 CN4118RFinalReport

thealgorithmremainsidentical.Forsimplicityandpracticalpurpose,derivationofthis
resultisperformedonthespecificassumptionoflefttorightflow.

Assumingagainthatinformationispropagatedonlyfromtheleftinterface.Ina
singletimestep ,theshadedareainFigure8representstheamountflowintotheCV
incelli. Zisaproductof and , .TheflowintotheCVisequaltotheshadedarea
multiplybythecrosssectionalarea,S.Theshadedareacanbefoundbysomegeometric
manipulation

Z , Z
shadedarea
2
Z Z
where , and Z Z .Substituting,

Z
FlowintoCV 1 S ,
2 Z

FlowoutoftheCVisderivedsimilarly.

Z
FlowoutofCV 1 S ,
2 Z

, ,
NetflowoutofCV S Z
1 , 1 ,
2 Z Z

VolumeintegralremainsthesameasinFigure6becausethesymmetricfunctiondoes
notaltertheareaboundedbeneathit.

Volumeintegralofquantity at S Z

Summing up the change in quantity in cell i and dividing by cell volume (S Z) to

obtainthenewcellcenteredaveragevalue, .

, , 1 , ,
1 , 1 ,
Z 2 Z Z

Rearrangingthetermstoobtainamorerecognizableform,

, , , , (38)

Z

Z , Z ,
where , 1 and , 1
Z Z

27
CN4118RFinalReport

Note that the coefficient , and , are not defined by their cell centered values .
This is an intentional simplification by assuming that the value of the coefficient does
not changesignificantly across the cell and thus can be represented by the first order
upwind method. This assumption would be invalid in concentrated system where the
velocity,representedbythecoefficient,canvarydrasticallyacrossacell.Theremedyis
tousefinergridswhensolvingsuchsystems.

The manner inwhich the slopesarerepresented differentiates various second

orderschemes.Threeobviousnumericalapproximationsoftheslopescometomind.2
pointsbackward(upwind),forward(downwind)andcenteredapproximationsproduce
theBeamWarming[43],LaxWendroff[44]andFrommschemesrespectively.

However, the linear secondorder schemes can produce oscillations near sharp
wave fronts. Their approximations to the slopes are static, in terms of their formulae,
globally. A dynamic approach compares the three approximations and selects the
appropriate one locally. This is the working principle of high resolutions schemes. In
high resolution schemes, can be the product of a delimiter function r and ,
where

andr

isinsertedtopreventdivisionbyzeroduringcomputationandshouldbeaverysmall
number (~1010). Most first and secondorder schemes introduced earlier can be
expressed in this form by changing the delimiter function. A summary of wellknown
schemesispresentedinFigure9.Adetailedexplanationonhowthedelimiterfunctions
workcanbefoundin[35].

5.4.5 Remainingtermsandequations

Forthediffusionterm,itiscommonpracticeinFVMtoadoptcentraldifferencing
approximation for the first spatial derivatives. For the first cell, forward difference is
used. Boundary condition is imposed on the last cell. Table 6 presents a summary of
numericaldifferentiationformulaeforFVM.

Asecondorderaccurateschemeinvolvingthediffusionterm,orsecondspatial
derivatives in general, was developed by John Crank and Phyllis Nicolson [45]. The

28
 CN4118RFinalReport

CrankNicolson scheme is an implicit method and therefore is not desirable for the
purposeofthispaper.

Upwind r 0 (a)

Beamwarming r r (b)

LaxWendroff r 1 (c)

1
Fromm r r 1 (d)
2

Superbee r max 0, min 1,2r , min 2, r (e)

r 1
MC r max 0, min , 2, 2r (f)
2

r |r |
vanLeer r (g)
1 |r |
Figure9Asummaryofdelimiterfunctions[35]

Table6NumericaldifferentiationforFVM
Firstderivative
Nodes
,
Z ,

2
1
Z

2N
Z
N 1 0
Note: , ,

(38)isderivedforpureconvection.(39)appliesforaconvectiondiffusionplus
sourceform.

, , , , , , , ,
Source (39)
Z Z

Source terms are implemented in the same manner as in FDM (35) using cell
centeredvalues(see(40)and(41)).Higheraccuracymethodsforupdatingthesource
terminvolvesolvingtheconvectiondiffusionequationandsourceequationatseparate

29
CN4118RFinalReport

timesteps.Theseareknownasfractionalstepmethods[35].Despitetheirusefulness,
thesemethodsarenotpracticaltothisprojectforreasonsmentionedinsection6.1.

Mass balance is performed by FDM and using FDMdefined nodes (42). This is
becausein(38)onlycellinterfacevaluesareused.Intheeventthatcellcenteredvalues
arerequired,linearinterpolationshouldbeperformed.

xA (40)
xA xA , 1 N

xA byA CT
xA * xA , xA * , 1N (41)
1 byA CT

xA (42)
u, u , Z , , 1, 2 N 1

5.4.6 Boundaryconditions

(43)givestheboundaryconditionforadsorption(yA 0fordesorption).

ZyA 2yA
yA yA yA , yA yA0 ,
Z 2 (43)
yA yA yA yA

5.4.7 Effectsof

ThesizeofthetimestepcanaffectthequalityofthesolutionobtainedinFVM.
Followingthealgorithmderivedinthissection,astrictupperlimitisimposedon ;it
cannotbelargerthanthetimetakenforinformationtopropagateacrosstwocells.The
solutioninaparticularcellisassumedtobedependentoninformationfromitsadjacent
cellonly.Asimpleinequality(44)embodiesthislimitation.

1 (44)
Z

Ifthecoefficient representsvelocity,theinequalitypreciselymeansthat cannotbe

largerthanthetimetakenforthewavetotravelmorethanthedistanceofonecell(Z).
The left hand side of the inequality in (44) is known as the Courant number (Co), an
importantparameterinCFD.

Unsurprisingly, the Courant number appeared in the high resolution equation

(38). The high resolution equation reduces to a firstorder upwind scheme if the

30
 CN4118RFinalReport

Z
Courantnumberisequalto1,orif ,regardlessofthetypeofdelimiterfunction

used.Thisrepresentsthescenariowherethesolutionmovesonecelltotherightwith
no changes to the cell centered average values at all (i.e. ). It will provide
the exact solution for constant . For problems with variable coefficients, a smaller
Courantnumberisusuallyusedtotakeadvantageofthedelimiterfunction.

5.4.8 EffectsofZ

Refiningthegridexcessively(Z 0)willreducetheupperlimitof in(44),

pushingthesolutionclosertoinstability.Iftheproblemisdominatedbydiffusion(i.e.
small Pe), an excessively fine gridwill also cause instabilityin thediffusion term as it
scalesaccordingto1 Z .Investigatingtheerrorsinrelationtothegridsizeisoneof
thepurposesinthenumericalexperiments.

5.5 Analyticalsolutions

5.5.1 Axiallydispersedplugflowwithnosourceterm

TheanalyticalsolutionforansemiinfinitemediumwasdevelopedbyOgataand
Banks[46]

yA 1
, (45)
yA0 2
2 2

fortheboundaryconditions

yA 0, yA0 andyA , 0

5.5.2 Axiallydispersedplugflowwithadsorption

Klinkenberg provided the analytical solution based on semiinfinite medium,

dilutefeedandlinearisotherm[47].Theboundaryconditionsaresimilartothatinthe
precedingsection.

yA 1 1 1
z,t
yA0 2 8 8
(46)
kKz 1 z
where and k t
u u

31
CN4118RFinalReport

5.5.3 Massbalance

The mass balance can be performed across the whole adsorber. The input and
outputinformationisintegratedovertimewhiletheinformationinsidetheadsorberis
integratedoveritslength.

AmountofcomponentA,

Fedatthestartofthecolumn yA0 CT u

Obtainedattheendofthecolumn yA,out CT u

L
Remaininginsidethecolumn yA CT

L
Adsorbed xA qs 1

Assuming that the adsorbents are completely used up,yA yA0 andxA xA0 . The
simulatorcanprovidetheoutputinformation.Puttingthetimeintegralonthelefthand
sideandtheothertwotermsontheright,

L 1 qs xA0
1 yA,out u 1 (47)
u CT yA0

32
 CN4118RFinalReport

PARTII:METHODOLOGY

6. NumericalExperiments

A major aim of the numerical experiments is to compare the performance

between FDM and FVM as formulated in Chapter 5. The series of experiments were
designed to begin from the most simple equation pure convection, through
progressivelycomplicatedonesbyaddingthediffusiontermsandthesourceterms.

6.1 FVMImplementationinMATLAB

Thissectionshalldescribetheimplementationof(39)to(43),withcellinterface
valuesapproximatedby(38)andTable6.Firstorder,secondorderandhighresolution
methods were tested in Experiment 1 (section 6.3.1) while only high resolution
methods were used in the remaining experiments. The software used was MATLAB
(version7.10.0.499R2010a).

MATLAB provide two ways to solve PDEs. Although there is a builtin PDE
toolbox,theFDMandFVMalgorithmsderivedinChapter5shouldbesolvedusingODE
methods. A builtin MATLAB ode solver (ode45), which uses the RungeKutta method,
canbeusedtosolvenonstiffdifferentialequations[48].

The solver automatically adjusts the time step to produce accurate results. It
acceptsformulationoftheODEsinafunctionhandlemfile2ofthefollowingform:

, iscomputedatthestartofeachsolverloopinthefunctionhandlemfile
before being passed into the main solver mfile, where is determined by
perturbation,afterwhichthenewsolutionisobtainedfromtheoldone.Itisapparent
that the high resolution scheme described by (38) cannot be used directly because
cannotbecalledintothefunctionhandlemfilebeforeitisdetermined.Thisisalsothe
reason why fractional time step methods involving the source term is difficult to
implement using a variable time step solver like ode45. To solve this problem, the
Courantnumberin(38)isassumedtobemuchsmallerthanone.Thevaluesatthecell
interfacesarethereforeapproximatedas:

2FunctionsusedinMATLABarestoredasmfiles.

33
CN4118RFinalReport

Z Z
, r and , r
2 2

The validity of this approximation is dependent on the solvers ability to use a

timestepthatsatisfiestheassumptionmade.

6.2 PerformanceIndices

Threeperformanceindicesusedintheexperimentsaredefinedhere.

6.2.1 Accuracyerrors

The solutions obtained via the two methods are compared to a benchmark
solution. The benchmark solution are provided by the analytical solution described in
section 5.5. Denoting numerical solution as , the benchmark solution as and h as
eitherthetimeorspaceinterval,theLerror,L2errorandthemaxerrorare

Lerror | | h

L error h

maxerror max| | h

Thethreeerrorsillustratedifferentaspectsofaccuracy.TheLerrorgivesequalweight
toallerrors.TheL2errorgivelargeweightstolargeerror,thereforeitprioritizesonthe
overall fitness of the numerical solution to the benchmark. The maxerror penalizes
oscillationsastheseerrorscanbeseveraltimesthemagnitudeofthetruncationerrors.

6.2.2 Massbalanceerror

Animportantaspectofadsorptionmodelingandsimulationistherigorousness
in which mass is conserved. Mass balance is carried out to make sure that the
simulation is as close to reality as possible. Quantities inside the column were
integratedoverthespacedomainwhilequantitiesmeasuredatthecolumnoutletwere
integratedwithtime.Itiscomparedtothatofthebenchmarksolutionandtherelative
error is recorded. Denoting the solution mass balance as MB and benchmark mass
balanceasMB*,theMBerroris

|MB MB|
MBerror
MB

34
 CN4118RFinalReport

where

MB z Z

For quantities inside the column, the integral is evaluated using trapezoidal rule for
FDM grids, while for FVM, the sum of the solution is simply multiplied by the grid
interval,sincethevaluedefinedbyaFVMnoderepresentstheaveragevalueinthecell.

6.2.3 Computationtimetaken

Thisperformanceindexisobtainedonlyinexperiment3.

Thetimetakentoruneachcodeiscompared.MATLABprovidesaprofilertool
whichbreaksdownthesimulationtimeintodetailedanalysisofeachfunctionused.The
timesclockedbytwoselectedfunctions(Table7)andthenumberoftimesodefunwas
calledareofinterest.

Table7Functionsselectedforprofilingandtheirdescriptions
Function Description

ode45 BuiltinMATLABodesolverusingRungeKuttamethod

odefun Functionhandleforode45;containsgoverningequations

The number of times odefun was called depends on the time step size
determinedbyode45andcanaffectthecomputationaltime.Itisputforwardherethat
this number will be larger for the FVM schemes as a penalty for higher accuracy,
especiallyincaseswherethewavefrontissharp.Repetitionsofthesameproblemwere
ranfivetimesandtheaveragetimereported.Thedatawerethencomparedonthebasis
ofequivalentnumberofgridsandequivalenterrorproduced.Therefore,althoughFVM
is penalized by having odefun called more times, it should provide computational
efficiencybyusinglessernumberofgrids.

The data for equivalent error produced was obtained by the following
procedure:
1) Gridsizesobtainedforalistoferrorvaluesbyinterpolationfromerrordata.
2) Computation time obtained for the set of grid sizes from the previous step by
interpolationfromcomputationaltimedata.
3) Computationtimematchedtoerrorproduced.

35
CN4118RFinalReport

6.3 Experiments

6.3.1 Experiment1

The plug flow model was tested here by setting Pe to be infinitely large. The
benchmarksolutionistheinputstepfunctiontranslatedby intothecolumn.Thegrid
sizewasvariedandthetotalsimulatedtimeis0.5 .

Theaimofthisexperimentistoestablishgeneraltrendsandtoevaluatevarious
schemesbasedonamodelwheretheanalyticalsolutionisexact.

6.3.2 Experiment2

Theaxiallydispersedplugflowmodelwasinvestigatedhere.Pefrom10to105
in logarithmic intervals was investigated. The benchmark solution was provided in
section5.5.1.ThedatawascollectedatZ=0.1tosatisfytheassumptionofsemiinfinite
medium.

Theaimofthisexperimentistofindouthowtheintroductionofthedispersion
term,functionofPe,canalterthetrendsobservedinExperiment1.Thehighresolution
schemeswerefurtherdifferentiatedbasedontheirrobustness.

6.3.3 Experiment3

Theaxiallydispersedplugflowmodelwithsourcetermwastestedhere.Pewas
kept large while the lumped parameter mass transfer coefficient (k) was varied. The
benchmarksolutionwasprovidedinsection5.5.2.ThedatawascollectedatZ=0.1to
satisfytheassumptionofsemiinfinitemedium.MB*followssection5.5.3andbyusing
solutionobtainedatimelongaftertheconcentrationfronthadbrokenthrough.Linear
isothermwasused.

One aim of this experiment is to complete the analysis of accuracy and mass
balancebylookingattheentiresetofgoverningequationsinafixedbedadsorber.

Another aim of this experiment is to illustrate visually the capability of FVM

schemesinmodelingsharpwavefrontandtheextenttowhichFDMfailedinthisaspect.
ParametersarechosensuchthatthePDEsrepresentextremecasesandareameasure
ofrobustness.

Finally, the computational time is looked at to seeif FVM schemes provide any
significantadvantageoverFDMintermsofcomputationalefficiency.

36
 CN4118RFinalReport

PARTIII:RESULTS&DISCUSSION3

7. Experiment1

Parametersused:

BedLength FeedVelocity TimeSimulated

L(cm) u0(cm/s) t(s)

100 1 50

7.1 FDM

Figure 10 demonstrates that FDM is not suitable to model sharp wave fronts.
This is because spatial derivatives near the top of the steep slope is overestimated,
resulting in overshoots and generate oscillations. The accuracy of the solution can be
improvedbyusingfinergrids.However,oscillationsstilloccurevenatveryfinegrids.

1 1
N=100 N=2000

0 0
yA/yA0at50seconds

1 1
N=500 N=5000

0 0

1 1
N=1000 Exact
Solution
0 0
0 20 40 60 80 100 0 20 40 60 80 100
Positionincolumn,z(cm)
Figure10FDMonplugflowmodelusingdifferentN

7.2 Firstorderupwindscheme

ThefirstorderupwindschemehasanadvantageoverFDMasitdoesnotcreate
oscillation(Figure11).However,itishighlyinaccuratewithcoarsetomoderatelyfine
grids.Also,asthewavefronttravelsthroughthecolumn,theslopebecomesgentlerand

3Allnumericaldatausedinthispaperisuploadedtohttp://tinyurl.com/FYPdata.

37
CN4118RFinalReport

seemstobehavelikedispersedplugflow,althoughitoughttoremainsharp.Thisisthe
dissipativeeffectoftheupwindscheme.

Concentrationprofilenearwavecenteratdifferenttimes
1

At10s

At30s
yA/yA0

At50s

ExactSolution

0
25 20 15 10 5 0 5 10 15 20 25

Distancefromwavecenter(cm)

Figure11Dissipativeeffectoffirstorderupwindscheme
FirstorderupwindschemeonplugflowmodelusingdifferentN
1

N=100
yA/yA0

N=500
Att =50sec N=2000
N=5000
ExactSolution

0
30 35 40 45 50 55 60 65 70

Positionincolumn,z(cm)

Figure12Effectsofgridsizeonfirstorderupwindscheme

Theaccuracyoftheupwindschemecanbeimprovedbyusingfinergridaswell
(Figure12).AtN=5000,thesolutionoftheupwindschemecanbecomparabletothat
ofFDMwithouttheoscillations.Thedilutionofthewavefrontcanbejudgebylooking
athowfarthesolutionspreadsfromthewavecenter(e.g.~40cmforN=100).

38
 CN4118RFinalReport

7.3 SecondorderandhighresolutionFVMschemes

The secondorder and high resolution schemes have shown significant

improvementsoverthefirstorderupwindscheme.Usingthesamenumberofgrids,the
sharpness of the wave front is captured; solutions are spread over a shorter distance
fromthewavecenter(~20cm)ascomparedtothatoftheupwindscheme.However,
allthreesecondorderschemesareoscillatory(Figure13(a)).

The downwind differencing LaxWendroff scheme has oscillation lagging the

wavecenter,similartoFDM.TheupwinddifferencingBeamWarmingschemeandthe
central differencing Fromm scheme both have oscillations at the front, although for
Frommscheme,theoscillationsarelesssevere.

Asforthehighresolutionschemes,thedifferencebetweenthemcanbeobserved
byfocusingontheareanearthewavecenter(Figure13(b)).TheSuperbeeschemegave
thebestaccuracy,followedbyMCandvanLeer.

7.4 Comparisonoferrors

Figure14,Figure15andFigure16showtheloglogplotsoftheerrorsagainst
gridsizeZ.Asamplecalculationfortheerrorswasgivenbelow.

Table8SampledataforfirstorderupwindschemewithN=10
Z yA/yA0 (yA/yA0)* (yA/yA0)* yA/yA0
0 1 1 0
5 0.9933 1 0.0067
15 0.9596 1 0.0404
25 0.8753 1 0.1247
35 0.7350 1 0.2650
45 0.5595 1 0.4405
55 0.3840 0 0.3840
65 0.2378 0 0.2378
75 0.1334 0 0.1334
85 0.0681 0 0.0681
95 0.0318 0 0.0318
100 0.0318 0 0.0318

39
CN4118RFinalReport

Secondorderschemesonplugflowmodel(N=100)

BeamWarming
yA/yA0

LaxWendroff
Fromm
At t=50sec
ExactSolution

0 10 20 30 40 50 60 70 80
0.2
Positionincolumn,z(cm)

(a)
Highresolutionschemesonplugflowmodel(N=100)
1

vanLeer
yA/yA0

Superbee
At t=50sec MC
ExactSolution

0
45 50 55
Positionincolumn,z(cm)
(b)
Figure13(a)SecondorderFVMschemes(b)HighresolutionFVMschemes

Lerror 0.1 0 0.0067 0.0404 0.1247 0.2650 0.4405 0.3840 0.2378

0.1334 0.0681 0.0318 0.0318 0.17642

L error 0.1 0.0067 0.0404 0.1247 0.2650 0.4405 0.3840

.
0.2378 0.1334 0.0681 0.0318 0.0318 0.07141

maxerror 0.1 0.4405 0.04405

40
 CN4118RFinalReport

The loglog plots summarize the findings made in the preceding sections. The
positiveslopesindicateapositiverelationshipbetweenthetwovariables.Thegenerally
lineartrendssuggestthattheyfollowapowerlawrelationship.Moreaccurateschemes
arefoundlowerintheplot.TheSuperbeeschemeisthemostaccurateacrossallthree
plots,followedbytheothertwohighresolutionschemes,thesecondorderschemesand
finally the firstorder upwind scheme. The LaxWendroff scheme and FDM share the
exactsameplot,confirmingthatbothalgorithmsareactuallyequivalent.

Themeritsofhighresolutionschemesarecleartosee;onlytheseschemeswould
be selected to test against the FDM in subsequent experiments. Although the overall
errors (Lerror) of upwind converge towards the FDM at fine grids, it is still far less
accurate than FDM if the overall fitness (L2error) is considered instead. It is much
moredifficulttousetheupwindschemetocapturethesharpwavefrontthattheother
schemescaneasilyproduce.

When using the maxerror (Figure 16), the ranking is reversed. Upwind is
deemedtobemoreaccuratethanmostschemesatthesmallestgridsize;onlySuperbee
is more accurate. Schemes which produce oscillatory solutions have the largest max
error.

7.5 Massconservation

Figure 17 shows the loglog plot of mass balance error against grid size Z. A
samplecalculationfortheMBerrorwasgivenbelow.

ForFVMdefinednodes,usingdatafromTable8,

yA,
MB Z 4.9778 10 49.778
yA0

MB 50forasimulatedtimeof50seconds.

|50 49.778|
MBerror 4.44 10
50

ForFDMdefinednodes,usingdatafromTable9,

yA, yA,
MB Z 5.01705 10 50.1705
2yA0

MB 50forasimulatedtimeof50seconds.

41
CN4118RFinalReport

|50 50.1705|
MBerror 3.41 10
50
Table9SampledataforFDMschemewithN=10
Z yA/yA0
0 1
10 1.0429
20 1.1310
30 1.0056
40 0.6932
50 0.3797
60 0.1710
70 0.0650
80 0.0219
90 0.0055
100 0.0025

log(Lerror)againstlog(Z)
0
5 4 3 2 1 0

(1)
(2)& (3)
log(Lerror)

2
(1) Upwind
(4) (2) FDM
(5) (3) LaxWendroff
(6) (7) (4) BeamWarming
(5) Fromm
(6) vanLeer
(7) MC
(8) Superbee
(8)
4
log(Z)

Figure14ComparisonofLerrorsbetweendifferentschemes

42
 CN4118RFinalReport

log(L2error)againstlog(Z)
1
4.5 3.5 2.5 1.5 0.5

2
(1)
log(L2error)

(1) Upwind
(2) FDM
(3) LaxWendroff
(4) BeamWarming
3
(5) Fromm
(2),(3)&(4)
(6) vanLeer
(5) (7) MC
(6)&(7) (8) (8) Superbee
4
log(Z)
Figure15ComparisonofL2errorsbetweendifferentschemes

log(maxerror)againstlog(Z)
1
4 3 2 1

(2)&(3)

2
log(maxerror)

(1) Upwind
(2) FDM
(3) LaxWendroff
(4) BeamWarming
3
(5) Fromm
(6) vanLeer
(8)
(7) MC
(8) Superbee
4
log(Z)
Figure16Comparisonofmaxerrorsbetweendifferentschemes

The data illustrate another aspect of FVM schemes mass conservation. FVM
schemeswhicharetotalvariationdiminishing(TVD)(i.e.firstorderupwind,vanLeer,
MC and Superbee) conserve mass rigorously [35]. The MBerrors of these schemes
quickly settle onto a plateau. Other methods, including the FDM, follow a downward

43
CN4118RFinalReport

trendasthegridrefines.ItcanbeseenfromFigure17thatthemassbalanceofFDMis
much less accurate than that of FVM, even though its degree of accuracy (~108) is
usuallyacceptable.Theimprovementatveryfinegridisalsominimal.Theshapesofthe
MBerrorloglogplotswouldbeunderstoodfurtherinExperiment2.

log(MBerror)againstlog(Z)
0
4 0
(3)

(5)
5
log(MBerror)

(2)
(1) Upwind
(2) FDM
(1)
(3) LaxWendroff
(4) (4) BeamWarming
10
(6) (5) Fromm
(8) (6) vanLeer
(7) MC
(7) (8) Superbee
15
log(Z)
Figure17ComparisonofLerrorsbetweendifferentschemes

8. Experiment2

Parametersused:

BedLength FeedVelocity TimeSimulated PecletNumber

L(cm) u0(cm/s) t(s) Pe

100 1 50 10to105

8.1 Generaltrend

TheeffectofPeontheaccuraciesisillustratedinFigure18(a),whichistruein
generalforalloftheschemesinvestigated.Thenumericalmethodsaremoreaccurate
whenPeissmaller.Thisistobeexpectedasdiffusionspreadsoutthewavefrontand
thenumericalmethodscanmakemoreaccurateapproximations.

ThedataforwhenPeisequalto10buckstheobservedtrend.Thisisprobably
due to assumptions made by the analytical solution not being fulfilled. The analytical
solution requires that the outflow boundary condition have no effect on the solution.
ThiscouldbeuntrueforthecasewherePe=10becauseatthestipulatedtime(t=50s),
44
 CN4118RFinalReport

the diffused wave front has reached the end of the column. Taking this into account,
dataobtainedforPe=10wasignored.

Interestingly, there are changes to the effect of Z on the errors. There is a

threshold grid size beyond which the pattern observed in Experiment 1 is disrupted
(Figure 18(b)). It most readily occurs when Pe is small (10 to 100) and N is at least
1000.Thethresholdgridsizeseemstobedifferentfordifferentschemes.Unfortunately,
as most data fall into the range indicated by Figure 18(b)(ii), this change is not well
understood. It could be due to algorithm errors in calculating the benchmark solution
becauseoftheimmenselysmallgridsize,oritcouldbeduetonumericalinstability.

8.2 Superbee

The Superbee scheme does not follow the general trend described above. The
loglogplotfollowsanirregulardownwardstrendandhasacharacteristicbump.This
bump occurs at different grid size for different values of Pe (Figure 19). It occurs at
largergridsizeasPedecreasesandasaresult,theaccuracyoftheschemesuffered.

The reason for this is because the Superbee scheme tends to sharpen smooth
slopes[35].ThistendencygaveititsadvantageinExperiment1wherethewavefrontis
a sharp discontinuity. However, as the diffusion effect becomes more dominant, the
smooth wave front is incorrectly sharpened. The extent of sharpening changes as the
gridisbeingrefined.Therefore,theSuperbeeschemeisunsuitableforsituationswhere
thewavefrontisspreadoutacommonoccurrenceinanadsorber.

8.3 EffectofPeonaccuracyofFDM

As mentioned earlier, the accuracies of all schemes improved as Pe decreases.

The extent of improvement by refining the grid is greater for FDM than for the FVM
schemes. The FVM schemes are still much more accurate than FDM for Pe at 105, 104
andmarginallyat103(Figure20(a)to(c)),FDMovertooktheFVMschemeswhenthe
Peis100(Figure20(d)).

This result is not restricted to the relative dominance of diffusion over

convection. With the introduction of the source term, the concentration wave front is

45
CN4118RFinalReport

EffectofPeonaccuracy
Inf 0
log(Error)

Pe=Inf
Pe=1e5
Pe=1e4
Pe=1e3
Pe=100
Pe=10

log(Z) Inf

(a)
EffectofZonthethreeerrors
Inf 0

(i) Rangewhere (ii) Rangewhere

overallerror refiningthegridstill
becomesstagnant; haspositiveeffect
slopesbecomeless onaccuracy
positive
log(Error)

Lerror
L^2error
maxerror

log(Z) Inf

(b)
Figure18Generaltrendobservedinloglogplotsforexperiment2

46
 CN4118RFinalReport

EffectofPeonaccuracyofSuperbeescheme
Inf 0

(b)
(a)
log(Error)

Pe=Inf
Pe=1e5
(d) Pe=1e4
(c) Pe=1e3
Pe=100

log(Z) Inf

Figure19TrendobservedforSuperbeescheme
Locationofbump:(a)Pe=1e5atZ=0.00167;(b)Pe=1e4at0.005;(c)Pe=1e3at0.0167;(d)
Pe=100at0.05;

expectedtospreadoutifthelumpedparametermasstransfercoefficientissmall,even
at extremely large Peclet numbers. As Figure 20(d) has shown, FDM can be more
accuratethantheFVMschemesundersuchcondition.

8.4 EffectofPeonmassbalances

The mass balances were better for FDM at smaller Pe as well (Figure 21). The
shape of the loglog plots for FVM schemes are unremarkably the same as that in
experiment 1; the errors rapidly settled onto a plateau and did not improve over the
range of grid size investigated. For FDM, the errors decreased steadily before settling
onto a plateau as well. The order of magnitude at which the plateaus are located
dependedonPe.Itrangesfrom102to104overtherangeofPeinvestigatedfortheFVM
schemes,comparedto102to1011forFDM.

In terms of practicality, ultrahigh order of accuracy (i.e. 1011) may not be

required.ThemoreimportantaspectisthattheFVMschemesarerobustenoughthat
theirmassbalanceerrorsareguaranteedtobewithinreasonablerangeforanygridsize
used.Thiscouldbeanimportantconsiderationifrobustnessisrequired.

47
CN4118RFinalReport

0 0
3 2.5 2 1.5 1 3 2.5 2 1.5 1
Pe=1e5 Pe=1e4
log(Lerror)

FDM
vanLeer
MC
Superbee
4 4
log(Z)

(a) (b)
0 0
3 2.5 2 1.5 1 3 2.5 2 1.5 1
Pe=1e3 Pe=1e2

5 6

(c) (d)
Figure20ImprovementinFDMoverFVMasPe decreases
Lerrorsareplottedhere;L2errorsandmaxerrorsshowsimilartrends.

2 2
3 2.5 2 1.5 1 3 2.5 2 1.5 1
FDM
vanLeer
log(MBerror)

Pe=1e4
Superbee
MC

Pe=1e5
6 4
log(Z)

(a) (b)
0 0
3 2.5 2 1.5 1 3 2.5 2 1.5 1

Pe=1e3 Pe=1e2
12 12

(c) (d)
Figure21Loglogplotsofmassbalanceerrorsagainstgridsizeforexperiment2

48
 CN4118RFinalReport

9. Experiment3

Parametersused:

Feed Peclet LumpedParameterMass Henrys

BedLength
Velocity Number TransferCoefficent Constant
L(cm)
u0(cm/s) Pe k(s1) K

100 1 1010 10to0.01 10

9.1 GeneralTrend

Thetrendobservedinthisexperiment(Figure22)issimilartothatobservedin
experiment 2. Errors become smaller as k decreases, confirming that the numerical
methodsaremoreaccurateastheconcentrationfrontspreadsout.Thesamenegative
effect of excessive spreading on the accuracy (akin to when Pe = 10 in experiment 2)
wasobservedinthisexperimentfork=0.01;dataobtainedfork=0.01wouldnotbe
presentedsubsequently.

AcomparisonbetweenthedifferentschemesisshowninFigure23.Asexpected,
the advantages of FVM schemes over FDM diminish as the wave front becomes more
spreadout,oraskbecomessmaller.However,evenatthesmallestkvalueinvestigated,
the FVM schemes (Superbee scheme excluded) are still more superior. The patterns
seeninFigure23aregenerallytrueforL2errorsandmaxerrorsaswell.

9.2 Effectsofkonmassbalances

Figure 24 shows the loglog plots of MBerrors against grid size. The data for
when k = 0.01 is included because the analytical solution developed for mass balance
operateondifferentassumptions(seesection5.5).Themassbalanceerrorforwhenk=
10ischaoticandcanbeattributedtotheverysharpwavefrontexitingthecolumn.As
theexitingwavefrontbecomemuchgentler,familiarpatternsemerged.

TheMBerrorsforFDMfollowedadownwardtrendasthegridwasrefined.The
firstlocalminimumseenon(c)(nearlog(Z)=1.2)isduetooscillationsatisfyingthe
analyticalsolutioncoincidentally.Theotherlocalminimumandalsotheonein(d)are
duetoFDMswitchingfromoverestimatingthemassbalancetounderestimatingit.

49
CN4118RFinalReport

Effectofkonaccuracy
Inf

0
log(Error)

k=10
k=1
k=0.1
k=0.01

log(Z) Inf

Figure22Generaltrendobservedinloglogplotsforexperiment3

k=10 k=1
log(Lerror)

FDM
vanLeer
MC
Superbee

3 2.5 2 1.5 3 2.5 2 1.5 0.5 1

0.5 1
log(Z)
(a) (b)

k=0.1

3 2.5 2 1.5 0.5 1

(c)
Figure23 LoglogplotsofL errorsagainstgridsizeforexperiment3

50
 CN4118RFinalReport

The MBerrors for the FVM schemes were observed to settle onto a plateau
quicklybeforebecomingmoreoscillatoryafteracertainthreshold gridsize.Thismay
be due to numerical instability arising from the small grid sizes. Under numerically
stablecondition,theFVMschemesareconsistentlyunderestimatingthemassbalance,
anotherpropertyofTVDschemes.

0 2
3.5 k=10 2.5 1.5 3.5 k=1 2.5 1.5
log(MBerror)

FDM
vanLeer
MC
Superbee 8 8
log(Z)

(a) (b)
2 3.5
3.5 k=0.1 2.5 1.5 3.5 k=0.01 2.5 1.5

6 6

(c) (d)
Figure24Loglogplotsofmassbalanceerrorsagainstgridsizeforexperiment3

9.3 Robustness

The robustness of FDM and FVM are tested by using parameters representing
extreme cases. Because of the closeness between the FVM schemes, data obtained via
the van Leer scheme is chosen to represent the FVM. The simulated solutions for 3
extreme cases are presented in Figure 25, along with the relevant parameters. The
descriptionforeachcaseisgivenbelow:

(a) Highlumpedparametermasstransfercoefficientforlinearisotherm

(b) Highlumpedparametermasstransfercoefficientforhighlynonlinearisotherm,
withsaturationconstantadjustedtoanextremelylowlevel

(c) Highlynonlinearisothermwithahighersaturationconstantthan(b)

51
CN4118RFinalReport

The solutions obtained via FDM have proved to be untenable in each of these three
cases whereas FVM could provide smooth solutions. Although it is complicated to
validatethesesolutionsanalytically,thoseobtainedviaFDMwerewhollyimpossiblein
thephysicalsense(i.e.gasphaseconcentrationhigherthanfeedconcentration).

In all three cases investigated above, the oscillations were sustained until
breakthrough.Thisisduetotheconcentrationfrontbeingsteepfortheentireduration.
In the event where the front will disperse to a gentle slope, the FDM can produce
oscillationaswellduringtheearliestperiodofsimulationwhenthefrontisstillsteep.
Therefore, even if the final breakthrough solution looks smooth, it may be inaccurate
duetotheoscillationswhichhadoccurredearlier.ThisisillustratedinFigure26,where
snapshots of the concentration front were taken at different times using (a) FDM and
(b)vanLeerschemes.

9.4 Computationtimetaken

9.4.1 Equivalentnumberofgrids

Threeplotsforeachvalueofkarepresentedhere:

(a) timesclockedbyode45

(b) numberoftimesodefunwascalled

(c) thefractionoftimeclockedbyodefunagainstode45

ThelastplotwouldgiveanindicationofhowthebehaviorofFVMschemeschangedask
increases(orasthewavefrontbecomesteeper).Figure27,Figure28andFigure29are
theplotsfork=0.1,1and10respectively.

Figure 27 shows that when the wave front is gentle, the two methods behave
similarly with FVM clocking in longer times than FDM in general. Instability for the
Superbee scheme occurred within the range of grid sizes investigated (Figure 27(b)),
whichisindicatedbyaspikeinthenumberoftimesodefunwascalledataround900
grids.

In Figure 28, the differences begin to become more obvious, with longer times
clockedbyFVMmainlyduetoinstabilityatlargegridsizes.However,thebehaviorsof
thedifferentschemesremainedsimilar(Figure28(c)).

52
 CN4118RFinalReport

0.015

0.01
vanLeer
yA

FDM
0.005

0
0 20 40 60 80 100 120

Positionincolumn,z(cm)

(a) k=100;b=2.5e4;z=1
0.015

0.01

vanLeer
FDM
0.005

0
0 20 40 60 80 100 120

(b)k=100;b=1e9;qs=1e9;z=1
0.025

0.02

0.015
vanLeer

0.01 FDM

0.005

0
0 20 40 60 80 100 120

(c)k=1;b=1e9;qs=1e6;z=1
Figure25Simulatedsolutionsforextremecases

53
CN4118RFinalReport

0.01

0.005
yA

0
0 20 40 60 80 100 120

Positionincolumn,z(cm)

(a)FDM
0.01

0.005
yA

0
0 20 40 60 80 100 120

Positionincolumn,z(cm)

(b)vanLeer
Figure26ComparisonbetweensnapshotsofFDMandvanLeerschemes

54
 CN4118RFinalReport

40
FDM
Timeclockedbyode45(s)
vanLeer
30
MC
Superbee
20

10

0
0 200 400 600 800 1000

Numberofgrids

(a) Timesclockedbyode45againstnumberofgrids
32000
Numberoftimesodefunwas

31500
called

31000

30500

30000
0 200 400 600 800 1000

Numberofgrids

(b)Numberoftimesodefunwascalledagainstnumberofgrids
0.88
Timesclockedbyodefunover
ode45

0.8

0.72
0 200 400 600 800 1000

Numberofgrids

(c)Timesclockedbyodefunoverode45
Figure27Computationaltimeplotsfork=0.1

55
CN4118RFinalReport

80
Timeclockedbyode45(s)
60
FDM
vanLeer
40
MC
Superbee
20

0
0 200 400 600 800 1000

Numberofgrids

(a) Timesclockedbyode45againstnumberofgrids
50000
Numberoftimesodefunwas
called

40000

30000
0 200 400 600 800 1000

Numberofgrids

(b)Numberoftimesodefunwascalledagainstnumberofgrids
0.88
Timesclockedbyodefunover
ode45

0.8

0.72
0 200 400 600 800 1000

Numberofgrids

(c)Timesclockedbyodefunoverode45
Figure28Computationaltimeplotsfork=1

56
 CN4118RFinalReport

1000
Timeclockedbyode45(s) FDM
vanLeer
MC
Superbee
500

0
0 200 400 600 800 1000

Numberofgrids

(a) Timesclockedbyode45againstnumberofgrids
210000
Numberoftimesodefunwas

190000
called

170000

150000
0 200 400 600 800 1000

Numberofgrids

(b)Numberoftimesodefunwascalledagainstnumberofgrids
0.88
Timesclockedbyodefunover
ode45

0.8

0.72
0 200 400 600 800 1000

Numberofgrids

(c)Timesclockedbyodefunoverode45
Figure29Computationaltimeplotsfork=10

57
CN4118RFinalReport

Thelargestdifferenceintimesclockedisobservedwhenthewavefrontissteep
(Figure29).ThetimesclockedforFVMschemesincreasedexponentiallywhilethatfor
FDM increased almost linearly. There is a marked difference in the number of times
odefunwascalledbytheFVMschemesrightfromthelowestnumberofgrids.Thisisto
be expected as ode45 have to use smaller time steps for the steep front to propagate
properly.ThisisnotnecessaryforFDMastheoscillationssmoothenthefrontandallow
larger time steps to be used. Also, the behaviors of FVM schemes deviated from FDM
(Figure29(c)).Withlargenumberofgrids,ahigherproportionoftimewasspentinthe
otherfunctionsapartfromodefun.Thisisanindicationthatthesolverisspendingmore
timeoptimizingthetimestepinsteadofdoingtheactualcomputinginodefun.

9.4.2 Equivalenterror

Thefiguresareobtainedbymatchingtheloglogplotsoferrorsagainstgridsize
withthetimesclockedateachparticulargridsize.Figure30summarizesthefindings.

ItisobviousthatFVMschemesarelessefficientthanFDMfork=0.01.However,
the deficiency is not significant in terms of absolute amount of time taken. It is put
forwardherethatthediscrepanciesareentirelyduetotheadditionalcodesrequiredby
the delimiter functions, which would not contribute to improved accuracy since the
wavefrontisgentle.Onthewhole,FVMschemesareabout25%lessefficientthanFDM.

Askincreases,theadvantagesofFVMbecomeprominent.Atk=1,allschemes
alreadyhavecomparablecomputationefficiencies(FVMhas~10%improvementover
FDMovermidrangeaccuracy).Athighaccuracy,FVMschemes(withtheexceptionof
Superbee)becomesignificantlymoreefficient(~3050%for3datapoints)thanFDM.
TheSuperbeeschemedisplayedspikesincomputationtimeathighaccuraciesinFigure
30(a)and(b).

Whenthewavefrontisverysteep(i.e.k=10),FVMschemesaredefinitelymore
efficient than FDM, with Superbee being the most efficient (up to 55% improvement
over FDM). The spikes seen in (a) and (b) did not appear in (c) for the range
investigated. At the highest accuracy investigated, van Leer scheme become less
efficientthanFDM,althoughtheothertwoschemesremainedmuchmoreefficientthan
FDM.

58
 CN4118RFinalReport

20
FDM
k=0.1 vanLeer
Computationtime(s)

MC
Superbee
10

0
0.5 0 0.5 1 1.5 2 2.5

log(Lerror)

(a)
50
Computationtime(s)

k=1

25

0
0.5 0 0.5 1 1.5 2 2.5

log(Lerror)

(b)
180
Computationtime(s)

k=10

90

0
0.2 0.7 1.2 1.7 2.2

log(Lerror)

(c)
Figure30Computationtimebyequivalenterrorproduced

59
CN4118RFinalReport

10. Improvements

The numerical experiments conducted in this project were limited in its scope.
Thereareafewdirectionsinwhichthisinvestigationcanproceed.

First, a different ODE solver can be used in MATLAB. MATLAB provide several
ODE solvers for different problem types (stiff or nonstiff). One of them, ode15s, can
solvestiffproblemsfasterthanode45atreasonableorderofaccuracy.Ideally,anODE
solvershouldbewrittenspeciallyfortheFVMsuchthatthetimestepsizecanbeused
in defining the time derivative (see section 6.1). This would allow a more accurate
representation of the algorithm and may reduce the minimum grid size where
numerical instability occurs. It would also allow fractional step methods to be
implementedforthesourceterms.Admittedly,writingarobustODEsolverinMATLAB
can be a monumental task (ode45 consist of 400 over lines of codes while ode15s has
900).

Second, this project limits itself to fundamental numerical methods, explicit

schemes and fixed grids, in solving the problems. Advance methods such as using
adaptive grids may be useful in adsorption problems where sharp fronts occur,
providing both accuracy and efficiency [27]. The implicit CrankNicolson method is
commonlyusedtosolvediffusionproblemswhileafractionalstepimplicitmethodcan
be applied to convectiondiffusionreaction problems [49]. These methods are left
unexploredinthisproject.

Third,thetrendsobservedforcomputationaltimetakenonlyconsiderkasthe
independentvariable.Thereareamultitudeofparametersinthegoverningequations
andcanallhaveaneffectontheefficiencyofthecodes.Atrulythoroughanddefinitive
assertionthatFVMismoreefficientthanFDMintheadsorptionproblemmusttakeall
parametersintoaccount.

Lastly,thenonisothermalcodesimplementedintheGUIwerenotvalidatedwith
experimentaldata.Thesolutionsobtainedviathosecodes,thoughnotpresentedhere,
are identical to the isothermal ones by configuring parameters to simulate isothermal
conditions.Howeverthisonlyvalidatesaportionofthecodes.Experimentaldatawith
smallfluctuationsintemperatureisrequiredtovalidatetheentirealgorithm.

60
 CN4118RFinalReport

PARTIV:DEVELOPMENTOFGUI

11. GUI

ThegraphicuserinterfacewascreatedusingMATLABGUIdesignenvironment
(GUIDE).GUIDEallowstheprogrammertocreateapplicationswhichcanbedistributed
royaltyfree with the MATLAB compiler runtime (MCR). The package can be installed
andtheapplicationwillworkoncomputerswithoutacopyofMATLAB.Thusfar,only
computersrunningonMicrosoftWindowsOSareverifiedtobeabletorunthepackage
successfully.

ApartfromusingGUIDE,analternativemethodtocreatingaGUIwasexplored.
MATLAB Builder NE allows the programmer to deploy .NET components written as
functions in MATLAB. These components can then be used in more advance
programming languages such as Microsoft Visual Basic or C#. However, interfacing
betweenMATLABandVisualBasichadbeenunsuccessfulandtheremedysoughtfrom
theonlineMATLABcommunitywasverytedious.Nonetheless,thismethodcouldbere
explored if and when MATLAB develops a more robust interface process in future
releases.

The motivation behind the creation of this GUI is to assist in the teaching of
fundamentalsinadsorption.Assuch,accuratedepictionofparameterchangesismore
important than producing highly accurate simulation results. The application should
also be able to project visual information directly using data obtained from the
simulation. As of the time of writing, the application allows the user to formulate an
adsorberusingequations(25)to(28)and(11),withthefollowingcapabilities:

saveandloadmultipletxtfilescontainingparametricinformation

selectandexportgasphaseconcentration,velocity,temperatureandsolidphase
loadingdatato.txtor.xlsformat

plot all aforementioned variables dynamically during simulation as column

interiordata

plot all aforementioned variables (excluding solid phase loading) at the end of
simulationasbreakthroughdata

speedcontrolduringdynamicplotting

61
CN4118RFinalReport

12. AlgorithmSpecifictoNonisothermalEquations

The simulator would be using the FVM van Leer scheme. Boundary conditions
for nonisothermal equations (29) and (30) are used to define values at the left and
rightinterfacesrespectively.

Forgasphaseconcentration:

ZyA 2yA
yA , yA yA0
Z 2

Forgasandsolidphasetemperature:

T T 4 5 Z
T ,
2 4 2 T T

Forinterstitialgasvelocity: 1

Forcorrespondingoutflowboundaries:yA yA , T T , Tw Tw

5
istheequivalentPecletnumberforT.
4 T T

Ta 1is one of the simplifying assumptions made. It is valid when the operating
temperature is near ambient temperature, or the ambient temperature is made to be
the same as the operating temperature. OnlyCT or g is allowed to change with
temperaturewhileallotherparameters,whichshouldbefunctionsoftemperature,are
assumedtobetemperatureindependent.

62
 CN4118RFinalReport

PARTV:CONCLUSION

Theresultshavesubstantiatedconventionalwisdomabouttheadvantagesthat
FVMprovideoverFDM.FVMissuperiorinbothaccuracyandefficiencywhenthereisa
sharp discontinuity. FVM also ensure mass balance robustly, regardless of the
parameters. On the other hand, FDM remains more efficientat handling smooth wave
fronts whereas in such cases, the computational overheads incurred by the more
complicatedFVMcodesarenotutilizedeffectively.FDMssimplicityinformulationcan
alsoallowanoviceprogrammertopickupthemethodquickly.

Several FVM schemes were investigated. The Superbee scheme is the most
suitable if the sharp discontinuity persists throughout the entire simulation. Van Leer
andMCschemesaresimilartooneanotherandaremoreusefulforgeneralcases.Many
advanceschemesmaybemoreusefultotheadsorptionproblem[26],althoughitmay
beimpracticaltoimplementthemasanindividualeffort.

The limitation of the MATLAB ODE solver is acknowledged but it must be

emphasized that the numerical experiments were designed to glimpse at the general
trendsandnottoquantifytheresultsdefinitively.Otherimprovementstothenumerical
experimentswerealsodiscussed.

Acknowledgment

TheauthorwouldliketothankProfessorShamsuzzamanFarooq,mentortothisproject,
for his patience and guidance, without which so much could not have been learnt in
suchlittletime.

63
CN4118RFinalReport

NotationsandAbbreviations

CFD ComputationalFluidDynamics
CMS CarbonMolecularSieve
FDM FiniteDifferenceMethod
FEM FiniteElementMethod
FVM FiniteVolumeMethod
MEA Monoethanolamine
PSA PressureSwingAdsorption
SMB SimulatedMovingBed
TSA ThermalSwingAdsorption
VSA VacuumSwingAdsorption
b, b0 LangmuirconstantandArrheniusparameter
Co Courantnumber
C, CT Gasphaseconcentrationandtotalconcentration
Cpg , Cps , Cpa Heatcapacitiesofgas,solidandadsorbedphases
Dcol Columndiameter
DL Longitudinaldispersioncoefficient
dp Particlesize
H Gibbsfreeenergyforadsorption
hi , ho Innerandouterconvectiveheattransfercoefficient
K HenrysconstantforLangmuirisotherm
Kw Wallconductiveheattransfercoefficient
Kz Longitudinalconductiveheattransfercoefficient
k Lumpedparametermasstransfercoefficient
ka , kd Rateconstantofadsorptionanddesorption
N Numberofcellsoramountofspecies
Pe Pecletnumber
P Pressure
q, qs , q Amountadsorbed,maximumcapacityandequivalentamountadsorbed
R Molargasconstant
r Argumentfordelimiterfunctions(seesection5.4.4)
ri , ro Innerandouterradiiofcolumn
S Columncrosssectionalarea
T,T ,T Gas/solidphase,wallandambienttemperatures
u Interstitialvelocity

Greekletters Subscripts
Bedvoidage 0 Initialcondition
,
g s Densityofgasandsolidphases A ComponentA
Fractionofoccupiedsites I Inertgascarrier
, Arbitraryquantity i Gridcelli
, Arbitrarycoefficient L,l Leftboundaryandleftinterface
, Seesection5.2.2 R,r Rightboundaryandrightinterface
Seesection5.2.3 z Longitudinaldirection
1 to8 Seesection5.2.3
, Seesection5.5.2 Superscripts
Seesection12 n Timestepn

Dimensionlessquantities(Variablesaccentedwithbars arenotlistedhere)
Time
xi Amountadsorbed
yi Gasphaseconcentration

64
 CN4118RFinalReport

References

1. Baron, G.V., Pressure Swing Adsorption: Latest Developments, Economic Aspects

and Modelling, in Gas Separation Technology, E.F. Vansant and R. Dewolfs,
Editors.1990,ElsevierSciencePublishersB.V.:Amsterdam.p.137148.
2. Humphrey, J.L. and G.E. Keller, SeparationProcessTechnology. 1997, New York:
McGrawHill.
3. Colls,J.,AirPollution.2002,London:SponPress.
4. Steinigeweg,S.andJ.Gmehling,Distillation,inGreenSeparationProcesses,C.A.M.
AfonsoandJ.G.Crespo,Editors.2005,WileyVCH:Weinheim.p.127154.
5. Castle, W.F., Air separation and liquefaction: recent developments and prospects
for the beginning of the new millennium. International Journal of Refrigeration,
2002.25:p.158172.
6. Geankoplis, C.J., Transport processes and separation process principles. 4th ed.
2003,NewJersey:PrenticeHall.
7. Kim, S., H.T. Kim, and B.C. Choi, Optimization of CO2 Absorption Process with
MEASolution.StudiesinSurfaceScienceandCatalysis,2004.153:p.429434.
8. Cousins,A.,L.T.Wardhaugh,andP.H.M.Feron,Analysisofcombinedprocessflow
sheetmodificationsforenergyefficientCO2capturefromfluegasesusingchemical
absorption.EnergyProcedia,2011.4:p.13311338.
9. Yu,Y.S.,etal.,AninnovativeprocessforsimultaneousremovalofCO2andSO2from
flue gas of a power plant by energy integration. Energy Conversion and
Management,2009.50:p.28852892.
10. Ohashi, Y., T. Ogawa, and N. Egami, Development of Carbon Dioxide Removal
SystemfromtheFlueGasofCoalFiredPowerPlant. Energy Procedia, 2011. 4: p.
2934.
11. Kemper, J., G. Ewert, and M. Grnewald, Absorption and regeneration
performance of novel reactive amine solvents for postcombustion CO2 capture.
EnergyProcedia,2011.4:p.232239.
12. Nunes,S.P.andK.V.Peinemann,MembraneTechnologyintheChemicalIndustry.
2006,Weinheim:WileyVCH.
13. Yampolski,Y.,I.Pinnau,andB.Freeman,eds.MaterialsScienceofMembranesfor
GasandVaporSeparation.2006,Wiley:NewJersey.
14. Stookey,D.J.,GasseparationMembraneApplications,inMembraneTechnology:in
theChemicalIndustry,S.P.NunesandK.V.Peinemann,Editors.2006,WileyVCH:
Weinheim.
15. Yang, R.T., Gas Separation by Adsorption Processes. 1987, Boston: Butterworth
Publishers.
16. Tien, C., Adsorption Calculations and Modeling. 1994, Boston: Butterworth
Heinemann.
17. Wankat,P.C.,SeparationProcessEngineering.2nded.2007,UpperSaddleRiver,
NJ:PrenticeHall.
18. Gray, P.G. and D.D. Do, MulticomponentSorptionPhenomenaofSO2andCO2on
Activated Carbon, in Gas Separation Technology, E.F. Vansant and R. Dewolfs,
Editors.1990,ElsevierSciencePublishersB.V.:Amsterdam.p.223230.
19. Chaffee, A.L., et al., CO2 capture by adsorption Materials and process
development.InternationalJournalofGreenhouseGasControl,2007.1:p.1118.
20. Fulker,R.D.,Adsorption,inProcessesforAirPollutionControl,G.Nonhebel,Editor.
1972,CRCPress:BocaRaton,FL.

65
CN4118RFinalReport

21. Basmadjian, D., The Little Adsorption Book: A Practical Guide for Engineers and
Scientists.1997,BocaRaton,FL:CRCPress.
22. Pilarczyk, E. and H.J. Schrter, NewPSAProcesseswithCarbonMolecularSieves
for Recovery of Carbon Dioxide andMethane, in Gas SeparationTechnology, E.F.
Vansant and R. Dewolfs, Editors. 1990, Elsevier Science Publishers B.V.:
Amsterdam.p.271280.
23. Sjostrom, S., et al., Pilot test results of postcombustion CO2 capture using solid
sorbents.EnergyProcedia,2011.4:p.15841592.
24. Yang, R.T., Adsorbents: Fundamentals and Applications. 2003, New Jersey: John
Wiley&Sons.
25. Liu, Y., et al., Recent developments in novel sorbents for flue gas clean up. Fuel
ProcessingTechnology,2010.91:p.11751197.
26. Versteeg, H.K. and W. Malalasekera, An introduction to computational fluid
dynamics:Thefinitevolumemethod.1995,Essex:LongmanScientific&Technical.
27. Chung, T.J., ComputationalFluidDynamics. 2nd ed. 2010, New York: Cambridge
UniversityPress.
28. Dudukovi, M.P. and H.S. Lamba, Solutionofmovingboundaryproblemsforgas
solid noncatalitic reactions by orthogonal collocation. Chemical Engineering
Science,1978.33:p.303314.
29. Ebrahimi, A.A., H.A. Ebrahim, and E. Jamshidi, Solving partial differential
equations of gassolid reactions by orthogonal collocation. Computers and
ChemicalEngineering,2008.32:p.17461759.
30. Lim,Y.I.,J.M.LeLann,andX.Joulia,Accuracy,temporalperformanceandstability
comparisons of discretization methods for the numerical solution of Partial
Differential Equations (PDEs) in the presence of steep moving fronts. Computers
andChemicalEngineering,2001.25:p.14831492.
31. Patankar,S.,Numericalheattransferandfluidflow.2980,NewYork:Hemisphere
Publishing.
32. Botte, G.G., J.A. Ritter, and R.E. White, Comparisonoffinitedifferenceandcontrol
volumemethodsforsolvingdifferentialequations. Chemical Engineering Science,
2000.24:p.26332654.
33. Vinokur, M., An Analysis of FiniteDifference and FiniteVolume Formulations of
ConservationLaws.JournalofComputationalPhysics,1989.81:p.152.
34. Stynes, M., Finite volume methods for convectiondiffusion problems. Journal of
ComputationalandAppliedMathematics,2995.63:p.8390.
35. LeVeque, R.J., Finite volume methods for hyperbolic problems. 2002, New York:
CambridgeUniversityPress.
36. Leonard,B.P.,Comparisonoftruncationerroroffinitedifferenceandfinitevolume
formulationofconvectionterms.AppliedMathematicsModelling,1994.18:p.46
50.
37. Rouboa,A.andE.Monteiro,Heattransferinmultiblockgridduringsolidification
Performance of Finite Differences and Finite Volume methods. Journal of
MaterialsProcessingTechnology,2008.204:p.451458.
38. Hoffmann, K.A. and S.T. Chiang, Finite Volume Method, in Computational Fluid
Dynamics.2000,EngineeringEducationSystem:Wichita.p.385417.
39. Langmuir,I.,J.Amer.Chem.Soc,1918.40:p.1361.
40. HerrmannGeppert,I.andU.Kramm.Gassorptionmeasurements.21stSep2010
[cited 2011 4th Jun]; Available from: http://www.helmholtz

66
 CN4118RFinalReport

berlin.de/forschung/enma/solarebrennstoffe/analytische
methoden/gassorptionsmessungen_en.html.
41. Glueckauf, E., K.H. Barker, and G.P. Kitt, Theory of chromatography VIII. The
separation of lithium isotopes by ion exchange and of neon isotopes by low
temperatureadsorptioncolumns Discussion of the Faraday Society, 1949. 7: p.
199213.
42. Bischoff,K.B.andO.Levenspiel,Fluiddispersiongeneralizationandcomparison
of mathematical models. Chemical Engineering Science, 1962. 17: p. 245255,
257264.
43. Warming, R.F. and R.M. Beam. Upwind secondorder difference schemes and
applications in unsteady aerodynamic flows. in Proc. AIAA 2nd Computational
FluidDynamicsConf.1975.Hartford,Conn.
44. Lax, P. and B. Wendroff, Systemsofconservationlaws. Communications on Pure
andAppliedMathematics,1960.13(2):p.217237.
45. Crank, J. and P. Nicolson, Apracticalmethodfornumericalevaluationofsolution
of partial differential equations of the heatconduction type. Advances in
ComputationalMathematics,1996.6(1):p.207226.
46. Ogata, A. and R.B. Banks, A solution of the differential equation of longitudinal
dispersioninporousmedia.USGeologicalSurveyProfessionalPaper411A,1961.
47. Klinkenberg, A., Numericalevaluationofequationsdescribingtransientheatand
masstransferinpackedsolids.Ind.Eng.Chem.,1948.40(10):p.19921994.
48. MathWorks. ODE Solvers. 2011 [cited 2011 6th Jun]; Available from:
http://www.mathworks.com/help/techdoc/math/f1662913.html#f1753130.
49. Tyson, R., L.G. Stern, and R.J. LeVeque, Fractional step methods applied to a
chemotaxismodel.JournalofMathematicalBiology,2000.41:p.455475.

67